ENH: sparse.linalg: Add the feature of spectral norm of sparse matrices

[zhaog6]

Reference issue

What does this implement/fix?

In this PR, the spectral norm (2-norm and -2-norm) of the matrices are added, the singular values are computed by efficient LOBPCG indirectly.

Additional information

[lobpcg]

This can be more easily done just by calling https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/_svds.py
that includes lobpcg as one of the solvers.

[zhaog6]

Dear Andrew @lobpcg
I had written the feature by calling svds before (made rebase and squeezed the commits history, so you don't see which). But the reason I changed my mind to directly call lobpcg instead is that, when I submit the code and execute some test cases in the SciPy CI/CR, a few special test cases will not pass because of the issue of the algebraic accuracy (for example, I stiill vaguely remember the following case,

A = [[-4, -3, -2],
     [-1,  0,  1],
     [2,   3,  4]];

the error between the smallest singular value solved by svds and solved by numpy.linalg.norm is ~5.3e-8, and the error of the smallest sing. val, solved by svds and the exact smallest sing. val. is also ~5.3e-8). I checked, it has nothing to do with adjusting the parameters tol and maxiter. Finally I located the problem, because the error of the smallest eigenvalue of A'A is ~2.8e-15 (the exact smallest eigenvalues of A'A and sing. val of A is zero). In this case, if we take the square root of the smallest eigenvalue of A'A, the error of the sing. val. will be O(1e-8), which is not zero and is still a certain distance O(1e-8) from numpy's calculation results or the exact value. So I used lobpcg instead of svds.

[lobpcg]

svds should pretty much do exactly the same thing (I wrote the part that adds lobpcg to it) as you do in this PR, so I am unsure why there might be a difference in the results, especially for simply computing the single largest singular value.

The case you provide above is for the smallest singular value, so I am also unsure how it relates to the spectral norm.

[zhaog6]

@lobpcg In addition to the above reasons, there is another reason. For scipy.sparse.linalg.norm, I just want to obtain the smallest and largest sing. value without introducing additional operations. So I didn't choose scipy.sparse.linalg.eigen.svds, although it's easy to call.

[zhaog6]

svds should pretty much do exactly the same thing (I wrote the part that adds lobpcg to it) as you do in this PR, so I am unsure why there might be a difference in the results, especially for simply computing the single largest singular value.

The case you provide above is for the smallest singular value, so I am also unsure how it relates to the spectral norm.

Dear Andrew,
In fact, on my server and laptop, svds is no problem for the alegbraic accuracy (it is zero exactly), But it reports the accuracy problem (it became 5.3e-8) on the SciPy Github Azure platform.

[lobpcg]

Why the smallest is needed for the norm?

[zhaog6]

@lobpcg The largest sing. value is no problem on any platform, but the smallest sing. val will have problems on some specific platforms.

[zhaog6]

Why the smallest is needed for the norm?

This is because NumPy defines a - 2 norm, SciPy need to be aligned with NumPy on some operations. But from computational mathematics textbook, -2 norm is not defined.

[lobpcg]

I recollect that svds does something special and strange with the smallest singular values, so I believe that there may be a difference that may justify a separate code for -2, indeed.

[zhaog6]

I recollect that svds does something special and strange with the smallest singular values, so I believe that there may be a difference that may justify a separate code for -2, indeed.

Dear Andrew @lobpcg
Okay. If the smallest eigenvalue of A'A is smaller 1e-14, Could it be set to zero? So it is still zero after taking square root.

[lobpcg]

Possible, look at the end of svds code.

[zhaog6]

Possible, look at the end of svds code.

@lobpcg I modified the file _svds.py and implementation of 2/-2 norm of sparse matrices, Please let me know if it is not a popular implementation way.

[lobpcg]

You need to make sure that everything works in single precision as well. So you cannot use fixed constants like 1e-14
You also need unit tests in both precisions.

[zhaog6]

You need to make sure that everything works in single precision as well. So you cannot use fixed constants like 1e-14 You also need unit tests in both precisions.

Okay, I try it.

[lobpcg]

For your absolute constant use a power of np.finfo(dtype).eps with dtype=np.float64 or np.float32 depending on the case

[zhaog6]

I resubmitted the PR, used 50 * np.finfo(...).eps in the codes (7d792a7), Could I do it this way?

[lobpcg]

The expression of the threshold should also depend on the shape and probably some norm of the input, like in the original code.

[zhaog6]

The expression of the threshold should also depend on the shape and probably some norm of the input, like in the original code.

Yes, the threshold 50 in the code also indeed depend on float32 or float64. For float64, the threshold is ~20 (there are some small fluctuations), whereas for float32, the threshold is probably 30-50 (affected by different platforms)

[lobpcg]

Not good, the expression should be the same for different dtypes. Try using a fixed power or log of eps instead of multiplying eps by different constants.

[zhaog6]

Dear Andrew, setting it to 50 is only to cut it off within the fluctuation range of machine accuracy，but the fluctuation seems doesn't have much regularity and there is no formula guide to rely on.

[zhaog6]

I don't have much experience in error fluctuation caused by different platforms and different accuracy (float32 & float64). For the fixed power or log of eps, I have no idea yet currently. What value do you think is more reasonable？Thanks.

[lobpcg]

My suggestion is to try obtaining both values ~20*np.finfo(np.float64).eps and ~50*np.finfo(np.float32).eps by a single formula f(x) with x=np.finfo(dtype).eps with dtype=np.float64 or np.float32. Try to find experimentally such a fixed function f(x) to handle both cases, if possible.

[zhaog6]

Okay, I'll try to find it, thanks.

[zhaog6]

Because I can't reproduce this problem on my laptop [macOS Big Sur-64bit] and server [Ubuntu-20.04-64bit] (the problem only occurs in Win32 environment in Azure Pipelined of SciPy CI/CR), I've turn it into a draft and continue to test it by submitting the draft. If there are any email reminders/bother, please ignore them until I'm ready to please you review again. Thanks.

[zhaog6]

Dear Andrew @lobpcg I make some tests varous platform in Docker, there is no rule for the fluctuation range of the result. However, I find a bug from NumPy when I check the call process in NumPy and SciPy (Notice: the scipy/sparse/linalg/tests/test_norm.py file compare the "SciPy norm" and "Numpy norm"). The bug occurs in the svd function in numpy/linalg/linalg.py
https://github.com/numpy/numpy/blob/ba3664766d1f8c6b2521b82edd717a2f2be974f2/numpy/linalg/linalg.py#L1659-L1661

For a np.float32 matrix, s in line 1660 is np.float64, and s in line 1661 is converted to np.float32, in other words, _umath_linalg_svd_n (gufunc) computed singular values with np.float64, then the type conversion (np.float64->np.float32) is performed. It is different from svds in SciPy (svds->lobpcg->eigh->drv). So SciPy truely implemented it with np.float32 when given matrix is np.float32, but NumPy doesn't seem to be. The above are only the clues I found when debugging, for reference only.

[lobpcg]

Please submit a bug report and the fix, since you have already found it!

[zhaog6]

Please submit a bug report and the fix, since you have already found it!

Okay, I'll submit the bug report and fix it later. Thanks.

[lobpcg]

For the smallest SVD, have a look at #11829