# How to compare the outputs of two algorithms computing SVD?

For example we have 2 algorithms from R: SVD and irlba and I want to compare them int terms of speed,memory and precision.

But I don't understand how to compare output of algorithms, they must be differ on some epsilon, but how to define epsilon?

And why sign of principal vectors can be different?

Here is some test code:

rows= 10
cols= 8
A <- matrix(runif(rows*cols),rows,cols)

pc=5
a <- irlba(A,nu=pc,nv=pc)
b <- svd(A,nu=pc,nv=pc)

#false
a$u == b$u
a$v == b$v

#show all vars
a$u b$u
a$v b$v
a$d b$d

#compare u
.list <- list(abs(a$u), abs(b$u))
Reduce('-', .list)

#compare v
.list <- list(abs(a$v), abs(b$v))
Reduce('-', .list)

#compare d
.list <- list(abs(a$d), abs(b$d[c(1:pc)]))
Reduce('-', .list)


output:

> #compare u
> .list <- list(abs(a$u), abs(b$u))
> Reduce('-', .list)
[,1]          [,2]          [,3]          [,4]         [,5]
[1,]  6.383782e-16 -5.551115e-17 -1.082467e-15  1.665335e-15 -0.123009729
[2,]  0.000000e+00  5.689893e-16 -3.053113e-16 -1.665335e-15 -0.117362383
[3,]  1.110223e-16 -5.481726e-16 -2.220446e-16  7.216450e-16  0.142286325
[4,]  3.885781e-16 -3.053113e-16 -1.526557e-16 -2.942091e-15  0.096784035
[5,]  4.996004e-16 -9.436896e-16  1.776357e-15 -1.956768e-15 -0.081504771
[6,] -2.775558e-16 -6.661338e-16  7.216450e-16  1.137979e-15  0.007450404
[7,] -5.551115e-16  4.440892e-16 -5.273559e-16  5.551115e-17  0.018699113
[8,] -5.551115e-16  7.216450e-16  4.440892e-16  1.665335e-16 -0.013035052
[9,] -4.440892e-16  1.054712e-15 -1.942890e-16  2.498002e-16  0.048095160
[10,]  3.330669e-16  3.330669e-16  7.771561e-16  6.106227e-16  0.044937739
>
> #compare v
> .list <- list(abs(a$v), abs(b$v))
> Reduce('-', .list)
[,1]          [,2]          [,3]          [,4]        [,5]
[1,] -1.887379e-15 -2.775558e-16 -9.020562e-17 -1.193490e-15 -0.02102106
[2,]  1.831868e-15 -1.276756e-15 -6.661338e-16  3.386180e-15 -0.09872305
[3,]  7.216450e-16 -3.885781e-15 -9.436896e-16  9.992007e-16  0.04699264
[4,] -2.053913e-15  2.331468e-15 -4.996004e-16  1.498801e-15  0.05876057
[5,]  8.881784e-16 -5.551115e-16 -8.881784e-16 -1.942890e-16  0.16146650
[6,] -2.053913e-15  8.881784e-16  5.551115e-16 -2.775558e-16  0.11420391
[7,]  1.498801e-15 -1.665335e-16  1.665335e-16  9.714451e-16 -0.13600871
[8,]  8.326673e-16 -2.831069e-15 -7.042977e-16 -2.997602e-15 -0.01507590
>
> #compare d
> .list <- list(abs(a$d), abs(b$d[c(1:pc)]))
> Reduce('-', .list)
[1] -8.881784e-16  0.000000e+00  8.881784e-16  2.220446e-16 -2.483630e-02


here is also some comparision https://github.com/graphlab-code/graphlab/issues/48

-
Often it's more meaningful to compare the outputs in terms of accuracy rather than precision. To that end, have you considered providing the algorithms various inputs for which you know the answers exactly? (It's pretty easy to find such inputs for SVD: begin with the desired output $U, \Sigma, V$ and use those to calculate what the input matrix must be.) The thread at stats.stackexchange.com/questions/34396 addresses the question about the signs of the components--and shows how and why they cannot even be consistent within the same software implementation. –  whuber Sep 16 at 20:12

Any SVD algorithm outputs a space $\pmb V$ of orthogonal vectors. Your best bet is in comparing those sub-spaces.

Consider the following measure. Let $\pmb I_m$ be the $m$ dimensional subspace returned by one of the SVD algorithm and $\pmb P_n$ the $n$ dimensional subspace returned by the other (in R, these are stored in the \$v object of the output). Then, with$k=\min(m,n)$, denote$\pmb M_k$the matrix: $$\pmb M_k=\pmb I_k'\pmb P_k\pmb P_k'\pmb I_k$$ where$\pmb I_k$is the sub-matrix formed of the first$k\leqslant m$columns of$\pmb I_m$(and similarly for$\pmb P_k$and$\pmb P_n$).$\pmb M_k$is$k$by$k$, symmetric and positive definite. Then, consider the scalar$m$: $$m=\frac{2}{\pi}\arccos(\sqrt{\lambda_k})$$ where$\lambda_k$is the the smallest eigenvalue of$\pmb M_k$. This measure is called the maxsub. The maxsub actually measures the (absolute value of the) largest angle between a vector in$\pmb I_k$and the vector most parallel to it that belongs to$\pmb P_k$. The factor$\frac{2}{\pi}$is just there normalize the measure to lie in$[0,1]$(see also [0] for more details). Contrary to the measure you propose (difference of the two matrix of singular vectors), the maxsub measure is invariant to orthogonal transformations (e.g. rotations) of the data, which is a natural property if the two algorithms you use to compute the orthogonal vectors also are (this is certainly true for the SVD algorithm, I do not know for irlba). • [0] Krzanowski, W. J. (1979). Between-Groups Comparison of Principal Components. Journal of the American Statistical Association, Vol. 74, No. 367, pp. 703-707. - I really like @user603's answer but a more quick-and-dirty method might also suit your needs. This approach will work fine for smaller dimensions but if you have to reconstruct a large matrix it will be computationally expensive. It's possible to just reconstruct your$A$matrix and look at the Frobenius norm between your original matrix and the reconstructed one. The Frobenius norm can be thought of as a distance between two matrices so the idea is that, under the accuracy criteria, you'd pick the method with the smallest Frobenius norm. Here's some R code that continues your example: # reconstruction of matrices recon_a = a$u %*% diag(a$d) %*% t(a$v)
recon_b = b$u %*% diag(b$d[1:pc]) %*% t(b$v) # calculate F norm for irlba reconstruction .list <- list(A, recon_a) norm(Reduce('-', .list), type = "F") # calculate F norm for svd reconstruction .list <- list(A, recon_b) norm(Reduce('-', .list), type = "F")  One note: if I'm not mistaken the SVD gives the best$k$-rank approximation - in terms of the Frobenius norm - to the matrix$A$so it should always win using this evaluation method. But since you have several other criteria (speed and memory) where irlba should be winning this will give you an idea of how badly irlba is losing in precision over a range of$k\$.

-
I also was thinking about reconstruction accuracy, but why you call it quick-and-dirty? Also I choosing not the more precise algorithm, but I compare "ideal" canonical slow implementation with more fast implementation which use some approximations, so it's accuracy will be less.It can be more or less accurate, but how to measure "degree of accuracy" and determine if it appropriate for real tasks?(I don't know how it's called maybe numerical stable?) –  mrgloom Jul 25 at 6:03
the main problem with this solution is that it is not rotation invariant. If irlba is rotation equivariant (and we know SVD is), then, this solution will give inconsistent results (you could construct two matrices that are equivalent up to a rotation but they won t yield the same ranking). –  user603 Sep 16 at 18:04