How does the svd function determine or fix the phase of singular vectors?

Hello everyone,

I have been experimenting with MATLAB's `svd()` and `eig()` functions to compute the singular value decomposition of complex matrices. I noticed that:

- The results from `svd()` and `eig()` are numerically very close, but not exactly identical.

- Sometimes the signs of the singular vectors match, while other times they differ by -1 (or a phase factor in the complex case).

I would like to understand:

1. How does MATLAB internally compute SVD?

2. When there is phase or sign ambiguity in the singular vectors, how does `svd()` standardize or normalize them?

Any insights or references on this would be greatly appreciated, as I want to use this knowledge in my research on tensor analysis and array signal processing.

Thank you very much!

Best regards,

Lv-lin Zhong

Answer 1

[Steve Eddins]

For a discussion of how MATLAB computes the SVD, see Cleve Moler's blog post, "Two Flavors of SVD," 23-Feb-2025.

The svd doc page implies that singular vectors are not "standardized." See the descriptions of the output arguments. For example, the description of U says: "Different machines and releases of MATLAB® can produce different singular vectors that are still numerically accurate. Corresponding columns in U and V can flip their signs, since this does not affect the value of the expression A = U*S*V'."

Answer 2

[John D'Errico]

First, signs are to some extent arbitrary. That is, for EIG, you can multiply an eigenvector by -1, and it is still an eigenvector. Nothing changes. For SVD, it is not so easy, but you can still change the signs of the vectors to some extent, but you will need to then also fiddle with the signs of the other singular vectors.

That is, if we have

A = U*S*V'

then trivially

A = (-U)*S*(-V')

Can you flip signs of a singular vector? Well, yes, if you are careful. Consider this example:

A = rand(5)

A = 5×5

0.4047 0.0550 0.4328 0.3928 0.4381 0.8994 0.3817 0.2035 0.8605 0.6195 0.8756 0.0026 0.9248 0.1412 0.2320 0.6362 0.9016 0.7670 0.5431 0.6451 0.7280 0.2173 0.1320 0.6447 0.1886

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

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

[U,S,V] = svd(A)

U = 5×5

-0.3045 -0.0928 0.0869 0.8064 0.4908 -0.5297 0.4942 0.2681 0.1957 -0.6042 -0.4152 -0.8055 0.3349 -0.1642 -0.1994 -0.5714 -0.0194 -0.7978 -0.1769 0.0738 -0.3576 0.3131 0.4147 -0.5032 0.5906

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

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

S = 5×5

2.6033 0 0 0 0 0 0.8526 0 0 0 0 0 0.6881 0 0 0 0 0 0.3028 0 0 0 0 0 0.0996

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

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

V = 5×5

-0.6096 -0.0970 0.5290 -0.3976 -0.4255 -0.3122 0.2720 -0.7574 -0.4963 -0.0926 -0.4260 -0.7718 -0.2256 0.1149 0.3984 -0.4513 0.5471 0.2125 0.1368 0.6581 -0.3818 0.1469 -0.2247 0.7508 -0.4674

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

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

As you can see, U,S,V can be used to reconstruct A.

U*S*V' - A

ans = 5×5

1.0e-15 * 0 0.4025 0.8327 0.6106 0.6661 0.1110 -0.2220 0.1943 0.3331 0.1110 -0.6661 0.6926 0.3331 0.1388 0.1388 -0.2220 0.3331 0.1110 0.4441 0.2220 0 -0.1388 0.1110 0.2220 0.0555

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

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

Now tweak the signs of the second singular vectors by -1.

U2 = U;U2(:,2) = -U2(:,2);
V2 = V;V2(:,2) = -V2(:,2);
U2*S*V2' - A

ans = 5×5

1.0e-15 * 0 0.4025 0.8327 0.6106 0.6661 0.1110 -0.2220 0.1943 0.3331 0.1110 -0.6661 0.6926 0.3331 0.1388 0.1388 -0.2220 0.3331 0.1110 0.4441 0.2220 0 -0.1388 0.1110 0.2220 0.0555

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

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

And the result is the same, all floating point trash. In fact, the trash is identical in both cases.

Next, you tell us the results are close, but not exact. SURPRISE! They are done using two different algorithms, and anytime you have a set of numbers computed using two distinct sequences of operations or by compeltely different algorithms, you should NEVER expect to see identical results. Must I repeat that? NEVER. In fact, it is more often a surprise if you do see the same results down to the least significant bit.

A trivial example of this is...

x = 0.3 - 0.2 - 0.1

x = -2.7756e-17

y = -0.2 - 0.1 + 0.3

y = -5.5511e-17

x and y should be identical, no? Even for such a simple thing, they are not, because these computations were done in double precision, not in infinite precision arithmetic.

How is SVD computed? Its been a long time since I wrote code for an SVD, but I recall the old LINPACK codes for the SVD first b-diagonalized the matrix, using Householder transformations. Then they killed off the upper diagonal elements using a sequence of Givens rotations. I'd guess that is not how the corrent code works since it has been roughly 40 years since I wrote code for an SVD, but that description is probably not terribly far off.

Finally, how does SVD standardize/normalize the eigenvectors? It does nothing of the sort. There is no need to artifically normalize them, since they are already unit vectors. That is, every operation done to create them is in the form of a product with a "rotation" matrix, and that does not change the norm of a vector.

If you want to understand more than that, my guess is MATLAB uses LAPACK codes to compute both EIG and SVD. So start reading the LAPACK documentation for those tools. The docs for EIG and SVD say a little, but they don't really state the methods used in any depth, since those algorithms might change in a future release.

Answer 3

[Christine Tobler]

Hi Lv-lin,

1.

MATLAB computes SVD internally through a completely separate algorithm than EIG, with both algorithms provided by the LAPACK library.

While in exact arithmetic, the svd can be computed through eig on A'*A, this would result in larger numerical errors. The algorithms for EIG and SVD are related, but not identical.

2. MATLAB applies no standardization to the singular vectors and eigenvectors in terms of phase or sign. In the svd, and for simple Hermitian eigenvalue problems, the eigenvectors are normalized in Euclidean norm and form an orthogonal basis.

That said, for the complex case, the LAPACK method applies some standardization for the phase - we don't document this, because in future there might be new, faster methods that choose a different standardization. My recommendation is that if you require some standardizeation of the phase, this should be applied to the outputs, to guarantee it will stay the same over time.

In complex non-Hermitian eig, the phase is chosen so that the element with largest magnitude for each eigenvector is real positive:

[U, ~] = eig(randn(4, 'like', 1i))
max(abs(U))

I find this a quite reasonable standard, which you can apply to U for all complexities and "symmetricities" as follows:

% Complex Hermitian has a different standardization
A = randn(4, 'like', 1i); A = A+A';

[U, ~] = eig(A)

% Find element with maximum absolute value
maxAbs = max(U, [], ComparisonMethod="abs")

% Multiply with the conjugate of the sign of this element
U = U .* conj(sign(maxAbs))

% Verify that now, the real positive element in each column has maximal
% magnitude
max(abs(U))

For SVD, you could choose either U and V, and apply a similar standardization. What standardization you apply (if any) depends on your use case. For many applications, no standardization is needed at all.