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.
- 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) )
- 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!