Check my work for Cholesky covar transform?

[Bob Carpenter]

Could someone check my work on the
Cholesky-factor covariance transform?

I don't need the code reviewed, just the
math, which is in the doc, the relevant
two pages of which I've attached.

- Bob

[Ben Goodrich]

It looks fine to me. I would continue to prefer the phrase "Covariance Matrix Inverse Transform" to your "Inverse Covariance Matrix Transform" because the latter is confusing if the reader mentally puts the parentheses like (Inverse Covariance Matrix) Transform instead of (Inverse) (Covariance Matrix Transform).

[Ben Goodrich]

Also, it might be good to mention in the documentation that a precision matrix can be composed in exactly the same way. This is better computationally because it avoids the need to do any matrix inversions for most models involving elliptical distributions. But, at the moment, users would have to do everything themselves in the .stan file because the C++ code does not have the necessary parameterizations for this.

Ben

On Tuesday, April 24, 2012 2:58:20 PM UTC-4, Bob Carpenter wrote:

[Bob Carpenter]

I'm indifferent, so I'll change it. And change
all the other ones so they're consistent.

Noun compounds are always problematic in terms
of both structure (where the implicit parens go)
and of meaning (what the relation between the nouns is).

I should really say "Inverse of the Covariance Matrix Transform"
or even "Inverse of the Transform of Covariance Matrices"
or even "Inverse of the Transform from Covariance Matrices to
Unbounded Column Vectors"! Maybe I should just write the
formula? Just kidding.

- Bob

[Bob Carpenter]

Thanks for looking at the math. I changed the
text to make it consistent everywhere.

I just realized I can use the Jacobian calculation
from Stan's auto-dif to check my work. I'll add that
test later tonight.

On 4/24/12 3:22 PM, Ben Goodrich wrote:
> Also, it might be good to mention in the documentation that a precision
> matrix can be composed in exactly the same way.

I agree. Should I change the data type form
cov_matrix to pos_def_matrix? I wasn't even sure
these were the same concept when I added the first
time (though we're really looking at only symmetric
pos def matrices).

Also, we can mention that users can rescale
a correlation matrix. I'll eventually replace
that transform with the hyperspherical transform
if it looks like it'll be more efficient.

> This is better
> computationally because it avoids the need to do any matrix inversions
> for most models involving elliptical distributions. But, at the moment,
> users would have to do everything themselves in the .stan file because
> the C++ code does not have the necessary parameterizations for this.

So maybe it can wait.

Or feel free to add the parameterizations.
I think the dust has settled enough around the
transforms and error checking that it's safe
to add again.

- Bob

[Ben Goodrich]

On Tuesday, April 24, 2012 3:38:39 PM UTC-4, Bob Carpenter wrote:

> Also, it might be good to mention in the documentation that a precision

> matrix can be composed in exactly the same way.

I agree.  Should I change the data type form
cov_matrix to pos_def_matrix?  I wasn't even sure
these were the same concept when I added the first
time (though we're really looking at only symmetric
pos def matrices).

Although the asymmetric positive definite matrices almost never come up in statistics, it is probably not good to have a data type of pos_def_matrix that implicitly also requires symmetry. We could have sym_pos_def_matrix but that is a bit long-winded. I'd be happy to stick with cov_matrix in the .stan file if we had a cov_matrix class in the C++ file that could represent a covariance matrix in a myriad of different ways. Then we could parameterize it via the Cholesky factor of the precision matrix without the user really knowing.

Also, we can mention that users can rescale
a correlation matrix.  I'll eventually replace
that transform with the hyperspherical transform
if it looks like it'll be more efficient.

I worked out that the LKJ and the onion are particular hyperspherical transforms. It is just that they are the only hyperspherical transforms that imply a known joint density for the correlation matrix. You could use another density on the hypersphere(s) and get some correlation matrix, but its joint density would just be (proportional to) the product of the (independent) hypersphere densities times a Jacobian determinant. Hard to see what would be the use of implementing that in stan until someone does it in such a way that the joint density becomes analytically attractive.

Ben

[Bob Carpenter]

On 4/24/12 4:50 PM, Ben Goodrich wrote:
> On Tuesday, April 24, 2012 3:38:39 PM UTC-4, Bob Carpenter wrote:
>
> > Also, it might be good to mention in the documentation that a
> precision
>
> > matrix can be composed in exactly the same way.
>
> I agree. Should I change the data type form
> cov_matrix to pos_def_matrix? I wasn't even sure
> these were the same concept when I added the first
> time (though we're really looking at only symmetric
> pos def matrices).
>
> Although the asymmetric positive definite matrices almost never come up
> in statistics, it is probably not good to have a data type of
> pos_def_matrix that implicitly also requires symmetry.

Because they tend to get used in the (y'Ay) setting, everything
I've been reading starts with "without loss of generality,
assume the pos def matrix is symmetric".

By the way, Matt found this book:

http://cs.nenu.edu.cn/Lunwen/2010050513555973499.pdf

on matrix differential calc -- it has all sorts of
great material on Jacobians and derivatives of
matrix operations.

> We could have
> sym_pos_def_matrix but that is a bit long-winded. I'd be happy to stick
> with cov_matrix in the .stan file if we had a cov_matrix class in the
> C++ file that could represent a covariance matrix in a myriad of
> different ways. Then we could parameterize it via the Cholesky factor of
> the precision matrix without the user really knowing.

I agree that having something like a covariance matrix
to cache factorizations makes sense. I'm still not sure
how to do it, though in terms of when things like Cholesky
factors get computed. If we allow general A[m,n] = x;
assignments, then we can't guarantee positive definiteness.
For operations that preserve positive definiteness, we
need to recompute Cholesky factors, etc.

> Also, we can mention that users can rescale
> a correlation matrix. I'll eventually replace
> that transform with the hyperspherical transform
> if it looks like it'll be more efficient.
>
> I worked out that the LKJ and the onion are particular hyperspherical
> transforms.

Interesting. I have to go into Eigen and figure out
what all the pull() and other ops you used were so I
can reverse engineer the transform! Or if it's easy
enough to write, you could add it to the manual (which
I still haven't broken out into chapters yet).

The transform I was considering would take the rows of
the Cholesky factor to have sums of squares equal to 1 (i.e.,
drawn from a unit hypersphere). That'll guarantee you get
a correlation matrix from LL'.

Is that something like what's going on for the
LKJ correlation inverse transform?

> It is just that they are the only hyperspherical transforms
> that imply a known joint density for the correlation matrix. You could
> use another density on the hypersphere(s) and get some correlation
> matrix, but its joint density would just be (proportional to) the
> product of the (independent) hypersphere densities times a Jacobian
> determinant.

> Hard to see what would be the use of implementing that in
> stan until someone does it in such a way that the joint density becomes
> analytically attractive.

I'm just thinking about supporting the corr_matrix data type
most efficiently in terms of inverse transforms and their
Jacobians. The transform only gets applied for initial values.

- Bob

[Bob Carpenter]

I'd done this in one other instance, but can't
recall if I mentioned it on the mailing list.

I validated my own work on the Jacobian determinant of
the Cholesky basis of covar matrices using

stan::agrad::jacobian()

and Eigen's determinant function. This gives a nice
programmatic way of checking one's work on transforms.

The example for covariance matrix constrining (inverse)
transform is in:

src/test/prob/transform_test.cpp

The Jacobian calculator in stan::agrad::jacobian() just
applies the Real^K to Real^K function, then runs K
partial derivative propagation passes, resetting the
adjoints in between each pass.

- Bob

[Ben Goodrich]

On Wednesday, April 25, 2012 11:42:23 AM UTC-4, Bob Carpenter wrote:

By the way, Matt found this book:

http://cs.nenu.edu.cn/Lunwen/2010050513555973499.pdf

on matrix differential calc -- it has all sorts of
great material on Jacobians and derivatives of
matrix operations.

Ah, a classic! Good to see that you can get it online and apparently legitimately.

 > We could have
> sym_pos_def_matrix but that is a bit long-winded. I'd be happy to stick
> with cov_matrix in the .stan file if we had a cov_matrix class in the
> C++ file that could represent a covariance matrix in a myriad of
> different ways. Then we could parameterize it via the Cholesky factor of
> the precision matrix without the user really knowing.

I agree that having something like a covariance matrix
to cache factorizations makes sense.  I'm still not sure
how to do it, though in terms of when things like Cholesky
factors get computed.

I would say the rule should be: If a function needs some representation of the covariance matrix that is not already present, that function should calculate and cache. That way, we just have to check whether something has been initialized before we use it.
 

If we allow general A[m,n] = x;
assignments, then we can't guarantee positive definiteness.

It's very difficult to allow such assignments in a covariance matrix. We can constrain some product-moment correlations to be zero by constraining the corresponding CPC appropriately, but that isn't implemented yet. Nor do we utilize any sparse matrix operations from eigen. What sort of A[m,n] = x assignments are you contemplating besides x being a block of zeros?
 

> I worked out that the LKJ and the onion are particular hyperspherical
> transforms.

Interesting.  I have to go into Eigen and figure out
what all the pull() and other ops you used were so I
can reverse engineer the transform!  

All that pull stuff in transform.hpp basically just reflects that fact that we get "K choose 2" CPCs in an Array and need to stuff them into the lower triangle of a Cholesky factor. As in,

L =
1
a x
b e x
c f  h x
d g i  j x

where x is the square root of the sum-of-squares in its row. So, first we pull a:d from the CPCs, and put them in L. Then remark the spot in CPCs that we pull from, pull e:g, and put them in L, etc. And calculate the diagonal values along the way.

The way CPCs work, the i'th row of L turns out to be a (nonuniform) draw on the surface of a hypersphere with i dimensions. If there were another distribution based on the hyperspheres-to-correlations transformation, it would probably be best to make the columns of L.transpose() be the hyperspheres.

> It is just that they are the only hyperspherical transforms
> that imply a known joint density for the correlation matrix. You could
> use another density on the hypersphere(s) and get some correlation
> matrix, but its joint density would just be (proportional to) the
> product of the (independent) hypersphere densities times a Jacobian
> determinant.

 > Hard to see what would be the use of implementing that in
> stan until someone does it in such a way that the joint density becomes
> analytically attractive.

I'm just thinking about supporting the corr_matrix data type
most efficiently in terms of inverse transforms and their
Jacobians.  The transform only gets applied for initial values.

We already have one function for transforming initial values in transform.hpp called factor_cov_matrix(). It could be based on that.

Ben

[Ben Goodrich]

On Wednesday, April 25, 2012 12:26:26 PM UTC-4, Ben Goodrich wrote:

All that pull stuff in transform.hpp basically just reflects that fact that we get "K choose 2" CPCs in an Array and need to stuff them into the lower triangle of a Cholesky factor. As in,

L =
1
a x
b e x
c f  h x
d g i  j x

where x is the square root of the sum-of-squares in its row. So, first we pull a:d from the CPCs, and put them in L. Then remark the spot in CPCs that we pull from, pull e:g, and put them in L, etc. And calculate the diagonal values along the way.

Err. I meant to say that the lower triangular cells of L are *proportional* to the CPCs. So, (beyond the first column) there is an intermediate step where you update the row-wise products.

Ben

[Bob Carpenter]

On 4/25/12 12:26 PM, Ben Goodrich wrote:
> On Wednesday, April 25, 2012 11:42:23 AM UTC-4, Bob Carpenter wrote:
>
> > We could have
> > sym_pos_def_matrix but that is a bit long-winded. I'd be happy to
> stick
> > with cov_matrix in the .stan file if we had a cov_matrix class in
> the
> > C++ file that could represent a covariance matrix in a myriad of
> > different ways. Then we could parameterize it via the Cholesky
> factor of
> > the precision matrix without the user really knowing.
>
> I agree that having something like a covariance matrix
> to cache factorizations makes sense. I'm still not sure
> how to do it, though in terms of when things like Cholesky
> factors get computed.
>
>
> I would say the rule should be: If a function needs some representation
> of the covariance matrix that is not already present, that function
> should calculate and cache. That way, we just have to check whether
> something has been initialized before we use it.

That makes sense.

> If we allow general A[m,n] = x;
> assignments, then we can't guarantee positive definiteness.
>
>
> It's very difficult to allow such assignments in a covariance matrix. We
> can constrain some product-moment correlations to be zero by
> constraining the corresponding CPC appropriately, but that isn't
> implemented yet. Nor do we utilize any sparse matrix operations from
> eigen. What sort of A[m,n] = x assignments are you contemplating besides
> x being a block of zeros?

Nothing in particular. I'm just imagining someone
doing something like rescaling one entry at a time.

Or just having the covariance matrix be on the left-hand
side of an assignment.

This is all related to whether we only allow operations
that maintain positive definiteness.

> > I worked out that the LKJ and the onion are particular
> hyperspherical
> > transforms.
>
> Interesting. I have to go into Eigen and figure out
> what all the pull() and other ops you used were so I
> can reverse engineer the transform!
>
>
> All that pull stuff in transform.hpp basically just reflects that fact
> that we get "K choose 2" CPCs in an Array and need to stuff them into
> the lower triangle of a Cholesky factor. As in,
>
> L =
> 1
> a x
> b e x
> c f h x
> d g i j x
>
> where x is the square root of the sum-of-squares in its row.

Nice. The hyperspherical transform is really easy if the
x values here (really 4 different ones) are positive.

Am I right in thinking we could also use a basis of the
form:

L =
0
a 0
b e 0
c f h 0
d g i j 0

where each of the row vectors is of length 1 (i.e., on
the hypersphere)? You still get 1s on the diagonal, a positive
definite result, and all values constrained between -1 and 1.
I'm just not sure it's an onto mapping given the way you did it.

> So, first
> we pull a:d from the CPCs, and put them in L. Then remark the spot in
> CPCs that we pull from, pull e:g, and put them in L, etc. And calculate
> the diagonal values along the way.

I figured it was something like that.

A more efficient way to do this for stan::agrad::var
(as things stand now) would be to use our own dot_product
ops to produce the result so the gradients are more efficient.
I'm going to roll that into the covariance matrix case
now that it's unit tested.

> The way CPCs work, the i'th row of L turns out to be a (nonuniform) draw
> on the surface of a hypersphere with i dimensions. If there were another
> distribution based on the hyperspheres-to-correlations transformation,
> it would probably be best to make the columns of L.transpose() be the
> hyperspheres.

Thanks -- I didn't realize the connection to CPCs.

>
> > It is just that they are the only hyperspherical transforms
> > that imply a known joint density for the correlation matrix. You
> could
> > use another density on the hypersphere(s) and get some correlation
> > matrix, but its joint density would just be (proportional to) the
> > product of the (independent) hypersphere densities times a Jacobian
> > determinant.
>
> > Hard to see what would be the use of implementing that in
> > stan until someone does it in such a way that the joint density
> becomes
> > analytically attractive.
>
> I'm just thinking about supporting the corr_matrix data type
> most efficiently in terms of inverse transforms and their
> Jacobians. The transform only gets applied for initial values.
>
>
> We already have one function for transforming initial values in
> transform.hpp called factor_cov_matrix(). It could be based on that.

Isn't that what's already going on now?

- Bob

[Ben Goodrich]

On Wednesday, April 25, 2012 2:19:00 PM UTC-4, Bob Carpenter wrote:

> eigen. What sort of A[m,n] = x assignments are you contemplating besides
> x being a block of zeros?

Nothing in particular.  I'm just imagining someone
doing something like rescaling one entry at a time.

Or just having the covariance matrix be on the left-hand
side of an assignment.

This is all related to whether we only allow operations
that maintain positive definiteness.

Fortunately, there are very few operations that maintain positive definiteness.

Am I right in thinking we could also use a basis of the
form:

L =
0
a 0
b e 0
c f h 0
d g i j 0

where each of the row vectors is of length 1 (i.e., on
the hypersphere)?  You still get 1s on the diagonal, a positive
definite result, and all values constrained between -1 and 1.
I'm just not sure it's an onto mapping given the way you did it.

L * L.transpose() is positive definite iff the diagonal entries of L are positive, so the above would yield a Sigma that is indefinite and singular.

A more efficient way to do this for stan::agrad::var
(as things stand now) would be to use our own dot_product
ops to produce the result so the gradients are more efficient.
I'm going to roll that into the covariance matrix case
now that it's unit tested.

If we are going to multiply L by L.transpose(), then we do want to

-- Use out own dot_product()
-- Not multiply anything by zero
-- Not dot_self() any row because we know it will be 1.0 with partial derivatives of zero
-- Utilize symmetry to only do "K choose 2" dot_products()

But then again, we don't want to multiply L by L.transpose() in practice.

> We already have one function for transforming initial values in
> transform.hpp called factor_cov_matrix(). It could be based on that.

Isn't that what's already going on now?

Possibly.

Ben

[Bob Carpenter]

Interesting. Thanks. I was just thinking that it
had the right number of degrees of freedom.
I'll have to work through why this is the case.

That matrix differential calc book I mentioned
is awesome, by the way. I'm glad to hear it's a classic.
I love the way they're working through calc based on
differentials.

> A more efficient way to do this for stan::agrad::var
> (as things stand now) would be to use our own dot_product
> ops to produce the result so the gradients are more efficient.
> I'm going to roll that into the covariance matrix case
> now that it's unit tested.
>
>
> If we are going to multiply L by L.transpose(), then we do want to
>
> -- Use out own dot_product()

Yes. We will need to add a non-allocating constructor
so we can do the allocations for the component vari
vectors just once.

> -- Not multiply anything by zero

I'll use the length-based pointer-based constructor to
make sure the structural zeros don't come into play.

> -- Not dot_self() any row because we know it will be 1.0 with partial
> derivatives of zero

Definitely -- I hadn't thought of that. I'll
have to think a bit more about how the gradients
need to propagate.

> -- Utilize symmetry to only do "K choose 2" dot_products()

Yes. We'll create a Matrix<var,Dynamic,Dynamic> on
output that has copies of the same underlying vari
on both sides.

This all makes me realize that the matrix multiplication
inside Eigen's solvers is going to be problematic as
we won't be able to get in and make it efficient without
a whole lot of template specializations. I totally think
that's worth doing in the long run, too.

> But then again, we don't want to multiply L by L.transpose() in practice.

I understand we may be able to skirt around this
issue eventually when we introduce a Cholesky-factor
data type. But for right now, the covariance matrix
inverse transform does the multiplication in order to
produce a covariance matrix. And yes, I know that's
bad because it gets factored again, but I really don't
want to start adding data types at this stage.

We don't even have discrete sampling or all the cumulative
densities implemented yet!

You can build a Cholesky factor up from scratch using
the built in positive-constraining inverse transform,
but it's very clunky.

- Bob

[Bob Carpenter]

Here's the top-level priorities as I see them:

1. Get Version 1.0 out the door.

1.A Remaining BUGS models.

1.B Finish the user's guide part of the manual (the basic
reference is done)

1.C Polish up the web site and installation/getting started
instructions so as not to chase away users with low
tolerance for pain.

Then, I think the next big priority is:

2. R interface

2.A data transfer in/out with objects in C++ (Rcpp)

2.B graphics to replace R2WinBUGS (ideally in ggplot2)

3. Convergence monitoring, ESS, etc.
really need this for the R interface

We want these two so that Andrew can start playing around with
it. It really doesn't exist until these are done.

Other big items on the to-do list that seem less critical
to me right away are:

4. Discrete samplers

5. Remaining cumulative densities (required for censoring/truncation)

6. NUTS III
-- need this for lighter and heavier tailed posteriors

7. More multivariate data types & transforms

7.A Tighter integration with Eigen ops, faster matrix ops

7.B New types and densities for Stan's language

8. API doc and user's guide

- Bob

[Ben Goodrich]

On Thursday, April 26, 2012 12:50:56 PM UTC-4, Bob Carpenter wrote:

You can build a Cholesky factor up from scratch using
the built in positive-constraining inverse transform,
but it's very clunky.

Should we have a "hacks" directory somewhere in stan where we collect snippets of .stan syntax to work around limitations of the C++ API? I guess all this is in the manual PDF and .tex file, but it would probably be more useful to have something in a format that users can copy, paste, and modify to their situation.

Ben

[Bob Carpenter]

I'd rather not call it "hacks", but I take your point.

When you say "limitations", you mean things like
not having vectorized forms of the distributions, using
covariance matrices, etc.?

There's also hacking around the limitations of the Stan
programming language itself. I should've added another
two to-do items:

8. Stan Functions

9. Conditionals w. if/then/else, while loops

I think the best thing to do is have working
programs. I'm adding source for everything in
the manual under stan/models.

We could add something to the web page, but
right now, it's very minimal. We took down the
Wiki section because the doc there was already
obsolete. But that'd be a natural place to
put it. Things should be more stable now, too.

- Bob

[Ben Goodrich]

On Thursday, April 26, 2012 1:37:54 PM UTC-4, Bob Carpenter wrote:

I'd rather not call it "hacks", but I take your point.

When you say "limitations", you mean things like
not having vectorized forms of the distributions, using
covariance matrices, etc.?

Yes, things like that. Perhaps also translations of BUGS hacks into hopefully less hackish stan constructs.

But, I don't think we are too far away from being able to vectorize the elliptical distributions. I think all we need at this point is a specialization of dot_self() that accepts a Matrix<T,Dynamic,Dynamic> and returns the column-wise sums of squares.

Ben

[Bob Carpenter]

On 4/26/12 7:01 PM, Ben Goodrich wrote:
> But, I don't think we are too far away from being able to vectorize
> the elliptical distributions. I think all we need at this point is a
> specialization of dot_self() that accepts a Matrix<T,Dynamic,Dynamic>
> and returns the column-wise sums of squares.

To elaborate a bit, that'd be a function with signature:

Matrix<T,Dynamic,1>
columns_dot_self(Matrix<T,Dynamic,Dynamic>);

that fills up the output vector with the dot
products of the columns?

That's certainly easy enough to write.

Go ahead and add it if it'd help. Or I can add
it if you'd rather.

> ... translations of BUGS hacks into

> hopefully less hackish stan constructs.

I'm working on these right now. I'm about
to finish up the cyclic DAG recoding with
bounds.

I need to add sections to the manual on some
other operations.

The thing that still looks like more of a hack
in Stan (and in JAGS) is missing data -- we still
need to break it out. I'm also about to write that
part of the manual.

Do you want
to add it or do you want me to?

- Bob

[Ben Goodrich]

I'll do columns_dot_self and you do missing data. -- Ben

[Ben Goodrich]

I am doing the columns_dot_self(), but to think out loud about missing data for stan 2.0 ...

Currently, we can deal with missing data in the dependent variable to some extent with the if_else() function because whatever arbitrary values are inserted to indicate missingness, they never get accessed by if_else(). But that workaround seems hard to implement for the case where the explanatory variables have missingness.

In general, once Eigen 3.1 is released with expanded sparse matrix support, it would be nice to use sparse matrices where appropriate. So, if we have a covariance matrix with a lot of fixed zeros or whatever, we could use that, save memory due to not allocating a bunch of zero vars and circumvent some of the difficulties of "constant vars" in an auto-diff context.

But could we do something similar to deal with missing data that would be like if_else() on steroids? Say you have a  SparseMatrixXd called X with missing values coded as zero (and thus not allocated) and a SparseMatrix<var,Dynamic,Dynamic> called X_mirror which is like the mirror image of X: observed data in X get coded as zero (and thus not allocated) in X_mirror while missing values in X are vars in X_mirror. Then, stan could sample the vars in X_mirror just like regular parameters, perhaps with an option to not write them to the .csv file.

The model section of the .stan file would have to be a bit more complicated because in addition to having priors over the missing data, you would need to do something like

for(i in 1:N) {
  mu = dot_product(X[i,], beta) + dot_product(X_mirror[i,], beta);
  // likelihood for y[i] | mu
}

And there would be a performance hit from using sparse operations, but I think it would be doable.

Am I missing something (besides data)?

Ben

[Bob Carpenter]

I would not recommend choosing fixed values and
treating them as "unknown". I'm going to modify
the BUGS models that do that to more efficient
forms.

The x + x_mirror thing might work. It still
looks about as clunky as the alternatives.

Instead, I suggested something like this in
the doc section I wrote yesterday on missing data:

data {
int N_obs;
int N_miss;
real y_obs[N_obs];
}
parameters {
real mu;
real(0,) sigma;
real y_miss[N_miss];
}
model {
for (n in 1:N_obs)
y_obs[n] ~ normal(mu,sigma);
for (n in 1:N_miss)
y_miss[n] ~ normal(mu,sigma);
}

Now if both y_obs and y_miss were sparse, they could
both have the same shape and have alternating missingness,
but there's not much point.

The less efficient alternative is to pack both into a
transformed parameter, as in this example, which computes
the same as the above,

data {
int K;
int N_obs;
int N_miss;
real y_obs[N_obs];
}
transformed data {
int N;
N <- N_obs + N_miss;
}
parameters {
real beta[K];
real(0,) sigma;
real y_miss[N_miss];
}
transformed parameters {
real y[N];
for (n in 1:N_obs)
y[n] <- y_obs[n];
for (n in 1:N_miss)
y[n + N_obs] <- y_miss[n];
}
model {
for (n in 1:N)
y[n] ~ normal(mu,sigma);
}

Sometimes you need to do that, for instance if y is
a multivariate parameter to another distribution and
you need it whole.

The second alternative is less efficient because

y[n] <- y_obs[n]

has to promote a double-based y_obs[n] to
a stan::agrad::var for y[n].

- Bob

[Ben Goodrich]

Now that Eigen has essentially full support for sparse matrices, can we revisit this idea for dealing with an arbitrary missingness pattern in the data? By arbitrary missingness pattern, I mean a datamatrix that looks something like this:

         x_1        y_1        y_2         y_3
1  0.1872073  1.3203964 -1.3392774          NA
2  0.5026465 -0.1830042  0.6394578 -0.07292972
3  0.7816282         NA -1.7383923 -0.62743554
4 -1.8949041 -0.1960061  0.2543012  0.18612357
5  0.8035139  1.7151699         NA -0.10530234
6  0.5682540         NA         NA -0.94662401

Stan v1.0.2 can currently proceed iff the data can easily be sorted into blocks of fully observed and fully missing data so that the model{} section can separately loop over each block. The frustrating thing about missingness among several variables is that sometimes no such sorting is possible, and even in special cases where it is possible, it is not particularly easy for a user to do. For example, you could pull rows 2 and 4 into a separate matrix, but what do you do with rows 1, 3, 5, and 6, which have both observed and NA cells?

With this situation in mind, I originally said:

On 5/1/12 11:12 PM, Ben Goodrich wrote: 
> But could we do something similar to deal with missing data that would be like if_else() on steroids? Say you have a
> SparseMatrixXd called X with missing values coded as zero (and thus not allocated) and a
> SparseMatrix<var,Dynamic,Dynamic> called X_mirror which is like the mirror image of X: observed data in X get coded as
> zero (and thus not allocated) in X_mirror while missing values in X are vars in X_mirror. Then, stan could sample the
> vars in X_mirror just like regular parameters, perhaps with an option to not write them to the .csv file.
>
> The model section of the .stan file would have to be a bit more complicated because in addition to having priors over
> the missing data, you would need to do something like
>
> for(i in 1:N) {
>   mu = dot_product(X[i,], beta) + dot_product(X_mirror[i,], beta);
>   // likelihood for y[i] | mu
> }
>
> And there would be a performance hit from using sparse operations, but I think it would be doable.

 

I think I failed to sufficiently emphasize the prevalence / importance / difficulty of the arbitrary missingness pattern, but in response Bob said (with my comments today now mixed in).

 
On Wednesday, May 2, 2012 11:01:51 AM UTC-4, Bob Carpenter wrote:

I would not recommend choosing fixed values and
treating them as "unknown".  I'm going to modify
the BUGS models that do that to more efficient
forms.

This is an inherent downside with my proposal, but I don't see an alternative if Stan cannot read in a NA value and Eigen can't sparse out a NA value. The user would have to shift the partially observed data so that the observed data had no genuine zeros, so that zero would indicate missingness. Although maybe we could have a function in rstan that would facilitate this.
 

The x + x_mirror thing might work.  It still
looks about as clunky as the alternatives.

In my mind, the user would never do x + x_mirror because that would entail both promotion and fill-in. But the user could do

mu <- dot_product(x, beta) + dot_product(x_mirror, beta);

which is just the sum of two vars and the promotion inside dot_product() is no worse than what we have to do anyway.  That is still clunky, but it is the only thing that comes to mind for an arbitrary missingness pattern in the data.
 

Instead, I suggested something like this in
the doc section I wrote yesterday on missing data:

data {
   int N_obs;
   int N_miss;
   real y_obs[N_obs];
}
parameters {
   real mu;
   real(0,) sigma;
   real y_miss[N_miss];
}
model {
   for (n in 1:N_obs)
     y_obs[n] ~ normal(mu,sigma);
   for (n in 1:N_miss)
     y_miss[n] ~ normal(mu,sigma);
}

 
In this case, being able to sort into blocks where y is observed (data) and y is missing (parameters) works fine, but it doesn't generalize to a situation where sorting isn't possible or easy.

Now if both y_obs and y_miss were sparse, they could
both have the same shape and have alternating missingness,
but there's not much point.

In Bob's example above, utilizing the sparse features in Eigen would indeed be overkill.
 

The less efficient alternative is to pack both into a
transformed parameter, as in this example, which computes
the same as the above,

I think this is going to be too inefficient in general (and a pain for the user to code)
 

data {
   int K;
   int N_obs;
   int N_miss;
   real y_obs[N_obs];
}
transformed data {
   int N;
   N <- N_obs + N_miss;
}
parameters {
   real beta[K];
   real(0,) sigma;
   real y_miss[N_miss];
}
transformed parameters {
   real y[N];
   for (n in 1:N_obs)
     y[n] <- y_obs[n];
   for (n in 1:N_miss)
       y[n + N_obs] <- y_miss[n];
}
model {
   for (n in 1:N)
     y[n] ~ normal(mu,sigma);
}

Sometimes you need to do that, for instance if y is
a multivariate parameter to another distribution and
you need it whole.

Yes, but that is a bit of a different situation than partially observed data that do not sort into blocks of fully observed and fully missing data.
 

The second alternative is less efficient because

   y[n] <- y_obs[n]

has to promote a double-based y_obs[n] to
a stan::agrad::var for y[n].

 
Right. The main point of using a sparse vector x (with missing data in x "coded" as zero) and another sparse vector x_mirror of parameters (with zeros in the places that x is observed) is that the zeros aren't allocated / accessed / promoted.

But we would need to introduce some special new syntax to make it as easy as possible on the user. Maybe something like

data {
 int<lower=1> N;
 int<lower=1> K;
 matrix[N,K,has_missingness] X;
 real y[N];
}
parameters {
  matrix[N,K,!X] X_mirror;
  real beta[K];
  real<lower=0> sigma;
  // other parameters
}
model {
  for (i in 1:N) {
    for (j in 1:K) {
      if(X[i,j] == 0) {
        // prior on X_mirror[i,j], which is often some model involving other parameters
      }
    }
  }
  // prior on beta and sigma and other parameters
  y ~ normal(X * beta + X_mirror * beta, sigma);
}

That is hackish and would require a lot of new C++, but it is scalable to arbitrary missingness patterns and consistent with the principles of not allocating more than you need to and not promoting more than you have to.

Ben

[Ben Goodrich]

On Monday, October 8, 2012 1:33:57 AM UTC-4, Ben Goodrich wrote:

Right. The main point of using a sparse vector x (with missing data in x "coded" as zero) and another sparse vector x_mirror of parameters (with zeros in the places that x is observed) is that the zeros aren't allocated / accessed / promoted.

And, perhaps most importantly, the zeros are not auto-diffed.

Ben

[Bob Carpenter]

On 10/8/12 1:33 AM, Ben Goodrich wrote:
> Now that Eigen has essentially full support for sparse matrices,
> can we revisit this idea for dealing with an arbitrary
> missingness pattern in the data?

Sure -- it's just a matter of prioritizing what you're
intersted in doing. All of the things you've put on the table
would be good for Stan:

-- caching covariance data structure
-- LKJ gradients
-- basic Eigen operator traits and expression templates
-- missing data with sparsity hack
-- sparse vectors

-- auto-dif unfolds for linear algebra (that may be me, not you)

Other than the LKJ gradients and auto-dif unfolds for
linear algebra, they're all fairly substantial projects
that involve the guts of Stan's parsing, data/param I/O, and
type inference (assignment, sampling, function calls, operators).

The really big change for the sparsity solution with mirrors would
be handling NA somehow (that's a very fundamental change) and
a whole bunch of I/O and type inference if the mirror solution can't
be reduced to a regular matrix op, which it can't be efficiently.

- Bob

[Ben Goodrich]

On Monday, October 8, 2012 2:17:41 PM UTC-4, Bob Carpenter wrote:

The really big change for the sparsity solution with mirrors would
be handling NA somehow (that's a very fundamental change) and
a whole bunch of I/O and type inference if the mirror solution can't
be reduced to a regular matrix op, which it can't be efficiently.

 
I understand that it would be a big undertaking, but I'm not sure I understand what you are saying here. Under my proposal, the user would have to shift the observed data so that none were zero, recode the NAs as zeros, and somehow indicate in the .stan file that the matrix should have sparse storage. I don't think we would need to directly handle the NAs as special values like R does.

For the X_mirror sparse matrix of missing-values-as-parameters, the basics of what we need are

-- a way to construct / fill it
-- X_mirror * coefficients
-- row(X_mirror, i)
-- col(X_mirror, j)
-- X_mirror(i,j)

and Eigen has all of those now (although X_mirror(i,j) is slow). But users would probably be doing things like

col(X_mirror,j) ~ normal(mu, sigma);

so, we would also have to overload the densities to ignore the zeros in the sparse vector.

There are probably other things that we will need, but there are not that many operations that I would expect to be used on X_mirror if sparse matrices were only being used for missing data. If we also wanted to have sparse covariance matrices and whatnot, then we would have to implement a lot more, but that could wait.

But, since it is a big undertaking no matter what, does anyone have other ideas for a convenient way for users to handle arbitrary missingness patterns?

Ben

[Marcus Brubaker]

On Mon, Oct 8, 2012 at 1:33 AM, Ben Goodrich <goodri...@gmail.com> wrote:

Now that Eigen has essentially full support for sparse matrices, can we revisit this idea for dealing with an arbitrary missingness pattern in the data? By arbitrary missingness pattern, I mean a datamatrix that looks something like this:

         x_1        y_1        y_2         y_3
1  0.1872073  1.3203964 -1.3392774          NA
2  0.5026465 -0.1830042  0.6394578 -0.07292972
3  0.7816282         NA -1.7383923 -0.62743554
4 -1.8949041 -0.1960061  0.2543012  0.18612357
5  0.8035139  1.7151699         NA -0.10530234
6  0.5682540         NA         NA -0.94662401

Stan v1.0.2 can currently proceed iff the data can easily be sorted into blocks of fully observed and fully missing data so that the model{} section can separately loop over each block. The frustrating thing about missingness among several variables is that sometimes no such sorting is possible, and even in special cases where it is possible, it is not particularly easy for a user to do. For example, you could pull rows 2 and 4 into a separate matrix, but what do you do with rows 1, 3, 5, and 6, which have both observed and NA cells?

Maybe I'm missing something, but if you store a set of indices stan can handle arbitrary missing data patterns, though perhaps it's a bit tedious.  E.G., if you had a DxD matrix M with some observed and some missing values:

data {
  int D;
  int N_obs;
  int x_obs[N_obs];
  int y_obs[N_obs];
  real y_obs[N_obs];

  int N_miss;
  int x_miss[N_miss];
  int y_miss[N_miss];
}

parameters {

  real y_miss[N_miss];
}

transformed parameters {

  matrix[D,D] M;
  for (n in 1:N_obs) {
    M[x_obs[n],y_obs[n]] = y_obs[n];
  }
  for (n in 1:N_miss) {
    M[x_miss[n],y_miss[n]] = y_miss[n];
  }
}

model {
  M ~ ...
}

So long as the indexes x_obs/y_obs and x_miss/y_miss are disjoint and complete then everything works out just fine.  It should even be possible to compute N_miss, x_miss and y_miss in the transformed data block with a little trickery which could be made even easier with a few functions.

Cheers,
Marcus