Patch for bug in agrad row() and col()
Marcus Brubaker
(Last email for the night, honest!)
Cheers,
Marcus
Bob Carpenter
My bad -- I just fixed the math/matrix row() and col()
and didn't fix the agrad ones.
We also need to keep the unit tests in synch with the
code, but I'll take care for that for the row/col
updates.
The problem is that we're really thinking of the C++ lib
as having 0-based indexing. The Stan code generator
does all the subtraction of 1 to get the MATLAB/BUGS/R-style
indexing from 1. Except for the row() and col() functions.
So what we really need is some other row() and col() functions.
So I'dd add row0() and col0() methods if we really wanted
to fully support 0-based indexing at the API level.
On our previous topic, efficiency concerns are different
in the agrad case than the double case. Anything that
saves creating new agrad::var instances is going to be a win
over using the clever/small template expressions. This will
mainly be an issue in the (Matrix<double> * var) methods,
because the promotion of the Matrix<double> to Matrix<var> is
so expensive. It's faster to just create the Matrix<var>
result directly.
- Bob
Marcus Brubaker
Thanks for spotting that.
My bad -- I just fixed the math/matrix row() and col()
and didn't fix the agrad ones.
We also need to keep the unit tests in synch with the
code, but I'll take care for that for the row/col
updates.
The problem is that we're really thinking of the C++ lib
as having 0-based indexing. The Stan code generator
does all the subtraction of 1 to get the MATLAB/BUGS/R-style
indexing from 1. Except for the row() and col() functions.
So what we really need is some other row() and col() functions.
So I'dd add row0() and col0() methods if we really wanted
to fully support 0-based indexing at the API level.
Yeah, I see that inconsistency may be problematic. Is it worth trying to shift so that only the modeling language uses 1-based indexing and just put a bunch of "-1"'s when outputing the code? This would make the library a lot more consistent for people who use the C++ directly. Alternatively, maybe put all the 1-based indexing functions in a separate file/namespace?
On our previous topic, efficiency concerns are different
in the agrad case than the double case. Anything that
saves creating new agrad::var instances is going to be a win
over using the clever/small template expressions. This will
mainly be an issue in the (Matrix<double> * var) methods,
because the promotion of the Matrix<double> to Matrix<var> is
so expensive. It's faster to just create the Matrix<var>
result directly.
Yeah, I'm beginning to see that. I did a little profiling of my model yesterday and found that about 60% of the time is spent in agrad::chainable::operator new() and the agrad::vari(double) constructor. I was able to cut that a bit by using a Cholesky decomposition in multi_normal_log (I'll send along this stuff a little later once it's cleaned up) but it still constitutes a huge chunk of computation.
Cheers,
Marcus
Ben Goodrich
Yeah, I'm beginning to see that. I did a little profiling of my model yesterday and found that about 60% of the time is spent in agrad::chainable::operator new() and the agrad::vari(double) constructor. I was able to cut that a bit by using a Cholesky decomposition in multi_normal_log (I'll send along this stuff a little later once it's cleaned up) but it still constitutes a huge chunk of computation.
The textbookish multi_normal_log() function is a self-inflicted efficiency disaster, and we had planned to change it. I actually had a branch to fix this when I was using git-svn and stan was using svn on the server. But it was breaking the build for some reason at the time so I couldn't push it and then it got disrupted by the movement of files.
The basic problem is that check_cov_matrix() does a Cholesky decomposition to check whether Sigma is positive definite, but that local computation gets lost when check_positive_definite() exits. So, the L matrix is not available to calculate the determinant or the inverse in the multinormal density. The cherry on top is that we already construct L in cov_matrix_constrain() and then post-multiply it by L' to get Sigma (to then decompose three times).
It is better to utilize multi_normal_cholesky_log() if possible.
Ben
Marcus Brubaker
You can also more efficiently (and stably) compute the log determinant from L. Perhaps the way to write it is to compute L first and then just use that to support the rest of the needed checks.
Perhaps even better, compute L and call multi_normal_cholesky_log(). Thoughts?
It is better to utilize multi_normal_cholesky_log() if possible.
It seems like just adding a chol() function would probably solve most of these problems. I'll look into that.
Cheers,
Marcus
Bob Carpenter
The only places it touches the API is in row() and col().
I was originally thinking I'd need to do something to declare
which arguments of which functions needed a -1 in them,
but it was simpler to just change row() and col() and put
in warnings in the doc. If we get lots of people wanting to
use our API, we may need to rethink.
Having multiple namespaces is trickier. I'd rather try to
keep that simpler. The old-school approach would IFDEF the
-1 with some property. The new-school approach would be to
add policy template args to everything. This all seems like
overkill.
The right solution in the longer term is to get rid of the
matrix library altogether in favor of adding
heterogeneous binary ops to Eigen itself (allowing
you to multiply Matrix<A> with Matrix<B> whenever you can
multiply an A by a B, in particular with A=double, B=stan::agrad::var),
we can get rid of our whole matrix lib, which is just there
to wrap stan::agrad::var so that we can use matrix ops.
No matter what we do the gradients are going to be the
biggest part of the computation. We've got our own arena-based
memory management for the new() calls. There's a ton of work
being done in the virtual chain methods, but it's either that
or pointers to functions, and either way, it's a deref for
every chaining function call. There's really no way to do it
all statically -- we need to store different kinds of expressions
on the derivative propagation stack.
Ben's got a bunch of stuff queued up for Cholesky decompositions
for multivariates, so you should definitely run what you want to
do by us -- it might already be done (at least conceptually).
I've been trying to control the feature creep so we can actually
get a working version out the door.
- Bob
Ben Goodrich
You can also more efficiently (and stably) compute the log determinant from L. Perhaps the way to write it is to compute L first and then just use that to support the rest of the needed checks.
Perhaps even better, compute L and call multi_normal_cholesky_log(). Thoughts?
That is exactly what I had done in my branch that got clobbered. I can redo it soon.
The bigger and more invasive thing is avoiding 1 and 2 in:
0. Construct L from the unbounded primitive parameters
1. Set Sigma := L * L'
2. Calculate L := Sigma.llt().matrixL()
3. Calculate the determinant and inverse of Sigma from L
Ben
Bob Carpenter
(Maybe my naming was too verbose? Now would be the
time to change it.)
But it only returns the decomposition, not the entire
LLT solver object. And it's not doing error checking
itself -- so that should probably be added.
In general, it's better to do the matrix ops in Eigen
directly where possible, because our library's explicitly
clobbering the expression templates in order to deal
properly with the auto-dif variables in mixed
var/double operations.
- Bob
On 3/9/12 4:06 PM, Marcus Brubaker wrote:
...
Marcus Brubaker
On Friday, March 9, 2012 4:06:36 PM UTC-5, Marcus Brubaker wrote:You can also more efficiently (and stably) compute the log determinant from L. Perhaps the way to write it is to compute L first and then just use that to support the rest of the needed checks.
Perhaps even better, compute L and call multi_normal_cholesky_log(). Thoughts?
That is exactly what I had done in my branch that got clobbered. I can redo it soon.
Actually, I just implemented something basically like that which seems to be working. I need to test it a bit more, but it's pretty close.
The bigger and more invasive thing is avoiding 1 and 2 in:
0. Construct L from the unbounded primitive parameters
1. Set Sigma := L * L'
2. Calculate L := Sigma.llt().matrixL()
3. Calculate the determinant and inverse of Sigma from L
Where are 1 and 2 being done exactly? I assume that happens only when Sigma is a parameter?
Cheers,
Marcus
Marcus Brubaker
There is a cholesky_decompose() method in matrix.hpp.
(Maybe my naming was too verbose? Now would be the
time to change it.)
Ah, I had missed that. The naming is just a matter of preference, I skew towards shorter names where they make sense but my fingers are just lazy.
But it only returns the decomposition, not the entire
LLT solver object. And it's not doing error checking
itself -- so that should probably be added.
In general, it's better to do the matrix ops in Eigen
directly where possible, because our library's explicitly
clobbering the expression templates in order to deal
properly with the auto-dif variables in mixed
var/double operations.
I ended up implementing it that way. I'll commit it once it is better tested.
Cheers,
Marcus
Bob Carpenter
> I ended up implementing it that way. I'll commit it once it
> is better tested.
Did you figure out the unit testing scheme?
That's where to put the tests.
We haven't been testing behavior of auto-dif,
but we are trying to test that the distributions
actually work with all the possible combinations
of double and agrad vars.
- Bob
Ben Goodrich
On Fri, Mar 9, 2012 at 5:39 PM, Ben Goodrich wrote:On Friday, March 9, 2012 4:06:36 PM UTC-5, Marcus Brubaker wrote:You can also more efficiently (and stably) compute the log determinant from L. Perhaps the way to write it is to compute L first and then just use that to support the rest of the needed checks.
Perhaps even better, compute L and call multi_normal_cholesky_log(). Thoughts?
That is exactly what I had done in my branch that got clobbered. I can redo it soon.
Actually, I just implemented something basically like that which seems to be working. I need to test it a bit more, but it's pretty close.
Good
The bigger and more invasive thing is avoiding 1 and 2 in:
0. Construct L from the unbounded primitive parameters
1. Set Sigma := L * L'
2. Calculate L := Sigma.llt().matrixL()
3. Calculate the determinant and inverse of Sigma from L
Where are 1 and 2 being done exactly? I assume that happens only when Sigma is a parameter?
Yes, 1 and 2 only apply when Sigma is a parameter. If Sigma is fixed data, we had emailed about inverting Sigma once at the outset (and optionally calculating the determinant) and passing the precision matrix to an appropriately named version of the multivariate normal log density parameterized in terms of the precision matrix. However, nothing came of that discussion yet.
The sequence of calls is in src/stan/prob/transform.hpp: cov_matrix_constrain() does some transformations calls read_cov_matrix(), which calls read_cov_L() to produce L and then returns L * L' . What is needed is a function possibly named L_constrain() that does some transformation and calls read_cov_L() to return L. Then, that L should get passed to mult_normal_cholesky_log() or any other elliptical distribution we may have in the future. Furthermore, if L_constrain() existed, then we could make cov_matrix_constrain() just call L_constrain() and return L * L'.
If you want to do that, go ahead. Otherwise, I'll get to it soon.
Ben
Ben Goodrich
The sequence of calls is in src/stan/prob/transform.hpp: cov_matrix_constrain() does some transformations calls read_cov_matrix()
Bob, why is the Jacobian adjustment commented out in:
template <typename T>
Matrix<T,Dynamic,Dynamic> cov_matrix_constrain(const Matrix<T,Dynamic,1>& x,
size_t k,
T& lp) {
size_t k_choose_2 = (k * (k - 1)) / 2;
Array<T,Dynamic,1> cpcs(k_choose_2);
int pos = 0;
for (size_t i = 0; i < k_choose_2; ++i)
cpcs[i] = corr_constrain(x[pos++]); // ,lp);
Array<T,Dynamic,1> sds(k);
for (size_t i = 0; i < k; ++i)
sds[i] = positive_constrain(x[pos++]); // ,lp);
return read_cov_matrix(cpcs, sds, lp);
}
Also, the last sections in transform.hpp have a lot of non-idiomatic Eigen. I think we can do something like
template <typename T>
Matrix<T,Dynamic,Dynamic> cov_matrix_constrain(const Matrix<T,Dynamic,1>& x,
size_t k,
T& lp) {
size_t k_choose_2 = (k * (k - 1)) / 2;
Array<T,Dynamic,1> cpcs = x.head(k_choose_2).array().tanh();
Array<T,Dynamic,1> sds = x.tail(k).array().exp();
return read_cov_matrix(cpcs, sds, lp);
}
Which brings me to the question of why we are using std vectors instead of Eigen Arrays to hold (parts of) the parameter vector? It seems like we do a lot of verbose relocations of things in one "vector" to another "vector".
Ben
Marcus Brubaker
But wouldn't that require special logic to determine if the covariance matrix is being passed to an elliptical distribution and then passing L instead to a different function.
Maybe a simpler solution is to introduce a cholesky_factor type which can be used explicitly in model building?
Cheers,
Marcus
Ben Goodrich
But wouldn't that require special logic to determine if the covariance matrix is being passed to an elliptical distribution and then passing L instead to a different function.
Maybe a simpler solution is to introduce a cholesky_factor type which can be used explicitly in model building?
This is more or less where the discussion ran out of steam last time. Part of the problem is that Eigen does not have a triangular matrix type or a symmetric matrix type, only triangular and self-adjoint "views" of a full matrix, but Doug Bates among other people is hoping that these get added soon. Bob is somewhat reluctant to add new types to stan. I think Andrew was somewhat reluctant to have default behavior that deviates from what people see in textbooks. I kind of liked the idea of having a covariance matrix class whose attributes included all the stuff that any function might need, such as its Cholesky factor, its inverse, a bunch of flags for validity checks, etc. But that would be inconsistent with Stan's mostly functional orientation. So that will have to get worked out. As it stands now, Stan users will likely gravitate toward multivariate models that are at best slow and at worst wrong.
Ben
Bob Carpenter
done so we'd get the right answers. I had reasoned
through them to the point where I thought you were
applying the Jacobians yourself, so I took them out
of the transform arguments.
There may still be mistakes, there. If so, could
you add unit tests that break before adding the
fixes?
Also, please feel free to refactor all of the
Eigen linear algebra. I repeat, there's no reason
for my textbook implementations other than ignorance!
And if we can add more tests for the derivatives of
some of these complicated things, even better.
The reason there are std::vector instances all over
is that I was trying to minimize the dependence
on Eigen. In retrospect, that was probably a mistake.
The copying is also the least of our problems
performance-wise in the context of Stan. But, if we
want to provide a good external library, then we
have to clean this stuff up.
An interim solution to ease the refactoring would be
to add overloaded versions of the functions for Eigen
vectors and matrices.
For vectors, we can lean on the fact that Eigen column
vectors, Eigen row vectors, and std::vector all
implement the .size(), and operator[](size_t i) concepts
(though actually, Eigen uses operator[](int) and
thus goes against the STL conventions). Given that,
we could use some fancier templating and
treat the argument/result as a template template (because
we need the item type, too, so need to pull out type
params for both std::vector<T> and T)
Matt's suggestion was to take some of these functions
even lower level with iterator/pointer arguments instead
of array arguments. That makes more sense if you don't
have to do Eigen operations on them.
The multivariate parameters for densities
are all Eigen vectors and matrices, which is where the
incompatibility comes in.
Within an expression in Eigen, you can efficiently Map
a std::vector<T> as an Eigen vector,
but you can't wrap std::vector<std::vector<T>>.
If we want to keep going on optimizing, the way I implemented
the compiled Stan data structures is inefficient as it
relies on std::vector<std::vector<double>> if you declare
double x[M,N]
That's wasteful in storage compared to an Eigen Matrix, which just
computes the indexing itself (using the row/col major policy).
But we'd need to generalize to x[M,N,K], and so on, to handle
all of the Stan types.
- Bob
Bob Carpenter
to get a stable version of Stan out the door. If you
look at what's involved in adding a new type, it's fairly
involved, but not impossible.
I like the idea of a covariance matrix class. Much of
what we were talking about, like caching the Cholesky
factor, could be handled behind the scenes. As I said before,
it would be more like exposing Eigen's LLT object.
As to speed, I want to emphasize that the best thing
we can do for speed is to partially evaluate all the
gradients in:
stan/agrad/matrix.hpp
The time for doing these ops dominates everything else,
and everywhere we promote a double to a var, we take
a HUGE hit in performance. Every time you evaluate
x * y
and x is a double matrix and y is a var matrix, we
create an entire var matrix needlessly, leading not
just to M x N copies, var constructor calls, but also
the overhead of access when we actually do the multiplication,
which itself spews more copies. This goes for vectors,
too.
I did it that way not because I didn't know how to do it
better, but because it would take a long time to implement
and test it all done the right way.
The poor choices in linear algebra operators were done out
of ignorance, and I again urge you to fix them. Or we
can sit down and fix them together.
Meanwhile, we need to write that R interface before we
can make a Stan user out of Andrew.
- Bob
:-)
>
> Ben
>
Bob Carpenter
of the ops. When you have a covariance matrix
that's being used for multiple observations,
the hard work only needs to be done once.
It's easy enough to do with simple vector data outcomes,
but harder to do with regression models. Either
the linear (or whatever) predictor needs to be
calculated into an array before calling a vectorized function,
or you need to do the implementation yourself using
vector/matrix ops (presumably what Ben meant by saying
users would do the wrong thing -- he's right, because
users do all sorts of crazy things, like my very
naive use of Eigen).
- Bob
Ben Goodrich
I fiddled with where the Jacobians were getting
done so we'd get the right answers. I had reasoned
through them to the point where I thought you were
applying the Jacobians yourself, so I took them out
of the transform arguments.There may still be mistakes, there. If so, could
you add unit tests that break before adding the
fixes?
If our previous agreement that "the function that performs the transformation is responsible for adjusting lp" is still operative, then cov_matrix_constrain() should be responsible for adjusting lp due to
A. mapping of unbounded parameters to canonical partial correlations on (-1,1)
B. mapping of unbounded parameters to positive standard deviations
Then, downstream functions would be responsible for adjusting lp due to
C. mapping of canonical partial correlations to Pearson correlations
D. mapping of Pearson correlations to covariances
A month or so ago, there was a bug due to me being off-by-one in part D, but I think you commented out correct implementations of parts A and B. In any event, inconsistency (with the same inits) between wishart and wishart2x2 in basic_distributions is sufficient reason to be concerned, right?
master!stan *> ./demo/wishart --seed=1331410622 --test_grad
STAN SAMPLING COMMAND
data =
init = random initialization
samples = samples.csv
append_samples = 0
seed = 1331410622 (user specified)
chain_id=1 (user specified)
iter = 2000
warmup = 1000
thin = 1
leapfrog_steps = -1
epsilon = -1
epsilon_pm = 0
epsilon_adapt_off = 0
delta = 0.5
gamma = 0.05
TEST GRADIENT MODE
Log probability=-6.64661
param idx value model finite diff error
0 0.589153 -0.529286 -0.529287 3.59687e-07
1 -1.45078 3.97253 3.97253 2.7128e-08
2 0.228394 0.842013 0.84201 3.15721e-06
master!stan *> ./demo/wishart2x2 --seed=1331410622 --test_grad
STAN SAMPLING COMMAND
data =
init = random initialization
samples = samples.csv
append_samples = 0
seed = 1331410622 (user specified)
chain_id=1 (user specified)
iter = 2000
warmup = 1000
thin = 1
leapfrog_steps = -1
epsilon = -1
epsilon_pm = 0
epsilon_adapt_off = 0
delta = 0.5
gamma = 0.05
TEST GRADIENT MODE
Log probability=-4.58141
param idx value model finite diff error
0 0.589153 -0.429512 -0.429513 3.44415e-07
1 -1.45078 1.9414 1.9414 2.98465e-08
2 0.228394 0.74342 0.743419 6.2906e-07
master!stan *>
Ben
Ben Goodrich
Meanwhile, we need to write that R interface before we
can make a Stan user out of Andrew.
Or do you still think we should consider having a 1.0 version that dumps the data and calls stan from the command line via system() ?
Ben
Andrew Gelman
But I say this in ignorance of all details of implementation.
A
Ben Goodrich
The Wishart-family approach in stan mostly makes sense from a syntax perspective, e.g.:
master!stan *> cat src/models/basic_distributions/wishart.stan
transformed data {
cov_matrix(2) S;
for (i in 1:2)
for (j in 1:2)
S[i,j] <- 0.0;
S[1,1] <- 2.0;
S[2,2] <- 0.5;
}
parameters {
cov_matrix(2) W;
}
model {
W ~ wishart(4, S);
}
The user has a Wishart prior and the Wishart prior implies some unnamed but calculatable prior beliefs about the primitive parameters that are transformed behind the scenes in cov_matrix(). That was a bit weird to program in C++, but is not wrong and not exposed to the user anyway.
But I anticipate users will tend to utilize Wishart-family priors on covariance matrices because that is what they are familiar with in the literature, rather than that is what is entailed by their prior beliefs. The LKJ approach allows the user to have a prior over each standard deviation and a joint prior (usually uniform, which is implied when eta == 1.0) over the correlation matrix.
transformed data {
real mu(2);
real(0,) sigma(2);
real(0,) eta;
mu[1] <- 0.0;
mu[2] <- 0.0;
sigma[1] <- 1.0;
sigma[2] <- 0.5;
eta <- 1.0;
}
parameters {
cov_matrix(2) W;
}
model {
W ~ lkj_cov(mu, sigma, eta);
}
But lkj_cov_log() is hard-coded to impose a lognormal prior on the standard deviations, which seems unnecessarily limiting. And there is the aforementioned efficiency nightmare where cov_matrix_constrain() first constructs L so that it can be post-multiplied by L' to get W and then lkj_cov_log() conducts a Cholesky decomposition of W to (re)obtain L and then evaluate the density.
To do things the right way, we need a function, perhaps called LKJ_mapping(), that is exposed to the user that takes CPCs and standard deviations and produces L such that W := L * L'. And we need a (trivial) C++ version of lkj_cor_log() that takes CPCs and eta. Then, the .stan file could look something like
parameters {
real(-1,1) CPCs[K_choose_2];
real(0,) SDs[K];
}
transformed parameters {
matrix L(K,K);
L <- LKJ_mapping(CPCs, SDs);
}
model {
for(i in 1:K)
// some prior on SDs[i]
CPCs ~ lkj_cor(1.0);
// maybe use L as part of a multinormal likelihood or something
}
Ben
Ben Goodrich
If we want to keep going on optimizing, the way I implemented
the compiled Stan data structures is inefficient as it
relies on std::vector<std::vector<double>> if you declaredouble x[M,N]
That's wasteful in storage compared to an Eigen Matrix, which just
computes the indexing itself (using the row/col major policy).
But we'd need to generalize to x[M,N,K], and so on, to handle
all of the Stan types.
Array< MatrixXd, Dynamic, 1 > x;
is supposed to be a valid declaration.
Ben
Bob Carpenter
> Another limitation of the way covariance matrices are handled is that we
> currently do not have a way for the user to utilize the LKJ approach in
> an intuitive and computationally efficient way. Granted, there are
> currently only 4 people in the world who think in terms of the LKJ
> approach and three of them are L, K, and J. But that is the best way to
> think about it in the stan context.
If you can convince Andrew of this and get him to
use his bully pulpit AND we put it in Stan, then I
think that number will be bigger than 4 before long.
From a user perspective, the LKJ version is pretty
easy to grasp. Easier than a scaled Wishart, at least.
> The Wishart-family approach in stan mostly makes sense from a syntax
> perspective, e.g.:
Right. Currently, it looks just like BUGS.
Except that we declare W as a covariance matrix.
> master!stan *> cat src/models/basic_distributions/wishart.stan
> transformed data {
> cov_matrix(2) S;
>
> for (i in 1:2)
> for (j in 1:2)
> S[i,j] <- 0.0;
>
> S[1,1] <- 2.0;
> S[2,2] <- 0.5;
A symmetric matrix data type would be nice.
Or a diagonal matrix data type for this case.
Also, you shouldn't be using a cov_matrix(2) data
type in transformed data. The type is only guaranteed
for reading in the parameter and not monkeying with it.
So I turned it off in this context. This is probably
going to confuse people.
> }
> parameters {
> cov_matrix(2) W;
> }
> model {
> W ~ wishart(4, S);
> }
>
> The user has a Wishart prior and the Wishart prior implies some unnamed
> but calculatable prior beliefs about the primitive parameters that are
> transformed behind the scenes in cov_matrix(). That was a bit weird to
> program in C++, but is not wrong and not exposed to the user anyway.
That was the plan! Looks like there may be a bug or two,
but conceptually it should work.
> But I anticipate users will tend to utilize Wishart-family priors on
> covariance matrices because that is what they are familiar with in the
> literature, rather than that is what is entailed by their prior beliefs.
Yes indeed. Andrew's very worried about this kind
of thing, too. I don't know if he's onboard with the LKJ
priors, but convincing him and building this into Stan could
go a long way toward convincing everyone else.
Now's the time to bend his ear about inclusion of material
in BDA 3, too :-)
> The LKJ approach allows the user to have a prior over each standard
> deviation and a joint prior (usually uniform, which is implied when eta
> == 1.0) over the correlation matrix.
>
> transformed data {
> real mu(2);
> real(0,) sigma(2);
> real(0,) eta;
>
> mu[1] <- 0.0;
> mu[2] <- 0.0;
> sigma[1] <- 1.0;
> sigma[2] <- 0.5;
> eta <- 1.0;
> }
> parameters {
> cov_matrix(2) W;
> }
> model {
> W ~ lkj_cov(mu, sigma, eta);
> }
>
> But lkj_cov_log() is hard-coded to impose a lognormal prior on the
> standard deviations, which seems unnecessarily limiting.
Agreed.
> And there is
> the aforementioned efficiency nightmare where cov_matrix_constrain()
> first constructs L so that it can be post-multiplied by L' to get W and
> then lkj_cov_log() conducts a Cholesky decomposition of W to (re)obtain
> L and then evaluate the density.
>
> To do things the right way, we need a function, perhaps called
> LKJ_mapping(), that is exposed to the user that takes CPCs and standard
> deviations and produces L such that W := L * L'. And we need a (trivial)
> C++ version of lkj_cor_log() that takes CPCs and eta. Then, the .stan
> file could look something like
>
> parameters {
> real(-1,1) CPCs[K_choose_2];
> real(0,) SDs[K];
> }
> transformed parameters {
> matrix L(K,K);
> L <- LKJ_mapping(CPCs, SDs);
> }
> model {
> for(i in 1:K)
> // some prior on SDs[i]
>
> CPCs ~ lkj_cor(1.0);
>
> // maybe use L as part of a multinormal likelihood or something
> }
That looks great. We don't even need to compute anything
for lkj_cor(alpha) when alpha = 1.0.
All that we need to do this in Stan are the four instances
of:
template <typename T1, typename T2>
Matrix<typename promote_args<T1,T2>::type,Dynamic,Dynamic>
LKJ_mapping(std::vector<T1> cpcs, std::vector<T2> sds);
and
typename promote_args<T1,T2>::type
lkj_cor(std::vector<T1> cpcs, T2 alpha);
That's for T1 in {double,var} and T2 in {double,var}.
In order to use Eigen vectors instead of std::vector in
the function, we'd need two new Stan types and declarations
like:
parameters {
corr_vector(K_choose_2) CPCs;
pos_vector(K) sds;
}
Or I could try to thread through constraints on vector
values. Right now, only non-vectors can have constraints.
- Bob
Bob Carpenter
On 3/10/12 5:29 PM, Ben Goodrich wrote:
> On Saturday, March 10, 2012 2:51:20 PM UTC-5, Bob Carpenter wrote:
>
> If we want to keep going on optimizing, the way I implemented
> the compiled Stan data structures is inefficient as it
> relies on std::vector<std::vector<double>> if you declare
>
> double x[M,N]
>
> That's wasteful in storage compared to an Eigen Matrix, which just
> computes the indexing itself (using the row/col major policy).
> But we'd need to generalize to x[M,N,K], and so on, to handle
> all of the Stan types.
>
> Do we just need three dimensional or do we need arbitrary dimensional?
Stan allows as many dimensions as you want.
> If it is the former, note that Eigen supports 3-dimensional Arrays in a
> kludgy fashion by allowing a Matrix or an(other) Array to be the "type"
> for an Array. In other words,
>
> Array< MatrixXd, Dynamic, 1 > x;
>
> is supposed to be a valid declaration.
Sure, but that's going to wash out just like
vector<vector<>> in terms of size (assuming they haven't
built specializations of Array for this case).
You could probably keep going and also have
Array<Array<Array...> but unless they have lots of
template specializations, it'll wash out just like
std::vector<std::vector...> with lots of substructures
redundantly storing length information and pointers
and managing memory.
If you stack them, you can probably use some of the
arithmetic on arrays or matrices, too, because the
template type implements arithmetic.
- Bob
Bob Carpenter
On 3/10/12 3:34 PM, Ben Goodrich wrote:
> On Saturday, March 10, 2012 2:51:20 PM UTC-5, Bob Carpenter wrote:
>
> I fiddled with where the Jacobians were getting
> done so we'd get the right answers. I had reasoned
> through them to the point where I thought you were
> applying the Jacobians yourself, so I took them out
> of the transform arguments.
>
> There may still be mistakes, there. If so, could
> you add unit tests that break before adding the
> fixes?
>
>
> If our previous agreement that "the function that performs the
> transformation is responsible for adjusting lp" is still operative,
Yes, the function
matrix cov_matrix_constrain(vector, lp);
needs to convert an ((N choose 2) + N)-dimensional
vector input, with unconstrained finite floating-point values,
into a positive-definite covariance matrix output and
add the log of the absolute Jacobian determinant of the
transform to lp. (Presumably there's a measure-zero subspace
of input vectors with non-positive-definite outputs.)
> then
> cov_matrix_constrain() should be responsible for adjusting lp due to
>
> A. mapping of unbounded parameters to canonical partial correlations on
> (-1,1)
> B. mapping of unbounded parameters to positive standard deviations
> Then, downstream functions would be responsible for adjusting lp due to
>
> C. mapping of canonical partial correlations to Pearson correlations
> D. mapping of Pearson correlations to covariances
template <typename T>
Matrix<T,Dynamic,Dynamic>
cov_matrix_constrain(const Matrix<T,Dynamic,1>& x,
size_t k,
T& lp) {
size_t k_choose_2 = (k * (k - 1)) / 2;
Array<T,Dynamic,1> cpcs(k_choose_2);
int pos = 0;
for (size_t i = 0; i < k_choose_2; ++i)
[1] cpcs[i] = corr_constrain(x[pos++]); // ,lp);
Array<T,Dynamic,1> sds(k);
for (size_t i = 0; i < k; ++i)
[2] sds[i] = positive_constrain(x[pos++]); // ,lp);
return read_cov_matrix(cpcs, sds, lp);
}
So if read_cov_matrix() isn't applying the Jacobians for
steps (A) and (B), we need to add the lp back to [1] and [2].
And read_cov_matrix() is going to be responsible for C and D.
> A month or so ago, there was a bug due to me being off-by-one in part D,
> but I think you commented out correct implementations of parts A and B.
I very well could have.
> In any event, inconsistency (with the same inits) between wishart and
> wishart2x2 in basic_distributions is sufficient reason to be concerned,
> right?
Indeed. Something's wrong somewhere if they're not producing
the same results. The param that gets printed is the unconstrained
version. So it's not surprising those are the same. But the
log prob and gradients should be the same, too.
Missing or duplicating part of the Jacobian is a likely culprit.
- Bob
Andrew Gelman
A