How does Matlabs svd() function calculate the value of V?

I have been attempting to recreate the singular value decomposition in Matlab without the use of the svd() function. I have been able to get the values of S and U, but I have had issues recreating the results of V. I have been using the definition of V= orthogonal basis of X' * X, where X has dimensions of M * N, represented in code as

V=orth(transpose(X) * X);

However, the output of this is a N * M matrix, while the expected dimensions (and the dimensions of V in the svd function) is M * M. I seem to be missing the last colum of the matrix consistantly. Am I misunderstanding the theory behind svd, or is there an issue with the orth() command?

[David Goodmanson]

Hi Caleb,

In general if n>m and A is an nxm matrix of orthonormal vectors (as might be produced by orth, for example), then

[A null(A')]

is an nxn unitary matrix.

[Steven Lord]

I have been attempting to recreate the singular value decomposition in Matlab without the use of the svd() function.

Why? Pedagogical reasons? To try to implement it faster than the built-in svd function? To try to compute the SVD on a system where you don't have MATLAB available (for example on a piece of hardware)?

In that last case, note that from the svd documentation page that the function supports the C/C++ Code Generation extended capability, though as the usage note states "Code generation uses a different SVD implementation than MATLAB uses. Because the singular value decomposition is not unique, left and right singular vectors might differ from those computed by MATLAB."

...

Am I misunderstanding the theory behind svd, or is there an issue with the orth() command?

As an aside, using orth is not going to satisfy your requirement to compute the singular value decomposition without using svd. Let's look at one of the lines in the implementation of the orth function.

dbtype 19 orth.m

19    [Q,s] = svd(A,'econ','vector');

[John D'Errico]

It makes a strange sense to me, in a mad, deluded way of course. Implement SVD using ORTH. Implement ORTH using SVD. You never need to write any actual code at all then. Infinite recurrences? What, me worry?

[Caleb]

I was attempting to do svd without the built in function as a test for myself to make sure I fully understood the theory behind it. However, I have found that may have been a bigger task then I was expecting.

[David Goodmanson]

Hi Caleb,

It's all to your credit to make your own svd, circular reasoning or not, if it helps to better understand what the svd is up to.

If the goal were to replace Matlab svd with your own version from here on out, that's of course not a good idea since you would be throwing away years of work on making the function as numerically accurate as possible. But it does not sound like that is what you have in mind. I'm not attempting here to 'explain' the svd but the code below does provide a numerical result (the svd not being unique in general).

Starting truly from scratch is a lengthy process and it makes sense to have Matlab do some of the work so I am assuming that you have a pxq matrix A with p>q and have used the idea that the singular values of A are the square roots of the eigenvalues of A'*A. Also I am assuming no repeated eigenvalues since that complicates things. So

p = 7;
q = 4;
A = randc(p,q)
[v lam] = eig(A'*A);
% find s1, the diagonal qxq part of s, sort the values, sort cols of v
% to match, use A = u*s1*v' to obtain the first q columns of u.
% the last, p-q columns of u don't count because they multiply the
% lower, all-zeros part of s
[sqrtlam ind] = sort(diag(sqrt(lam)),'descend');
s1 = diag(sqrtlam);
v = v(:,ind)
u = A*v/s1;
s = [s1;zeros(p-q,q)]

Now you need the last p-q columns of u, as you mentioned. These span the p-q vector space but are otherwise arbitrary. You can obtain those at once with

unull = null(u');
u = [u unull];

but you can do more of the work by means of gram_schmidt orthogonalization

for k = 1:p-q
  ucol = randc(p,1);
  ucol = ucol - u*(u'*ucol);
  ucol = ucol/norm(ucol);
  u = [u ucol];
end
% check
max(max(abs(A - u*s*v')))
unitry(u)
        % ans =  9.5505e-16
        % ans =  1.4433e-15

function y = randc(varargin)
% uniform random numbers on the complex region [-1, 1] x [-i ,i]
y = 2*rand(varargin{:}) + i*2*rand(varargin{:}) - (1+i);
end

function delta = unitry(u)
% find the amount by which u is short of being unitary.
% more foolproof than a yes/no result vs some specified tolerance. 
I = eye(size(u));
delta = max(max([abs(u'*u - I) abs(u*u' - I)])); 
end

[Paul]

Hi David,

Regarding "Also I am assuming no repeated eigenvalues since that complicates things."

Which matrix is assumed to have no repeated eigenvalues and why would that complicate things?

[David Goodmanson]

Hi Paul,

I was referring to the eigenvalues in

[v lam] = eig(A'*A);

As you know, if an eigenvalue is repeated there is not a unique eigenvector for each eigenvalue. (Here an eigenvector multiplied an aribitrary phase factor is still assumed to be the same unique eigenvector.) The eigenvectors just have to span the space of dimension corresponding to the degeneracy of the eigenvalue. So even if the columns of v are sorted by eigenvalue size, v is not unique. That means that the step

u = A*v/s;

has some arbitraryness to it. Maybe that step works out all right anyway, but at the time I just didn't want to think about it. So I added the caveat about complication. Since the svd is not unique in any case, does it seem to you that the step above will work regardless?

[Paul]

Hi David,

IIRC, the eigenvalues of a symmetric matrix are never defective and eigenvectors are an n-dimensional basis.

[Christine Tobler]

Yes, the eigenvectors of a symmetric matrix form an orthonormal basis, so there should be no problem with multiple eigenvalues.

Problems are likely if there are very small or zero singular values (dividing by zero with s1, and the possibility of eig(A'*A) returning very small negative eigenvalues).

Not to detract from the explanation, which is useful to describe the connection between eig and svd. This just means in practice, it's better to call svd directly.

[David Goodmanson]

I didn't doubt that A'*A has a full set of eigenvectors and is non-defective, as both of you have pointed out. It's just that when there are no repeated eigenvalues, each eigenvector is pinned down and well defined. When an eigenvalue is repeated q times, the associated eigenvectors are free to roam around in a vector space of dimension q. It probably works just fine anyway, but in the step

u = A*v/s;

I was/am not sure that the calculation of u is good every time.

[Paul]

Hi David,

According to the doc page at eig - Eigenvalues and eigenvectors - MATLAB, "If A is real symmetric, Hermitian, or skew-Hermitian, then the right eigenvectors V [returned by eig] are orthonormal."

Does that assuage concerns about v?

I added the part in brackets because, it seems to me, that the eigenvectors associated with repeated eigenvalues of a symmetric matrix don't have to be orthonormal, though it looks like eig does enforce that for the single input case, e.g.,

Z = zeros(2);
[V,D] = eig(Z)

V = 2×2

1 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 2×2

0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[V,D] = eig(eye(2))

V = 2×2

1 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 2×2

1 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I am curious as to how eig() works for symmetric matrices. Does it use issymmetric (or equivalent) on the input and take different paths depending on the result? Also, is M = A'*A guaranteed to always be symmetric in floating point (I thought that Matlab will magically ensure that M is symmetric, but I can't find anything on point in the doc).

The two-input form of eig is more interesting (I readily confess to not being that familiar with the generalized eigenvalue problem)

[V,D] = eig(Z,Z,'chol')

V = 2×2

1 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 2×2

NaN 0 0 NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[V,D] = eig(Z,Z,'qz')

V = 2×2

1 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 2×2

NaN 0 0 NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I don't know what to make of the D-matrix, insofar as any finite D and any finite, non-zero V will satisfy Z*V = Z*V*D

[David Goodmanson]

Hi Paul, as I said I didn't doubt that Matlab eig (with one input) produces an orthonormal matrix for V, whether the eigenvalues are repeated or not. I wasn't positive if in the repeated case, the resulting u from

u = A*v/s

would always be correct. But now I have verified it both by proof and example. So the initial caveat I had is unnecessary.

In what you posted above, is Z still hanging around as zeros(2) from the first line? Then it will be hard to construct anything out of eig(Z,Z).

Answer 1

[John D'Errico]

UM, You seriously do not want to use X'*X to compute anything about the SVD. Yes, mathematically, it is the same. But this is not the case computationally!!!!!!

How does MATLAB compute the singular vectors? It probably calls code from LAPACK, as I recall. When I wrote my own SVD code a million years ago, based on the old LINPACK algorithms, the idea was to bidiagonalize the matrix using Householder transformations. Then I recall I needed to do Givens rotations to finish things off, killing off the off-diagonal elements one at a time. (Its been a long time since I needed an SVD, so that is just my fading memory.) You can find the basic algorithm I used back then in the original LINPACK manual.

The point is, at no time do you EVER want to form X'*X. The problem is, once you do that, you destroy any information available down in the lower significant bits of your array. And this is why, when you tried to do it, those small singular values (and their associated singular vectors) are complete crap. Worthless.

[James Tursa]

E.g., some LAPACK documentation for SVD computation routines can be found on netlib:

https://www.netlib.org/lapack/explore-html/d1/d90/group__svd__top.html

Answer 2

[Christine Tobler]

Hi Caleb,

Happy to discuss the svd with you!

That being said, as others here have mentioned, svd is a building block and not easy to implement. Things like orth, rank, null, pinv all are built on top of the SVD, and so wouldn't be an independent way of computing the SVD.

There are some ways to write up the SVD through symmetricEIG, but then EIG is similarly a building block and hard to implement - so it's a bit like going in circles.

What is your final goal? Would you like to understand how the SVD is computed? Or are you looking at a specific problem where you are currently using the SVD?

[Caleb]

I was attempting to ensure that I understood the singular value theory and what each component was. However, as I now see from your comment and others, this might be a bigger task than I initially assumed.