Properties of SVD of a hermitian matrix not holding at single precision

Greetings,

I came accross something I did not expect and I was hoping for some insight.

I am working on a project that includes dozens of very large and dense hermitian matrices. To try and reduce the amount of memory used by my code, I am storing these matrices at single precision and I am also attempting to use a low-rank approximation of these matrices during execution.

Since these matrices are hermitian, I am anticipating the following to be true:

1. I should be able to express these matrices as, where S and V come from the singular value decomposition of A.( [~, S, V] = svd(A) )
2. For a rank k approximation of A,

Ak = V(:,1:k) * S(1:k,1:k) * V(:,1:k)';

That the frobenius norm of A-Ak should be equal to the k+1 singular value (that is, norm(A-Ak) == S(k+1) )

What I am finding is that 2) is not holding true if I run the svd() function on matrix A, when it is in single precision. However, if I convert the arrays to double before running the SVD, 2) holds true. It also holds true if I convert S and V back to single precision after running the SVD function.

Is this really a precission issue? The smallest singular value of my matrice is 1e-6 so I wasn't expecting a single precission array to be an issue. Any insight would be greatly appreciated.

Thank you in advance!

[David Goodmanson]

HI Michael,

What about step 1?

m = [2 3 4; 3 5 1; 4 1 8]
[~, s, v] = svd(m)
m1 = v*s*v.'

m =
     2     3     4
     3     5     1
     4     1     8
s =
   10.9161         0         0
         0    5.0000         0
         0         0    0.9161
v =
   -0.4818    0.1690   -0.8598
   -0.3780    0.8452    0.3780
   -0.7906   -0.5071    0.3433

   m1 =                          % definitely not equal to m
    3.3544    2.4046    3.4591
    2.4046    5.2617    1.2378
    3.4591    1.2378    8.2160

[Dana]

David brings up a good point: in the SVD [U,S,V]=svd(A), I believe the property U=V holds only if A is positive-definite (which is not the case for m in David's example). Michael, is it the case in your application that, in addition to being Hermitian, all of your matrices are also positive-definite? If not, you should be using the approximation Ak = U(:,1:k) * S(1:k,1:k) * V(:,1:k)' (i.e., U in place of the first V).

That still doesn't explain why you got the expected answer with doubles but not with singles, though. Can you provide a minimum working example showing the problem?

[Michael]

Hi Dana and David,

Thank you so much for your time!

David, my sincere apologies, the matrices I am working with are also positive definite which I didn't think of mentioning in my post. To my knowledge, 1) only holds for positive definite hermitian matrices.

>> m = [7 3 4; 3 5 1; 4 1 8]

m =

     7     3     4
     3     5     1
     4     1     8

>> chol(m)

ans =

    2.6458    1.1339    1.5119
         0    1.9272   -0.3706
         0         0    2.3616

>> ishermitian(m)

ans =

  logical

   1

>> [~, S, V] = svd(m);
>> m1 = V*S*V'

m1 =

    7.0000    3.0000    4.0000
    3.0000    5.0000    1.0000
    4.0000    1.0000    8.0000

Dana, I will look into recreating it with a minimum working example.

Thank you both again!

[Dana]

Michael, you don't need to upload your specific matrix at this point. If this is actually a problem, it should be reproducible with an arbitrary (randomly generated) matrix. Start with posting the code that compares the results between single and double and finds them to be different.

[David Goodmanson]

Hi Michael,

you have

norm(A-Ak) == S(k+1) 

However,

n = 5;
A = diag(n:-1:1)     % an hermitian positive definite matrix
[U S V] = svd(A);    % naturally both U and V are the identity,
Sdiag = diag(S)      % and S is n:-1:1

k = 3
Ak = V(:,1:k) * S(1:k,1:k) * V(:,1:k)'
B = A-Ak;
fn = norm(B(:))      % frobenius norm
Sdiag(k+1)           % not equal to fn unless k = n-1

[Michael]

Hi David,

Thanks for bringing this to my attention. I had not read my original research as carefully as I should have. The spectral norm of Ak- A will be equal to the k+1 singular value, not the frobenius norm. This is indeed the case if you make:

fn = norm(B)

in your example. This gives me food for thought... thank you again!

M

[David Goodmanson]

Hi Michael,

now that the problem is well defined, from your results it appears that you have a prima facie case that single precision is not enough precision.

I would not put that much import into the 1e-6 singular value. The question is, compared to what? If the matrix has a large condition number then it could easily happen that single precision does not get the job done.

[Bruno Luong]

Stupid question but do you take care of M numerically hermitian by running

A=0.5*(A+A')

before anything?

As David rightly pointout, spectral norm equality holds, not Frobenius norm (the later has little meaning interpretation).

Also don't expect SVD to return U == V, even for Hermitian matrix A, there is always an arbitrary phase-shift (sign) in SVD so that

U*diag(epx(1i*theta)) ~= V

with theta arbitrary (real) phase.

To get S and V such that V*S*V' = A

run

[V,S]=eig(A)

Answer 1

[Christine Tobler]

I'd suggest just calling EIG

[V, S] = eig(A)

instead of calling the SVD. If A is (exactly!) symmetric on input, this will return S and V such that V*S*V' == A, and you can check if A is numerically S.P.D. by seeing if S contains any negative numbers (if A is symmetric positive semi-definite, there are likely to be some negative up to round-off values). With the SVD, such a small number would be returned as a positive, and the associated vectors in U and V would have different signs.

If there are values on S's diagonal that are negative at round-off level, you can decide if you have to stop the algorithm, or if it's all right to just set these diagonal values as zero.