help with non-square Cholesky factors

[Bob Carpenter]

I remember you wanted this, but I still don't even
understand what is a non-square Cholesky factor is.

QUESTION FOR EVERYBODY:

1. Do we need non-square Cholesky factors? (I'm assuming
the answer is yes to this, though I don't know why.)

2. Which syntax:

a. two data types

cholesky_factor[K]
rectangular_cholesky_factor[M,N]

b. one data type

cholesky_factor[M,N]

c. one data type overloaded for dimensions

cholesky_factor[K]
cholesky_factor[M,N]

b. is the easiest to doc assuming the answer to 1. is yes.

I couldn't find anything in the Wikipedia under Cholesky
factorization or triangular matrices.

What are the constraints? Specifically, if the Cholesky factor
is M x N, what are:

* constraints on M and N

* elements constrained to be zero

* elements constrained to be positive

For square ones, my understanding is that M == N, L[m,n] == 0 if n > m,
and L[m,m] > 0.

From the example, I gather M >= N, L[m,n] == 0 if n > m, and L[m,m] > 0.
Is that right?

- Bob

On 7/28/13 9:49 PM, bgoodri wrote:
> Remember that we should allow Cholesky factors to be non-square. E.g. in R
>
> |> L <- matrix(0, 5, 2)
>> L[lower.tri(L)] <- rnorm(7)
>> diag(L) <- rexp(2)
>> L
> [,1] [,2]
> [1,] 0.2732593 0.0000000
> [2,] -0.1338806 1.3074583
> [3,] 0.4872495 -0.1464854
> [4,] 1.0683007 0.4246451
> [5,] -1.3734227 1.6973773
> |
>
> Maybe we want to call a non-square Cholesky factor something other than a Cholesky factor, but I think it would be fine
> to just call them all Cholesky factors regardless of whether they are square.
>
> —
> Reply to this email directly or view it on GitHub
> <https://github.com/stan-dev/stan/commit/79e98c990b6e47bc82dd8cdc20d03647f3e6766e#commitcomment-3734554>.
>

[Ben Goodrich]

On Monday, July 29, 2013 1:34:42 PM UTC-4, Bob Carpenter wrote:

I remember you wanted this, but I still don't even
understand what is a non-square Cholesky factor is.

If L is MxN with M > N, then

Sigma = L * L'

is a MxM covariance matrix that happens to be singular (actually that is true regardless of the constraints on L). Ostensibly, you could take the LDL' factorization of that Sigma (Eigen actually does pivoting for numerical stability so it would be PLDL'P' where P is a permutation matrix). If you construct the matrix

L * sqrt(D)

then you get a MxM lower-triangular matrix whose last M - N columns are all zero. You could then drop the zero columns to get a MxN Cholesky factor.

QUESTION FOR EVERYBODY:

1.  Do we need non-square Cholesky factors?  (I'm assuming
the answer is yes to this, though I don't know why.)

Yes, they are used in factor models (and some other economics stuff) where you have the functional form

Sigma = Lambda * Lambda' + Theta

where Lambda is lower-triangular with positive diagonal elements but not necessarily square and Theta is diagonal with positive diagonal elements. Then, you can write the likelihood of a sample covariance matrix given Sigma.
 

2.  Which syntax:

     a.  two data types

         cholesky_factor[K]
         rectangular_cholesky_factor[M,N]

     b.  one data type

         cholesky_factor[M,N]

     c.  one data type overloaded for dimensions

         cholesky_factor[K]
         cholesky_factor[M,N]

b. is the easiest to doc assuming the answer to 1. is yes.

I think c is best but b is okay, particularly if the parser could say something like "you put cholesky_factor[K], maybe you meant cholesky_factor[K,K]?"
 

I couldn't find anything in the Wikipedia under Cholesky
factorization or triangular matrices.

What are the constraints?  Specifically, if the Cholesky factor
is M x N, what are:

   * constraints on M and N

I think M >= N is okay. I have never seen the M < N case actually used, and it would be weird because there would be more free elements in L than in Sigma.
 

   * elements constrained to be zero

All cells where the column index is greater than the row index.

   * elements constrained to be positive

All cells where the column index is equal to the row index.

For square ones, my understanding is that M == N, L[m,n] == 0 if n > m,
and L[m,m] > 0. 

 From the example, I gather M >= N, L[m,n] == 0 if n > m, and L[m,m] > 0.
Is that right?

Yes.

Ben

[Bob Carpenter]

Thanks.

I'm going to try to get the square version working then
upgrade.

A lot of code gets touched by the declared data types and I
haven't added one in ages, so needed to refresh myself on
the layout anyway.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Ben Goodrich]

On Monday, July 29, 2013 2:41:16 PM UTC-4, Bob Carpenter wrote:

A lot of code gets touched by the declared data types and I
haven't added one in ages, so needed to refresh myself on
the layout anyway.

In which case, we probably need a wiki page on adding data types.

Another thing that we need is the Cholesky factor such that

L * L'

is a correlation matrix. Attached is a R function that does so in the square case so that it doesn't get lost.

Ben

[Bob Carpenter]

On 7/29/13 2:47 PM, Ben Goodrich wrote:
> On Monday, July 29, 2013 2:41:16 PM UTC-4, Bob Carpenter wrote:
>
> A lot of code gets touched by the declared data types and I
> haven't added one in ages, so needed to refresh myself on
> the layout anyway.
>
>
> In which case, we probably need a wiki page on adding data types.

:-) It'll have to get in line behind the other doc projects
in the works (like the Stan paper, RStan paper and auto-diff doc,
all of which I'd rather do first). None of the relevant code is
documented conceptually, so this would require
documenting the readers, writers and transforms as a basis.

It's probably a higher priority to document how to write a distribution,
too. That's a lot more complicated with the vectorization.

You can find the places where code needs to go by grepping for
cov_matrix.

First I have to finish testing the square case. Then I'll try
to get the non-square case working.

> Another thing that we need is the Cholesky factor such that
>
> |
> L * L'
> |
>
> is a correlation matrix. Attached is a R function that does so in the square case so that it doesn't get lost.

We'll need to work out the Jacobian.

And the non-square case if you want that.

We'll also need a name for the correlation matrix case when I get
around to coding it. Should we anticipate this type and adopt
a suitable naming convention to begin with for the cov case?

cholesky_factor_rectangular_cov[M,N]
cholesky_factor_cov[M]

cholesky_factor_rectangular_corr[M,N]
cholesky_factor_corr[M]

I'm still not sure about overloading the syntax for square and
rectangular versions --- I'm worried it'll be confusing.

This might all have to wait until after 2.0 --- it's turning
into a huger project than the couple of days I anticipated.

I was hoping to also get to the lp__ replacement, but maybe that
and correlation matrices Cholesky factors can wait until after 2.0.

I was also hoping to get forward-mode auto-diff done for 2.0,
but that's probably also going to have to wait.

- Bob

[Ben Goodrich]

On Monday, July 29, 2013 3:51:09 PM UTC-4, Bob Carpenter wrote:

> Another thing that we need is the Cholesky factor such that
>
> |
> L * L'
> |
>
> is a correlation matrix. Attached is a R function that does so in the square case so that it doesn't get lost.

We'll need to work out the Jacobian.

It is pretty easy, or at least easier than the Jacobians in transform.hpp.
 

And the non-square case if you want that.

Yeah, not hard.
 

We'll also need a name for the correlation matrix case when I get
around to coding it.  Should we anticipate this type and adopt
a suitable naming convention to begin with for the cov case?

   cholesky_factor_rectangular_cov[M,N]
   cholesky_factor_cov[M]

   cholesky_factor_rectangular_corr[M,N]
   cholesky_factor_corr[M]

I would just go with

cholesky_factor_cov[M,N]
cholesky_factor_corr[M,N]

if the parser can help the user that only specifies [M].

Ben

[Bob Carpenter]

On 7/29/13 4:03 PM, Ben Goodrich wrote:
> On Monday, July 29, 2013 3:51:09 PM UTC-4, Bob Carpenter wrote:
>
> > Another thing that we need is the Cholesky factor such that
> >
> > |
> > L * L'
> > |
> >
> > is a correlation matrix. Attached is a R function that does so in the square case so that it doesn't get lost.
>
> We'll need to work out the Jacobian.
>
>
> It is pretty easy, or at least easier than the Jacobians in transform.hpp.

Can we re-use part of the correlation matrix transform that's
already in there? If so, it should be much easier. Basically,
the Cholesky factor was just like stopping short on the covariance
matrix transforms.

> And the non-square case if you want that.
>
>
> Yeah, not hard.

Great.

> We'll also need a name for the correlation matrix case when I get
> around to coding it. Should we anticipate this type and adopt
> a suitable naming convention to begin with for the cov case?
>
> cholesky_factor_rectangular_cov[M,N]
> cholesky_factor_cov[M]
>
> cholesky_factor_rectangular_corr[M,N]
> cholesky_factor_corr[M]
>
>
> I would just go with
>
> cholesky_factor_cov[M,N]
> cholesky_factor_corr[M,N]
>
> if the parser can help the user that only specifies [M].

That seems cleanest.

It's no problem for the parser to have a separate case for a single
one that repeats the variable. The question is whether the users
will find it confusing. Probably no more confusing than having two
different names, so let's go with this approach. I can rename
the one I have when I generalize to rectangular.

The square covariance matrix branch is feature complete, unit
tested, and pushed --- I haven't tried any models, but it gives
a good example of what you need to do to add a new data type.

I'm loving this Git branching stuff. Took a long time to grok,
but now I don't know how I ever worked without it.

- Bob

[Ben Goodrich]

On Monday, July 29, 2013 5:04:14 PM UTC-4, Bob Carpenter wrote:

On 7/29/13 4:03 PM, Ben Goodrich wrote:
> It is pretty easy, or at least easier than the Jacobians in transform.hpp.

Can we re-use part of the correlation matrix transform that's
already in there?  If so, it should be much easier.  Basically,
the Cholesky factor was just like stopping short on the covariance
matrix transforms.

Not that much, but we could change the transformations (mostly based on the idea of CPCs) that we have now to the simpler transformations (mostly based on spherical distributions) that I used in onion.R since I have fully figured out how both relate to each other and to the LKJ paper. The onion way is better.

Ben

[Bob Carpenter]

On 7/29/13 2:47 PM, Ben Goodrich wrote:

...

> Another thing that we need is the Cholesky factor such that
>
> |
> L * L'
> |
>
> is a correlation matrix.

> Attached is a R function that does so in the square case so that it doesn't get lost.

Thanks.

Could you write the correlation matrix version of the three
cholesky_factor_cov functions in src/stan/prob/transform.hpp?
These are:

1. constrain (vector + dims to matrix)

2. constraint + Jacobian adjust (vector + dims to matrix, increment log prob)

3. unconstrain (matrix to vector)

There are examples of the tests I used including error conditions
in the corresponding test file.

Everything else is just a direct copy of the cholesky_factor_cov case
that I can knock off in no time.

- Bob

[Ben Goodrich]

On Wednesday, July 31, 2013 12:01:09 PM UTC-4, Bob Carpenter wrote:

Could you write the correlation matrix version of the three
cholesky_factor_cov functions in src/stan/prob/transform.hpp?
These are:

   1.  constrain  (vector + dims to matrix)

   2.  constraint + Jacobian adjust (vector + dims to matrix, increment log prob)

   3.  unconstrain (matrix to vector)

I think the best thing for Stan 2.0 is to use what we have now in transform.hpp for Cholesky factors of correlation matrices, and maybe we'll figure out something with the onion method later. We already have

1. read_corr_L(CPCs, K)

2. read_corr_L(CPCs, K, log_prob)

so I think I just need to add 3. This just entails pulling out the guts of factor_cov_matrix().

 

There are examples of the tests I used including error conditions
in the corresponding test file.

Everything else is just a direct copy of the cholesky_factor_cov case
that I can knock off in no time.

OK. Do you know the SHAs of the commits that did the cholesky_factor_cov case? Things are harder to find after a feature branch gets deleted.

Another thing that would be good to have for 2.0 is versions of wishart_log() and inv_wishart_log() that take the L such that L * L' is distrubuted accordingly, which will entail another Jacobian adjustment but it is one that we already know. Similarly, a version of lkj_corr() that takes the L such that L * L' is distributed LKJ.

Ben