Inconsistent behavior of linalg.slogdet()

[GBR-613]

There are two almost equal arrays generated by the same pipeline in different systems. (All differences between respective items are about 1e-16 and apparently caused by differences in hardware.)
I would expect result of linalg.slogdet() for them very similar, respectively.
However the result for one is about -49.04, and the result for another is -inf.
The npy files are attached.

Reproducing code example:

import numpy as np


a  = np.load("a.npy", allow_pickle=True)
s  = np.load("s.npy", allow_pickle=True)
_, det = np.linalg.slogdet(a)
print(det)
_, det = np.linalg.slogdet(s)
print(det)

# The section below illustrates how difference
# between these two arrays is small:
# nothing is printed
for i in range(a.shape[0]):
    for j in range(a.shape[1]):
        d = abs(a[i][j] - s[i][j])
        if d > 0.00000000000001:
            print("%d, %d" % (i, j))
            print(a[i][i])
            print(s[i][i])
            print(d)

Error message:

The code below prints the following:
-inf
-49.041742650541245
Expected result is about -49.04 in both cases.
(FYI: this inconsistency leads to 20-30% difference in final results in our project on different systems.)

NumPy/Python version information:

1.20.2 3.9.4 (default, May 4 2021, 16:03:03)
[GCC 9.3.0]

The same behaviour is with Python 3.7.10 on both Windows and Linux.

[GBR-613]

npy.zip
Since Github does not allow npy attachments, they are added to the attached zip.

[charris]

The slogdet is very small, so my suspicion is that the matrices are almost singular, hence subject to roundoff error. What are the singular values?

[GBR-613]

@charris not sure I got your point. If determinant of one of the matrices is -49 and these matrices are similar, then they are hardly singular.
I would say they are almost symmetric. But why it may be a subject to roundoff error?
You can easily convert npy file into human-readable format:

import numpy as np
a = np.load("a.pny",  allow_pickle=True)
f = open("a.txt", "w")
f.write(str(a))
f.close()

[bashtage]

Which BLAS are you using? On MKL I see

-inf
-inf

[bashtage]

The smallest eigenvalue is numerically not different from 0. Small differences can make a non-trivial difference to this noise eigenvalue.

np.linalg.eigvals(a)
Out[15]: 
array([ 4.87687978e+00,  2.18850804e+00,  1.84140287e+00,  1.55259050e+00,
        1.26482134e+00,  9.19017650e-01,  8.71262541e-01,  8.47400689e-01,
        6.49832078e-01,  2.98709429e-01,  2.81602249e-01,  1.91985822e-01,
        1.29640613e-01,  1.18498928e-02,  2.14567464e-02, -7.43559363e-17])

np.linalg.eigvals(s)
Out[16]: 
array([ 4.87687978e+00,  2.18850804e+00,  1.84140287e+00,  1.55259050e+00,
        1.26482134e+00,  9.19017650e-01,  8.71262541e-01,  8.47400689e-01,
        6.49832078e-01,  2.98709429e-01,  2.81602249e-01,  1.91985822e-01,
        1.29640613e-01,  1.18498928e-02,  2.14567464e-02, -2.67468955e-17])

Both np.linalg.matrix_rank and the LAPACK function pstrf, which can be called using SciPy, agree that these are rank 15 arrays, and so -inf is the right answer.

[GBR-613]

@bashtage numpy.show_config() is not informative now. I only can tell you that one of the systems is Windows 10 with processor Intel(R) Core(TM) i7-9850H, no GPU, and I used default 64-bit installer from python.org. But I can not guarantee that it is MKL.

Results of np.linalg.eigvals() in my system:

[4.87687978e+00 2.18850804e+00 1.84140287e+00 1.55259050e+00
1.26482134e+00 9.19017650e-01 8.71262541e-01 8.47400689e-01
6.49832078e-01 2.98709429e-01 2.81602249e-01 1.91985822e-01
1.29640613e-01 1.18498928e-02 2.14567464e-02 1.72487350e-17]

[4.87687978e+00 2.18850804e+00 1.84140287e+00 1.55259050e+00
1.26482134e+00 9.19017650e-01 8.71262541e-01 8.47400689e-01
6.49832078e-01 2.98709429e-01 2.81602249e-01 1.91985822e-01
1.29640613e-01 1.18498928e-02 2.14567464e-02 4.58920394e-17]

[GBR-613]

@bashtage OK, then there is an amendment to the issue:
Expected result is -inf in both cases.

[GBR-613]

@bashtage @charris Interesting update:
I have to confess that my Linux box is a VM running on AWS with Intel 64-bit processor. (BTW there is an option for free tier, so you may use it to reproduce the problem.)
AWS also has an option to create VM with ARM processor. I did, compiled last Python3.9 from sources, installed last numpy (which is 1.20.3) and ran the test script. Result is:
-49.041742650541245
-48.3485954699813
This is actually the result I would expect. (Although np.linalg.matrix_rank(a) returns 15 in ARM processor, too!)

Next, I tried to run np.linalg.eigvals(a)
Result:
numpy.linalg.LinAlgError: Eigenvalues did not converge
The same error is with np.linalg.eigh(a)
Issue #16744 recommends upgrading numpy to 1.19.5 to solve the problem, but now I am in the bigger version!
What do you think?

[bashtage]

You are using OpenBLAS. That bug only afffects Windows so isn't your problem in AWS. I think this is just an edge case. If you think your covariance may not be positive definite then you should check and use an appropriate correction. I think this is just squarely in the area of numerical precision limits which is why you are seeing these mixed results when you try different platforms or even different algorithms.

For example:

np.prod(np.linalg.eigvals(a)), np.log(np.abs(np.prod(np.linalg.eigvals(a))))

 (-6.735614088448808e-22, -48.74946306147741)

and since a and s are symmetric,

np.prod(np.linalg.eigvalsh(a)), np.log(np.abs(np.prod(np.linalg.eigvalsh(a))))

(4.085687572498557e-22, -49.24938201531954)

but

np.linalg.inv(a)

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-56-b24b020fe220> in <module>
----> 1 np.linalg.inv(a)

<__array_function__ internals> in inv(*args, **kwargs)

C:\Anaconda\lib\site-packages\numpy\linalg\linalg.py in inv(a)
    543     signature = 'D->D' if isComplexType(t) else 'd->d'
    544     extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 545     ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    546     return wrap(ainv.astype(result_t, copy=False))
    547 

C:\Anaconda\lib\site-packages\numpy\linalg\linalg.py in _raise_linalgerror_singular(err, flag)
     86 
     87 def _raise_linalgerror_singular(err, flag):
---> 88     raise LinAlgError("Singular matrix")
     89 
     90 def _raise_linalgerror_nonposdef(err, flag):

LinAlgError: Singular matrix