Testing test_components_from_metric_edge_cases breaks around line 130 on some environments

[RobinVogel]

This might be a particularly weird and not primordial issue.

Description

One of the tests fails on a newly created environment, depending on the platform I test on.

Steps/Code to Reproduce

With Python 3.7.5 installed with Anaconda:

conda create -n test_env numpy scipy pytest scikit-learn
conda activate test_env
pytest test/test_components_metric_conversion.py

Expected Results

Tests run without failure on any platform.

Actual Results

Fails around line 130 of test/test_components_metric_conversion.py for on the two platform tested.

============================================== FAILURES 
[...]
      assert np.linalg.slogdet(M)[1] > 1  # (just to show that the computed
      # determinant is far from null)
      with pytest.raises(LinAlgError) as err_msg:
>       np.linalg.cholesky(M)
E       Failed: DID NOT RAISE <class 'numpy.linalg.LinAlgError'>

test/test_components_metric_conversion.py:130: Failed

Versions

The environment in which it fails:

Linux-4.18.0-25-generic-x86_64-with-debian-buster-sid
Python 3.7.4 (default, Aug 13 2019, 20:35:49) 
[GCC 7.3.0]
NumPy 1.17.2
SciPy 1.3.1
Scikit-Learn 0.21.3
Metric-Learn 0.5.0

The environment in which it passes:

Linux-5.0.0-32-generic-x86_64-with-debian-buster-sid
Python 3.7.4 (default, Aug 13 2019, 20:35:49) 
[GCC 7.3.0]
NumPy 1.17.2
SciPy 1.3.1
Scikit-Learn 0.21.3
Metric-Learn 0.5.0

So the only difference is the platform.

[RobinVogel]

After some investigation it seems that the matrices generated by ortho_group.rvs are slightly different between the two computers. I tested for equality after saving those with np.save and generating with ortho_group.rvs(6, random_stage=rng) on the two machines. It is illustrated below.

One may chose to modify the test to make it robust to numerical instability, or understand the discrepancies in the generation of the matrices by scipy.

Illustration: both matrices were generated with:

import numpy as np
from scipy.stats import ortho_group

rng = np.random.RandomState(42)
P = ortho_group.rvs(6, random_state=rng)

and saved with np.save, loaded with np.load. They are different, as seen below:

In [1]: import numpy as np                                                                             

In [2]: P1 = np.load("P_laptop.npy"); P2 = np.load("P_desktop.npy")                                    

In [3]: P1.shape, P2.shape                                                                             
Out[3]: ((6, 6), (6, 6))

In [4]: for p1, p2 in zip(P1[P1 != P2], P2[P1 != P2]): print(p1); print(p2); print("")                 
-0.18193963255076545
-0.18193963255076542

0.2701832074139061
0.27018320741390633

-0.6994139929371295
-0.6994139929371296

0.4774097779144954
0.4774097779144953

0.1515489690153527
0.15154896901535272

-0.048868922516862756
-0.048868922516862784

-0.31821384787190377
-0.31821384787190365

0.32640756874271126
0.3264075687427112

-0.6554276151588394
-0.6554276151588395

0.45848551023953366
0.4584855102395337

[bellet]

Thanks for digging into this. I think we should definitely make the test more robust. I do not immediately recall the rationale behind this part of the test - @wdevazelhes maybe you have an idea to make this test more robust or can explain quickly what was intended?

[wdevazelhes]

Thanks for digging into this. I think we should definitely make the test more robust. I do not immediately recall the rationale behind this part of the test - @wdevazelhes maybe you have an idea to make this test more robust or can explain quickly what was intended?

This part of the test was to ensure that the convention for definiteness is coherent with the one in numpy (i.e. an eigenvalue is considered null if its abs value is lower than the maximum abs eigenvalue times the lenght of the matrix times a machine epsilon, see here:

metric-learn/metric_learn/_util.py

Line 352 in a7ede57

tol : positive `float`, optional

, or see matrix_rank in numpy). The example in the test ensures that for instance one does not check definiteness by computing the determinant, because sometimes the determinant is not null (like in the tested example) (which would mean definiteness in theory), but according to the numpy convention the matrix is not definite.

But it's not really clear to me yet why it fails on one machine and not the other, even if the orthogonal matrices are different... It could mean that maybe thenp.finfo(float).eps are also different in the two machines, but I don't see why they would... Otherwise a possible thing: in the test we assert undefiniteness using np.linalg.cholesky (looking at the error thrown). I guess I assumed it would be coherent with the convention in matrix_rank, but maybe it is not: i.e. maybe np.linalg.cholesky has its own way to compute definiteness that can throw no error when the matrix_rank convention should (it clearly should here because on my machine np.finfo(float).eps = 2.2204460492503131e-16, and clearly 1e-20 < 2.2204460492503131e-16 * 7 * 1e5 = 1.5543122344752192e-10)

[bellet]

After a quick look at np.linalg.cholesky it is not super clear how semi-definiteness is checked there..

[bellet]

It would be nice to fix this because it now makes the CI check fail. I can't reproduce the error on my local machine though. @wdevazelhes do you see a simple way to make this test more robust?

[wdevazelhes]

I'm not sure about it, but maybe we could replace the

    with pytest.raises(LinAlgError) as err_msg:
      np.linalg.cholesky(M)

by

assert np.linalg.rank(M) < M.shape[0]

Indeed, now that I re-read the test I don't really understand it :p In the sense that no matter if the matrix is definite or semi-definite according to any convention (even poor ones like the determinant) we have a way to deal with it and the result in these cases are already tested in other tests from test_components_from_metric_edge_cases (like the one below in this comment that checks that components_from_metric works for semidefinite matrices), so I don't really see why this test is here.
I guess we could still keep in "bonus", just to test that our cholesky decomposition works fine in some possible weird cases, like the one with a high determinant but undefinite according to some reasonable convention, like np.linalg.matrix_rank, hence the proposed modification
Now an error thrown by np.linalg.cholesky could also be a valid convention for undefiniteness, so to test really everything we could try to find a case where np.linalg.cholesky returns an Error AND the matrix has a high determinant, but I don't really know how to search for such a case, because as you said it's not clear how np.linalg.cholesky checks for definiteness, and moreover as said this test seems "bonus" i.e. less essential than say the following test (that checks that components_from_metric works fine for TRULY undefinite matrices i.e. with real zeros):

metric-learn/test/test_components_metric_conversion.py

Lines 107 to 112 in 86a5208

	# low-rank matrix (with zeros)
	M = np.zeros((7, 7))
	small_random = rng.randn(3, 3)
	M[:3, :3] = small_random.T.dot(small_random)
	L = components_from_metric(M)
	assert_allclose(L.T.dot(L), M)

We could also remove the test because I'm not sure it's a good idea to keep "bonus" tests, probably a test should check precisely some functionality... I don't know

However, this made me think that we should probably add another test__check_sdp_from_eigen_numy_convention in test_utils that here would be much more explicit and useful: it would test that the way we check for SDP is good (e.g. not the determinant), and follows the convention from numpy.linalg.rank.

I'll open a PR for that so it's easier to discuss there