Facing a unbalanced matrix, such as 3*5000, why someone reduce the computational cost by svd(A*A') .

Computing singular value decomposition is the main computational cost in many algorithms . For a matrix A(m*n) ,if m is much larger than n , one can compute the SVD of A*A',and then get an approximate SVD of by simple operations to reduce the computational cost. How does it reduce the computational cost?

The code is as follows.

function [S, V, D, Sigma2] = MySVD(A)
[m, n] = size(A);
if 2*m < n
    AAT = A*A';
    [S, Sigma2, D] = svd(AAT);
    Sigma2 = diag(Sigma2);
    V = sqrt(Sigma2);
    tol = max(size(A)) * eps(max(V));
    R = sum(V > tol);
    V = V(1:R);
    S = S(:,1:R);
    D = A'*S*diag(1./V);
    V = diag(V);
    return;
end
if m > 2*n
    [S, V, D, Sigma2] = MySVD(A');
    mid = D;
    D = S;
    S = mid;
    return;
end
[S,V,D] = svd(A);
Sigma2 = diag(V).^2;

Answer 1

[Christine Tobler]

The computational complexity of svd (when using the 'econ' option, which is very necessary for matrices that are far from square) is O(max(m,n)*min(m,n)^2). This can be seen as applying QR to A in a first step, and then computing the SVD of the R matrix coming out of this QR operation (which is relatively cheap in comparison). No need for you to call QR in MATLAB code, this is done internally in the SVD command already.

Instead computing A*A' / A'*A (whichever is smaller based on A's size) has the same computational complexity as above. However, matrix multiply is a cheaper operation than QR, so in practice it's going to be less expensive to apply this.

But this computation is also less accurate: Computing A*A' squares the condition number, so the SVD is applied to a matrix that is worse conditioned and will be less accurate than A applied to the original matrix.

[Bjorn Gustavsson]

But if "an arbitrary matrix" is that are that "far from square" would surely span the small-dimensioned space rather well? Wouldn't that hint that the condition-number of the right A*A'/A'*A often is reasonably OK?

[Christine Tobler]

Well, that's going to depend on the application. The code above is specifically checking for the rank R based on a tolerance, which I take to mean they sometimes have low-rank cases happen.

I'd think the tol there should be chosen as sqrt(eps) instead of just eps, since that's about the loss in accuracy to be expected when using A'*A instead of A. If that's a problem will be dependent on the application, it's irrelevant if a low-rank approximation up to tolerance sqrt(eps) is sufficient.

Here's an example:

A = randn(1e3, 20) * diag(logspace(0, -16, 20)) * randn(20);
semilogy([svd(A), sqrt(svd(A'*A))])
title("Singular values of A")
yline(sqrt(eps)*norm(A))
legend("Computed directly", "Computed based on A'*A", "sqrt(eps)*norm(A)")

[鹏程 张]

Thank you for your answer! But I have one more question .

Assumming that we get U, LandV by the command `[U L V] = svd(A*A')`. How can we use a formula for U, LandV to represent A ?

[Christine Tobler]

It's easier if you use EIG instead of SVD:

A = randn(3, 10);
[U, D] = eig(A*A');
norm(U*D*U' - A*A')

ans = 4.2690e-15

% The SVD of A is U*S*V', so inserting this gives us
% A*A' = U*S*V'*V*S*U' = U*S*(V'*V)*S*U' = U*S*S*U'
% So the U of SVD(A) and the U of EIG(A*A') can match, and D = S*S.
% Here we assume all diagonal elements of D are non-negative, 
% which is only true numerically if A is well-conditioned.
S = sqrt(D);

% Now let's go back to A = U*S*V'. Given we know U and S now, and assuming
% S has a positive diagonal (that is, again, that A is well-conditioned),
% we can reformulate this to
% U'*A = S*V' <=> inv(S)*U'*A = V' <=> A'*U*inv(S) = V
V = (A'*U)/S;
norm(A - U*S*V')

ans = 7.6329e-16

[鹏程 张]

Thank you for your answer again, which initially solve my confusion!