Matrix Low Rank Approximation using Matlab

Question

Consider a 256 x 256 matrix A. I'm familiar with how to calculate low rank approximations of A using the SVD.

Typically after using [U S V] = svd(A), I would use Ak = U(:,1:k)*S(1:k,1:k)*V(:,1:k)'; to get the rank k approximation of A.

My question is how do I create a vector E such that, E(k) = norm(A-Ak) for k=1,2,3.....,256. That is E is a column vector of 256 elements each of which is norm(A-Ak)

Answer 1

The answer is simply

diag(S)

Why?

There's a theorem¹ that says that the error between a matrix A and its rank-k approximation Ak has a (spectral) norm² given by the k+1-th singular value of A. That is, the error is given by the first not used singular value. Isn't it a wonderful³ result?

Example:

>> A = randn(8,8);
>> [U S V] = svd(A);
>> k = 5;
>> Ak = U(:,1:k)*S(1:k,1:k)*V(:,1:k)'; %'// rank-5 approximation 
>> norm(A-Ak) %// its associated error norm
ans =
    1.0590

>> k = 6;
>> Ak = U(:,1:k)*S(1:k,1:k)*V(:,1:k)'; %'// rank-6 approximation
>> norm(A-Ak) %// its associated error norm
ans =
    0.3924

>> diag(S).' %'// all error norms
ans =
    4.5528    3.2398    2.5863    2.2031    1.4252    1.0590    0.3924    0.1021

---

_{¹ Actually I didn't know about such theorem until a few minutes ago. I just computed norm(A-Ak) and noticed the resulting value was in S. Then I thought there must be a theorem that established this.}

_{² Thanks to @AlgebraicPavel for the correction.}

_{³ "Algebra is generous; she often gives more than is asked of her." — Jean le Rond d'Alembert}