Instead of [U, S, V] = svd(M)
, which tries to build a matrix U
that is 49152 by 49152 (= 18 GB !), do svd(M, 'econ')
. That returns the “economy-class” SVD, where U
will be 52 by 52, S
is 52 by 52, and V
is also 52 by 52.
cov(M)
will remove each dimension’s mean and evaluate the inner product, giving you a 52 by 52 covariance matrix. You can implement your own version of cov
, called mycov
, as
function [C] = mycov(M)
M = bsxfun(@minus, M, mean(M, 1)); % subtract each dimension’s mean over all observations
C = M' * M / size(M, 1);
(You can verify this works by looking at mycov(randn(49152, 52))
, which should be close to eye(52)
, since each element of that array is IID-Gaussian.)
There’s a lot of magical linear algebraic properties and relationships between the SVD and EVD (i.e., singular value vs eigenvalue decompositions): because the covariance matrix cov(M)
is a Hermitian matrix, it’s left- and right-singular vectors are the same, and in fact also cov(M)
’s eigenvectors. Furthermore, cov(M)
’s singular values are also its eigenvalues: so svd(cov(M))
is just an expensive way to get eig(cov(M))
, up to ±1 and reordering.
As @rayryeng explains at length, usually people look at svd(M, 'econ')
because they want eig(cov(M))
without needing to evaluate cov(M)
, because you never want to compute cov(M)
: it’s numerically unstable. I recently wrote an answer that showed, in Python, how to compute eig(cov(M))
using svd(M2, 'econ')
, where M2
is the 0-mean version of M
, used in the practical application of color-to-grayscale mapping, which might help you get more context.