generating a pseduo-random positive definite matrix

Question

I wanted to test a simple Cholesky code I wrote in C++. So I am generating a random lower-triangular L and multiplying by its transpose to generate A.

A = L * Lt;

But my code fails to factor A. So I tried this in Matlab:

N=200; L=tril(rand(N, N)); A=L*L'; [lc,p]=chol(A,'lower'); p

This outputs non-zero p which means Matlab also fails to factor A. I am guessing the randomness generates rank-deficient matrices. Am I right?

Update:

I forgot to mention that the following Matlab code seems to work as pointed out by Malife below:

N=200; L=rand(N, N); A=L*L'; [lc,p]=chol(A,'lower'); p

The difference is L is lower-triangular in the first code and not the second one. Why should that matter?

I also tried the following with scipy after reading A simple algorithm for generating positive-semidefinite matrices:

from scipy import random, linalg
A = random.rand(100, 100)
B = A*A.transpose()
linalg.cholesky(B)

But it errors out with:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.7/dist-packages/scipy/linalg/decomp_cholesky.py", line 66, in cholesky
    c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True)
  File "/usr/lib/python2.7/dist-packages/scipy/linalg/decomp_cholesky.py", line 24, in _cholesky
    raise LinAlgError("%d-th leading minor not positive definite" % info)
numpy.linalg.linalg.LinAlgError: 2-th leading minor not positive definite

I don't understand why that's happening with scipy. Any ideas?

Thanks,
Nilesh.

Answer 1

The problem is not with the cholesky factorization. The problem is with the random matrix L. rand(N,N) is much better conditioned than tril(rand(N,N)). To see this, compare cond(rand(N,N)) to cond(tril(rand(N,N))). I got something like 1e3 for the first and 1e19 for the second, so the conditioning number of the second matrix is much higher and computations will be less stable numerically. This will result in getting some small negative eigenvalues in the ill-conditioned case - to see this look at the eigenvalues using eig(), some small ones will be negative.

So I would suggest to use rand(N,N) to generate a numerically stable random matrix.

BTW if you are interested in the theory of why this happens, you can look at this paper:

http://epubs.siam.org/doi/abs/10.1137/S0895479896312869

Answer 2

As has been said before, eigen values of a triangular matrix lie on the diagonal. Hence, by doing

L=tril(rand(n))

you made sure that eig(L) only yield positive values. You can improve the condition number of L*L' by adding a large enough positive number to the diagonal, e.g.

L=L+n*eye(n)

and L*L' is positive definite and well conditioned:

> cond(L*L')

ans =

1.8400