eig() gives a negative eigenvalue for a positive semi-definite matrix

Based on the following posted questions on StackExchange:

https://math.stackexchange.com/questions/1463140/proof-for-why-a-matrix-multiplied-by-its-transpose-is-positive-semidefinite https://math.stackexchange.com/questions/473432/positive-definite-and-transpose

I think it is safe to conclude that a rectangular matrix A times its transpose results in a square matrix that is positive semi-definite.

Here is my problem:

A = [-0.0243, 0.0053; 0.0103, 0.0033; 0.0171, 0.0011];
S = A * A';
eig(S);

One of the eigenvalues is negative, though it is very close to zero.

Where am I wrong? Thank you!

Answer 1

[Christine Tobler]

While A*A' is mathematically symmetric positive semi-definite, its numeric representation will not be - it's more like A*A' + (round-off error). This error makes the smallest eigenvalue close to but not exactly zero, which is as likely to be positive as negative.

An alternative would be to compute the singular value decomposition, A = U*S*V', which will give you the eigenvalue decomposition:

 A*A' = U*S*V'*V*S*U' = U*S^2*U'.

Here, S is guaranteed to be non-negative, and it's also more robust because it doesn't have the numeric error from computing A*A'.

[Zachary L]

Thank you!

I kind of understand your point. But my main concern is that eig(S) will yield negative values, and this prevents me to do chol(S). I think you are right that singular decomposition is more robust, but it still can't get rid of getting negative eigenvalues, for example:

A = [-0.0000243, 0.0000053; 0.0000103, 0.0000033; 0.0000171, 0.0000011];
[U,S] = svd(A);
M = U * S * S' * U';
eig(M);

[Christine Tobler]

Ah, so you're not looking to compute the eigenvalues, but the Cholesky decomposition? The problem here is that Cholesky doesn't work for semi-definite - it actually requires the matrix to be positive definite.

You can still compute a decomposition of A*A' into a product of two triangular matrices:

 >> A = [-0.0243, 0.0053; 0.0103, 0.0033; 0.0171, 0.0011];
 >> [Q, R] = qr(A')
 Q =
   -0.9770    0.2131
    0.2131    0.9770
 R =
    0.0249   -0.0094   -0.0165
         0    0.0054    0.0047

 >> norm(A*A' - R'*R)
 ans =
   2.7205e-19

The matrix Q can be ignored. While R is not square, you can always just add a nonzero row at the bottom if this is a problem.

One advantage here is that you get exact zeros on the diagonal of R for the exactly zero eigenvalues of A*A'.

Answer 2

[Walter Roberson]

A*A' is rank 2, not rank 3, which is numerically unstable to process.

You can

eig(sym(A) * sym(A)')

Keep in mind though that there are no exact binary representations for any of the decimal values you list for A, and the round-off error can have significant effect in such cases.

[Zachary L]

Thank you!

Now I understand the main source of error is rounding

Now suppose I have a positive semi-definite square matrix S with small entries and I want to do

chol(S);

I tried:

chol(sym(S));

But it is very slow when S is at large scale, for example:

B = rand(90,100);
C = B*B';
eig(sym(C));

[Walter Roberson]

You should probably be doing

Bs = sym(B);
C = Bs * Bs';
eig(C)

to reduce rounding errors.

You might want to reduce digits() to speed up the computation.