• singular values distance = squared L2 norm of the difference between singular values from standard SVD and approximated ones from randomized_svd
• Frobenius distance = Frobenius norm of the difference between original input matrix and matrix product USV, output of randomized_svd
• noise component = parameter in [0,1] input of `make_low_rank_matrix
• the size of the dots in the last plots is proportional to the number of power iterations used for that run --growing to the right

So far, plots suggest 2 may be the choice. Larger number may even degrade quality of the approximation. The flat behaviour of the last two plots come from very small datasets.

Ping @ogrisel

This looks great although the plots on diabetes and abalone are not very interesting as they extract a single component. Maybe you should force to extract at least 5 components to make them more relevant.

scipy.sparse.linalg.svds does decrease a little time consumption, but the singular values and singular vectors returned are sorted differently from randomized_SVD and scipy.linald.svd. I would stick with the original, there is not much gain here anyway.

I have coded a last test. Looking to the 2 new plots, I think we can be confident that 2 is a good default value. Input was (500, 5000). Rank is the rank of the input matrix, n_comps is the parameter of randomized_svd.

For completeness, here the other plots updated. (There is no change on the small datasets, even computing 5 components.)

This is great. Thanks for this empirical study. +1 for using 2 power iterations in randomized_svd by default.

I wonder if we should change the default value for RandomizedPCA and TruncatedSVD from 3 to 2 as well:

• on the plus side: this would give a non-negligible speed up with similar accuracy on most datasets,
• on the negative side: it can cause a slight change of behavior that users upgrading from a previous version of scikit-learn might not expect.

Any opinion by others, e.g. @larsmans?

How do we compare against fbpca? https://github.com/facebook/fbpca ? That also just uses power iterations, right?

Another thing which would be nice to investigate if you have time: our implementation of range_finder uses Algorithm 4.3 but according to the paper the procedure is numerically unstable. The paper suggests a subspace iteration method in Algorithm 4.4. It would be interesting to compare both algorithms in terms of accuracy and computational time.

How do we compare against fbpca? https://github.com/facebook/fbpca ? That also just uses power iterations, right?

It does. And it also uses 2 power iterations by default.

I ran a few tests. Sklearn seems slower. Performance is similar, but sklearn is slightly better here. A difference is that fbpca's performance does not deteriorate when the number of power iterations is "too" large. (I did not looked into the code.)

Code on this gist. Graphs attached.

I wonder if we should change the default value for RandomizedPCA and TruncatedSVD from 3 to 2 as well:

What people think about this? @ogrisel @larsmans I am happy to have a look into it.

Sklearn seems slower. Performance is similar, but sklearn is slightly better here.

I don't understand this, could you please rephrase?

Also could you please annotate the scatter plot with the number of power iterations for each dot?

See for instance: http://stackoverflow.com/questions/5147112/matplotlib-how-to-put-individual-tags-for-a-scatter-plot (no need for arrows)

Just commenting on the two plots. Sklearn runs slower than fbpca. But performance in term of matrix factorization error is slighly better than fbpca --until sklearn starts to degrade, for lfw_people.

Also could you please annotate the scatter plot with the number of power iterations for each dot?

Sure.

The difference is how the power iterations are computed, see:

https://github.com/facebook/fbpca/blob/master/fbpca.py#L1562
vs
https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/utils/extmath.py#L226

A LU factorization lu(Q, permute_l=True) before each power iteration.

It would be worth checking the paper: maybe our code does not do it because of an oversight on our part. It seems to make the power iteration more stable: too many power iterations do not seem to cause a degradation in the quality of the approximation when the LU steps are there.

Maybe @ajtulloch would be interested in following this discussion :)

Plots updates

I think fbpca switches "modes" depending on the number of components - IIRC
from looking at the code they switch to doing "full" pca and truncating if
the number of components is >X% of the total. Are we comparing only to the
randomized solver of fbpca or to the external interface (which may switch
modes internally)?

On Mon, Aug 31, 2015 at 10:32 AM, Giorgio Patrini notifications@github.com
wrote:

Plots updates
[image: figure_1]
https://cloud.githubusercontent.com/assets/2871319/9581136/d6fa1710-4ffd-11e5-89da-de7cef990135.png
[image: figure_2]
https://cloud.githubusercontent.com/assets/2871319/9581137/d6fe918c-4ffd-11e5-804e-94f612e59d16.png

—
Reply to this email directly or view it on GitHub
#5141 (comment)
.

and should we also switch modes internally?

I think fbpca switches "modes" depending on the number of components - IIRC from looking at the code they switch to doing "full" pca and truncating if the number of components is >X% of the total. Are we comparing only to the randomized solver of fbpca or to the external interface (which may switch modes internally)?

This is only the case when n_components is very close min(n_samples, n_features). We could implement that strategy in sklearn as well (as a safety belt) but in practice I don't think it matters as on large enough data you can only do PCA with n_components << min(n_samples, n_features) (otherwise the SVD is too expensive).

Looking into the code, the 3 differences I see are:

• mode switching when n_comp > 4/5 min(n_samples, n_features).

Are we comparing only to the randomized solver of fbpca or to the external interface (which may switch modes internally)?

I am comparing with the external interface pca but I do not think there is any difference, except passing through n_comp > 4/5 min(n_samples, n_features). This was always false in what I ran.

• a LU factorization between every multiplication of A and A^T in the power iterations. This may be the reason why sklearn runs slightly faster. And this is very likely the reason of the increasing error after many power iterations, replying to @mblondel. From Halko et al.:

Unfortunately, when Algorithm 4.3 is executed in floating point arithmetic, rounding errors will extinguish all information pertaining to singular modes associated with singular values that are small compared with A . (Roughly, if machine precision is μ, then all information associated with singular values smaller
than μ 1/(2q+1) A is lost.)

The paper suggests QR factorizations anyway.

• fbpca does not factor the cases m>n and m<n. I guess the support to complex input does not make it so straightforward, as it is in sklearn that does .T at the beginning.

This may be the reason why sklearn runs slightly faster.

Actually on LFW, sklearn is slightly slower, not faster. Maybe the final QR is less expensive when the LU is done in the intermediate steps? It would be great to confirm with an experiment by adding an option to add the LU steps in our codebase.

And this is very likely the reason of the increasing error after many power iterations.

This looks very likely indeed.

The implementation in dask.linalg by @marianotepper does vanilla power iterations (no LU):

https://github.com/blaze/dask/blob/master/dask/array/linalg.py#L225

It would be interesting to compare QR and LU decompositions.

Another thing I am curious is whether it would be worth using a stopping criterion for the loop in the range finder. Halko et al. suggest ||A - QQ^TA||. This can however be expensive to compute.

Cool, great discussion! LMK if we (Mark Tygert and I wrote fbpca) can contribute in any way. There's a fairly detailed evaluation at http://tygert.com/implement.pdf if that's useful as well.

Thanks for the very useful reference @ajtulloch! I am going to make a benchmark file for playing with different versions of the power iterations. I will keep comparing with your fbpca.

whats_new.rst updated.

lgtm.

This lgtm as well

The what's new entries should be put under 0.18.

The what's new entries should be put under 0.18.

Sure. I was hoping for a merge in 0.17 ;)
I will update. I am also working on speeding up the benchmarks with the computation of approximate spectral norms as discussed above.

Also, we should have a common default policy for the power iterations (n_iter and n_oversampling) with the other PCA open PR ##5299. Don't we?

whats_new updated, with some other minors.

I did not implement the faster spectral norm estimate. That measure depends on random initialization, which make benchmarking trickier (need many runs and averaged performance etc.).

Done here.

+1 for merging this now and re-adjusting the default params for consistency in the PCA collapse PR if needed.

🍋cello!

🍻

:craft beer: