[MRG] Bugfix for not-yet-confirmed issue #12863: arpack returns singular values in ascending order, the opposite was supposed in sklearn

[gykovacs]

Bugfix for not-yet-confirmed issue #12863: arpack returns singular values in ascending order, the opposite was supposed in sklearn

Reference Issues/PRs

Fixes [MRG] not-yet-confirmed issue #12863

What does this implement/fix? Explain your changes.

In the sklearn code it was assumed that the scipy wrapper for arpack singular value decomposition returns singular values (and corresponding vectors) in descending order. This is not the case. The ARPACK documentation (https://www.caam.rice.edu/software/ARPACK/UG/node136.html#SECTION001210000000000000000
) clearly states that the singular values are returned in ascending order.

The bug manifests in differing results depending on the solver:

import numpy as np
from sklearn import cluster

A= np.array([[-2, -4, 2], [-2, 1, 2], [4, 2, 5]])

sc= cluster.bicluster.SpectralCoclustering(n_clusters= 2, 
                                           svd_method='randomized')

sc.fit(A)
print(sc.column_labels_)

gives

[0 0 1]

, but

sc= cluster.bicluster.SpectralCoclustering(n_clusters= 2, 
                                           svd_method='arpack')
sc.fit(A)
print(sc.column_labels_)

gives

[0 0 0]

This bugfix makes the arpack call return the same result as the randomized sv based solution.

Any other comments?

[jnothman]

Please add a test.

[gykovacs]

I'm working on it!

[gykovacs]

I have added a test checking if the labels implied by the singular values computed by 'arpack' and 'randomized' methods are equal (the svs are not available in the BaseSpectral objects, so the direct comparison of the order of svs cannot be tested).

I also tested the test, it passes with my updates to BaseSpectral._svd but fails without it, so the test matrix seems to be a good choice to drive the testing.

However, one CI job fails: the one related to OSX. Interestingly, another svd calculation fails, the numbers do not match the expectations in the 12th digit. In my best understanding this has nothing to do with my bugfix, this mismatch is based on an svd call made directly to scipy at line 163 of decomposition/truncated_svd.py. I'm wondering how to proceed. Shall I maybe decrease the required precision in this test to 11 digits?

[jnothman]

I've not looked into the failure, but I think it must be due to your change, even if your change is right.

[jnothman]

Yes, I think it should be okay to increase the tolerance.

[jnothman]

(The rpca in that test code is a bit weird... both are using 'arpack', not 'randomized')

[gykovacs]

I did some experimentation, let me summarize the results:

1. the test I added (in test_bicluster) makes a completely different test (test_singuler_values in test_truncated_svd.py) fail.

2. the common in the two tests is the arpack call

3. seemingly, the scipy arpack based svd is stateful

4. I also tried fixing the test "test_singular_value" by changing one of the algorithms to "randomized", the result is that the some singular values fail matching at the first digit following the decimal point:

    apca = TruncatedSVD(n_components=2, algorithm='arpack',
                        random_state=rng).fit(X)
    rpca = TruncatedSVD(n_components=2, algorithm='randomized',
                        random_state=rng).fit(X)
   

  assert_array_almost_equal(apca.singular_values_, rpca.singular_values_, 12)

E AssertionError:
E Arrays are not almost equal to 12 decimals
E
E (mismatch 100.0%)
E x: array([17.574469336247, 17.470009192113])
E y: array([17.522384124448, 17.359669704603])

My consequence is that "test_singular_values" needs to be fixed, however, I find changing the precision to 0 digits a too heavy intervention, although this difference in precision is what I would expect when a randomized method is involved.

Would it make sense to play around with the parameters (like the size of the sample matrix) and try to find a combination which ensures at least 3-4 digits matching?

[jnothman]

Are you sure it's the change to the test, not the change to the implementation, that resulted in the failure?

[jnothman]

I suspect that when that test comparing rpca to apca was adapted from test_pca, it was seen as inappropriate to compare apca and rpca, and somehow randomized was changed to arpack without thinking it through... not sure though.

[jnothman]

So I don't think we should be testing randomized there, at least not in this PR.

[gykovacs]

I did some further checks, let me summarize briefly:

1. originally, I added a fix to the bug, and all tests passed (https://travis-ci.org/scikit-learn/scikit-
   learn/builds/474089316?utm_source=github_status&utm_medium=notification)

2. then, I added a test checking the labels in coclustering (which are implied by the ordering of svds), and surprisingly, a different test failed (https://travis-ci.org/scikit-learn/scikit-learn/builds/482023389?utm_source=github_status&utm_medium=notification)

3. then, I put a SkipException into the test I added (to prevent from running), and all tests passed again (with my bugfix being part of the package), which suggests that not the bugfix, but the execution of the additional test causes the failure of that other test (https://travis-ci.org/scikit-learn/scikit-learn/builds/482182493?utm_source=github_status&utm_medium=notification)

4. One more observation on why scipy arpack seems to be stateful: as you already mentioned, the failing test is not correct, as rpca uses arpack just like apca. Both of these objects are initialized by the same random state and parameters, thus, the objects are identical and should give identical results. Contrarily, the test fails because of a tiny difference in the 12th digit. In my interpretation there should be no way of having any difference between the results unless some execution path maintains some state or uncontrolled randomness, and my guess is that it comes from scipy, as it is the only common point between the two tests.

Now, I have the bugfix in the code, my test is added, the required precision of the failing (arpack vs arpack) test is decreased to 11 digits, and things seem to work fine, all tests are passing. Maybe it would worth reporting the failure of the arpack vs arpack comparison in the truncated svd test (they should match to any precision) as a separate issue.

Any comments are welcome!

[jnothman]

Thanks. I've confirmed the test failing at master.

Please add an entry to the change log at doc/whats_new/v0.21.rst. Like the other entries there, please reference this pull request with :issue: and credit yourself (and other contributors if applicable) with :user:

The arpack state issue might be worth raising separately.

[gykovacs]

Done! Thank you for the guidance! I'll raise the arpack state issue and also try to take it soon.

[cmarmo]

Hi @gykovacs , if you are still interested in working on that, do you mind fixing conflicts? Thanks for your work and your patience!

[gykovacs]

Yeah, let me look into it!