# Making MATLAB svd robust to transpose operation

I'm playing with MATLAB's svd function to compute the svd of

 [ 1     4     7    10
2     5     8    11
3     6     9    12 ]


When I type [U1, ~, ~] = svd(X), I get

U1 =

-0.5045    0.7608    0.4082
-0.5745    0.0571   -0.8165
-0.6445   -0.6465    0.4082


But when I compute the svd of the transpose of X with [~, ~, U2] = svd(X'), I get

U2 =

0.5045    0.7608    0.4082
0.5745    0.0571   -0.8165
0.6445   -0.6465    0.4082


The first singular vectors seem to be pointing to the opposite directions but the others are the same. I know that svd is not unique and the solution is correct because the first component of V1 and V2 are pointing opposite directions as well. But, I would expect MATLAB to return the same answers. I thought to add a postprocessing and check the singular vector pairs and turn their directions to make them consistent under transpose but I couldn't find a mathematically reasonable way of doing it.

Do you know why MATLAB (or LAPACK as it says MATLAB uses LAPACK) computes this way? Do you have any suggestions how to make it consistent? Thanks.

Update: the reason I ask for it is that I wanted to see what happens if I apply higher order singular value decomposition (HOSVD) to 2 dimensional matrices. In theory, it holds. Using the same mathematical notation, SVD is formulated as follows:

Every $\mathbf{X} \in \mathbb{C}^{I_1 \times I_2}$ can be written as the product

.$$\mathbf{X} = \mathbf{U}^{(1)} \cdot \mathbf{S} \cdot \mathbf{U}^{(2)^H} = \mathbf{S} \times_1 \mathbf{U}^{(1)} \times_2 \mathbf{U}^{(2)}$$

in which

$\mathbf{U}^{(1)} = \left[ \begin{array}{cccc} U_{1}^{(1)} & U_{2}^{(1)} & \dots & U_{I_1}^{(1)} \end{array} \right] \in \mathbb{C}^{I_1\times I_1}$ is unitary.

$\mathbf{U}^{(2)} = \left[ \begin{array}{cccc} U_{1}^{(2)} & U_{2}^{(2)} & \dots & U_{I_I}^{(2)} \end{array} \right] \in \mathbb{C}^{I_2\times I_2}$ is unitary.

$\mathbf{S} \in \mathbb{C}^{I_1\times I_2}$ has the the following properties:

• pseudodiagonality: $\newcommand{\diag}{\mathop{\mathrm{diag}}} \mathbf{S} = \diag\left(\ \sigma_1, \sigma_2, \dots, \sigma_{\min(I_1,I_2)} \right)$

• ordering: $\sigma_1 \geq \sigma_2 \geq \ldots \geq \sigma_{\min(I_1, I_2)} \geq 0$

So, to compute it I found the right and singular vectors as if they were computed independently because I want to extend it to multidimensions later. To compute them independently, I calculated the eigenvalues of $\mathbf{XX'}$ and $\mathbf{X'X}$ as follows:

[U{1}, ~] = eig(X*X');
U{1} = fliplr(U{1});
[U{2}, ~] = eig(X'*X);
U{2} = fliplr(U{2});


but they seem to be as I explained above. Actually, it doesn't effect what I need but I would like to make the singular values positive. This way, they might become negative since the right and left singular vectors are not as I expected.

Update 2: Alternatively, I used the code with svd() function as below and it gives similar results. With the 'econ' option, only min(nDim1,nDim2) singular vectors are returned and I set the remaining as 0 vector. The economy option (i.e. computing the non-zero vectors makes it practical when the higher order matrices [n-way arrays] are turned into matrices and svd is applied to them.)

[U{1}, ~, ~] = svd(X, 'econ');
U{1} = [U{1} zeros(size(U{1},1), size(U{1},1)-size(U{1},2))];
[U{2}, ~, ~] = svd(X', 'econ');
U{2} = [U{2} zeros(size(U{2},1), size(U{2},1)-size(U{2},2))];

-
Have you already tried inspecting the other matrix of singular vectors and looking at the sign pattern? – J. M. Oct 11 '10 at 8:30
Yes, I've. It seems to show a similar pattern. The first singular vectors are in opposite signs but the second and the third ones are the same as expected. – İsmail Arı Oct 11 '10 at 12:29
There is inherent sign ambiguity in the singular vectors. There is even more ambiguity if there are repeated singular values (think of the identity matrix as an extreme case). – Darsh Ranjan Oct 14 '10 at 1:11
What do you need this normalization for? Please tell us your complete problem, maybe it is something that is better achieved at another step. – Federico Poloni Oct 14 '10 at 9:15
@Darsh Yes, you are definitely right. I don't know how it is computed in the background but at least for this one, it would be great to have the same results. @Federico I added the reason in the question. – İsmail Arı Oct 14 '10 at 10:25

You may be interested in the following article which addresses this issue:

R. Bro, E. Acar and T. G. Kolda. Resolving the sign ambiguity in the singular value decomposition. Journal of Chemometrics 22(2):135-140, February 2008. (doi:10.1002/cem.1122)

http://www.mathworks.com/matlabcentral/fileexchange/22118-sign-correction-in-svd-and-pca

Best wishes, Tammy Kolda

-
I have checked most of your papers but missed this one. Thank you for your kind response. – İsmail Arı Nov 17 '10 at 9:31

In Matlab, there is a sign ambiguity (i.e. $\pm1$). In LAPACK, it is worse, the ambiguity is $e^{i \phi}$, which Matlab removes later on.

Maybe you can take each vector of $U$ individual, find the coefficient with the largest absolute value, and make sure that this coefficient is always positive.

This problem has never bothered me before though.

-
When I find the largest valued vector, they I think your suggestion is to multiply the directions with -1. But again the vectors do not fit as in the example solution in the question. Do you have any other suggestion to do after finding this large valued vector. – İsmail Arı Oct 14 '10 at 10:22
@Ismail: did you remember to sort the eigenvectors returned by eig()? Just to be sure, you never mentioned it... – Federico Poloni Oct 14 '10 at 11:07
Sorry, I forgot to mention. Yes, I did. I updated the code in the question. – İsmail Arı Oct 14 '10 at 13:08

If you compute the (HO)SVD by solving the two eigenproblems, you always need to do some post-processing to sort the singular values and make them positive.

I'd suggest moving the normalization issue to the comparison operator: to compare your results with the SVD (i.e., that $U=\tilde U$, where the two are the same matrix computed with two different methods, up to unitary diagonal factors), a possibility could be computing the scalar products $\left\langle U(:,k),\tilde U(:,k)\right\rangle$ and checking if their absolute value is $1(\pm \epsilon)$.

-
Actually, the scalar products always be +/-1. I couldn't understand clearly how to do a normalization on comparison operator. The trick is to not to compute all svd but solve it partially. I tried and added another piece of code into the question where I used svd() function with economic usage, where it finds only min(M,N) singular vectors. The significant vectors again do not fit in sign. – İsmail Arı Oct 14 '10 at 13:18