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.

In the numpy (which is used by scipy), multiplication of arrays with `*` is element-wise, not algebraic. For the algebraic product, use the `dot` function. So you should write `B = dot(A, A.transpose())`, or `B = dot(A, A.T)`. — Warren Weckesser, Sep 08 '12 at 00:17
I'd generate a standard Gaussian A and do A * A^t. This will yield a Wishart distributed matrix, much better conditioned for this kind of experiments (eigenvalues have xhi^2 distribution IIRC). — Alexandre C., Sep 08 '12 at 10:27

score 6 · Accepted Answer · edited May 15 '14 at 20:58

6

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

edited May 15 '14 at 20:58

user79950

259
3
11

answered Sep 07 '12 at 19:09

Bitwise

7,577
6
33
50

Great! I did not think of that, so although the algorithm itself is stable, the ill-conditioned data causes the problem. Pretty interesting phenomenon described in the paper, although I'll need more background to understand it entirely. Thanks! Another question is, does the second method(rand(N, N)) provably generate SPD matrices? I am thinking there is a (extremely?) small probability that it will generate a rank-deficient matrix. Any thoughts? I don't know what's going on with SciPy though, the condition numbers seem OK for the generated random matrix (1e3). – Nilesh Sep 07 '12 at 21:42
It is actually pretty simple to see that with the second method you can still get an ill-conditioned matrix and get negative eigenvalues. Since it generates a random matrix, you can basically get any matrix (within 0,1) including ill-conditioned ones, and specifically you can get ones which are very similar to the ones which tril(rand(N,N)) gives you (with a small probability). For example: the matrix tril(rand(200,200))+0.0001 can probably be generated by rand(200,200) but will be ill-conditioned. – Bitwise Sep 15 '12 at 18:16

angainor · Answer 2 · 2012-09-13T18:26:48.417

1

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

edited Sep 13 '12 at 18:26

answered Sep 08 '12 at 10:09

angainor

11,760
2
36
56

score 0 · Answer 3 · answered Sep 07 '12 at 17:49

0

To generate a random positive definite matrix in MATLAB your code should read:

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

And you should indeed have the eigenvalues be greater than zero and p be zero.

answered Sep 07 '12 at 17:49

Malife

303
2
6

yes, I forgot to mention that in the original question, now updated. Thanks! But I don't understand why A cannot be factored if L is lower-triangular as in the first Matlab code sample? – Nilesh Sep 07 '12 at 18:34

score 0 · Answer 4 · answered Sep 07 '12 at 23:17

You ask about the lower triangular case. Lets see what happens, and why there are problems. This is often a good thing to do, to look at a test case.

For a simple 5x5 matrix,

L = tril(rand(5))
L =
      0.72194            0            0            0            0
     0.027804      0.78422            0            0            0
      0.26607     0.097189      0.77554            0            0
      0.96157      0.71437      0.98738      0.66828            0
     0.024571     0.046486      0.94515      0.38009     0.087634

eig(L)
ans =
     0.087634
      0.66828
      0.77554
      0.78422
      0.72194

Of course, the eigenvalues of a triangular matrix are just the diagonal elements. Since the elements generated by rand are always between 0 and 1, on average they will be roughly 1/2. Perhaps looking at the distribution of the determinant of L will help. Better is to consider the distribution of log(det(L)). Since the determinant will be simply the product of the diagonal elements, the log is the sum of the logs of the diagonal elements. (Yes, I know the determinant is a poor measure of singularity, but the distribution of log(det(L)) is easily computed and I'm feeling too lazy to think about the distribution of the condition number.)

Ah, but the negative log of a uniform random variable is an exponential variate, in this case an exponential with lambda = 1. The sum of the logs of a set of n uniform random numbers from the interval (0,1) will by the central limit theorem be Gaussian. The mean of that sum will be -n. Therefore the determinant of a lower triangular nxn matrix generated by such a scheme will be exp(-n). When n is 200, MATLAB tells me that

exp(-200)
ans =
   1.3839e-87

Thus for a matrix of any appreciable size, we can see that it will be poorly conditioned. Worse, when you form the product L*L', it will generally be numerically singular. The same arguments apply to the condition number. Thus, for even a 20x20 matrix, see that the condition number of such a lower triangular matrix is fairly large. Then when we form the matrix L*L', the condition will be squared as expected.

L = tril(rand(20));

cond(L)
ans =
   1.9066e+07

cond(L*L')
ans =
   3.6325e+14

See how much better things are for a full matrix.

A = rand(20);

cond(A)
ans =
       253.74

cond(A*A')
ans =
        64384

generating a pseduo-random positive definite matrix

4 Answers4