How to compare the outputs of two algorithms computing SVD?

Question

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

Answer 1

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] and [1] 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.
• [1] Åke Björck and Gene H. Golub (1973). Numerical methods for computing angles between linear subspaces. Mathematics of Computation, 1973, Vol.27, pp.579-594.

Answer 2

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$.