BUG: complex determinants broken with openblas 0.3.22 #18208

mkoeppe · 2023-03-28T22:16:01Z

Describe your issue.

scipy.linalg.det is broken when scipy is built from source with the just-released openblas 0.3.22.

Tested on macOS x86_64.

Reproducing Code Example

>>> import numpy, scipy
>>> A = numpy.array([[1,2,4],[5,3,9],[7,8,6]], dtype='complex')
>>> A
array([[1.+0.j, 2.+0.j, 4.+0.j],
       [5.+0.j, 3.+0.j, 9.+0.j],
       [7.+0.j, 8.+0.j, 6.+0.j]])
>>> scipy.linalg.det(A)
0j      # wrong
>>> Ar = numpy.array([[1,2,4],[5,3,9],[7,8,6]], dtype='float')
>>> scipy.linalg.det(Ar)
88.0

Error message

N/A

SciPy/NumPy/Python version and system information

>>> import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info); scipy.show_config()
1.10.1 1.24.2 sys.version_info(major=3, minor=11, micro=2, releaselevel='final', serial=0)
/Users/mkoeppe/s/sage/sage-rebasing/worktree-algebraic-2018-spring/local/var/lib/sage/venv-python3.11/lib/python3.11/site-packages/scipy/__config__.py:140: UserWarning: Install `pyyaml` for better output
  warnings.warn("Install `pyyaml` for better output", stacklevel=1)
{
  "Compilers": {
    "c": {
      "name": "clang",
      "linker": "ld64",
      "version": "14.0.0",
      "commands": "gcc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "0.29.32",
      "commands": "cython"
    },
    "c++": {
      "name": "clang",
      "linker": "ld64",
      "version": "14.0.0",
      "commands": "g++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld64",
      "version": "12.2.0",
      "commands": "gfortran"
    },
    "pythran": {
      "version": "0.12.1",
      "include directory": "/Users/mkoeppe/s/sage/sage-rebasing/worktree-algebraic-2018-spring/local/var/lib/sage/venv-python3.11/lib/python3.11/site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "darwin"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "darwin"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "openblas",
      "found": true,
      "version": "0.3.22",
      "detection method": "pkgconfig",
      "include directory": "/usr/local/Cellar/openblas/0.3.22/include",
      "lib directory": "/usr/local/Cellar/openblas/0.3.22/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP=1 NEHALEM MAX_THREADS=56",
      "pc file directory": "/usr/local/opt/openblas/lib/pkgconfig"
    },
    "lapack": {
      "name": "openblas",
      "found": true,
      "version": "0.3.22",
      "detection method": "pkgconfig",
      "include directory": "/usr/local/Cellar/openblas/0.3.22/include",
      "lib directory": "/usr/local/Cellar/openblas/0.3.22/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP=1 NEHALEM MAX_THREADS=56",
      "pc file directory": "/usr/local/opt/openblas/lib/pkgconfig"
    }
  },
  "Python Information": {
    "path": "/Users/mkoeppe/s/sage/sage-rebasing/worktree-algebraic-2018-spring/local/var/lib/sage/venv-python3.11/bin/python3",
    "version": "3.11"
  }
}

The text was updated successfully, but these errors were encountered:

ilayn · 2023-03-28T22:30:01Z

Are the NumPy versions the same? Asking to eliminate f2py sneaking in with the Fortran stuff.

ilayn · 2023-03-28T22:35:29Z

I can only see these changes numpy/numpy#22025 & OpenMathLib/OpenBLAS#3924 that might be relevant on OpenBLAS regarding ?getrf routines.

By the way, why are we even using fortran files for these things ? After so many years still finding out new quirks in linalg folder 😃

mkoeppe · 2023-03-28T22:37:43Z

I've tested with several numpy versions, including the published wheels; does not seem to make a difference.

isuruf · 2023-03-28T22:38:33Z

See also OpenMathLib/OpenBLAS#3976

andyfaff · 2023-03-28T22:55:39Z

I would hope that examples like the OP presented would be present in a test suite. We should definitely add something like this to the scipy test suite.

mkoeppe · 2023-03-28T23:03:08Z

Yes, we caught it using SageMath's test suite.

dimpase · 2023-03-28T23:18:43Z

LAPACK does not have a function to compute determinants. This means that one of LAPACK factorisation routines, e.g. LU decomposition, or whatever scipy is using here, got buggy.

Can scipy people tell us how they do it, what are they calling ? Perhaps experimenting with calling various matrix decompositions will tell more.

andyfaff · 2023-03-29T03:52:59Z

The code says:

The determinant is computed via LU factorization, LAPACK routine z/dgetrf.

The callstack is something like:

det --> zdet_r --> zdet_c --> call to zgetrf

ilayn · 2023-03-29T08:31:06Z

I'll fix that today or tomorrow evening and get rid of the custom fortran files and call ?getrf directly.. They are not doing anything substantial. Complex tests are also wrong it seems so I can add those fixes too.

ilayn · 2023-03-29T16:29:28Z

apologies for hijacking the issue but could anyone reading (especially for conda users with MKL.), can run the following snippet and post the plot here? I kept it comparable to the existing implementation with the excessive checks and whatnot;

# pip install perfplot
import perfplot
import numpy as np
import scipy.linalg as la
from scipy.linalg._decomp import _asarray_validated
from scipy.linalg import get_lapack_funcs
from scipy.linalg._misc import _datacopied


def mydet(a, overwrite_a=False):
    a1 = _asarray_validated(a, check_finite=True)
    if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
        raise ValueError('expected square matrix')
    overwrite_a = overwrite_a or _datacopied(a1, a)
    trf = get_lapack_funcs('getrf', dtype=a1.dtype)
    # det(a) = det(a.T) but can overwrite if dtype is compatible to (s,d,c,z)
    if a1.flags['C_CONTIGUOUS']:
        lu, piv, info = trf(a1.T, overwrite_a=overwrite_a)
    else:
        lu, piv, info = trf(a1, overwrite_a=overwrite_a)

    if info < 0:
        raise ValueError(f'illegal value in {-info}-th argument of internal '
                         'det.getrf')
    a_det = 1.
    # piv is 0-indexed
    for x in range(a1.shape[0]):
        a_det *= lu[x, x]*(1. if piv[x] == x else -1.)
    return a_det


perfplot.show(
    setup=lambda n: np.random.rand(n, n),
    # setup=lambda n: np.random.rand(n, n) + np.random.rand(n, n)*-2j,
    kernels=[la.det, mydet],
    labels=["SciPy", "this PR"],
    n_range=[3*x for x in range(1, 100, 3)],
    logx="auto",
    logy="True",
    equality_check=np.allclose,
    xlabel="input size (n)",
)

My machine is again smoking inappropriate stuff whenever partying with LAPACK, pure python is almost the same performance (which is understandable since most time is spent on getrf) but just wanted to make sure .

andyfaff · 2023-03-29T21:16:34Z

Looks like the differences are larger. Real on top, complex on bottom
macosx_arm64, scipy 1.10.1, python3.9

j-bowhay · 2023-03-29T22:18:50Z

real:

complex:

In [4]: import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info); scipy.show_config()
1.11.0.dev0+1761.02fb953 1.23.5 sys.version_info(major=3, minor=10, micro=8, releaselevel='final', serial=0)
/home/jakeb/development/scipy/build-install/lib/python3.10/site-packages/scipy/__config__.py:140: UserWarning: Install `pyyaml` for better output
  warnings.warn("Install `pyyaml` for better output", stacklevel=1)
{
  "Compilers": {
    "c": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "11.3.0",
      "commands": "/home/jakeb/miniconda3/envs/scipy-dev/bin/x86_64-conda-linux-gnu-cc"
    },
    "cython": {
      "name": "cython",
      "linker": "cython",
      "version": "0.29.33",
      "commands": "cython"
    },
    "c++": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "11.3.0",
      "commands": "/home/jakeb/miniconda3/envs/scipy-dev/bin/x86_64-conda-linux-gnu-c++"
    },
    "fortran": {
      "name": "gcc",
      "linker": "ld.bfd",
      "version": "11.3.0",
      "commands": "/home/jakeb/miniconda3/envs/scipy-dev/bin/x86_64-conda-linux-gnu-gfortran"
    },
    "pythran": {
      "version": "0.12.1",
      "include directory": "../../../miniconda3/envs/scipy-dev/lib/python3.10/site-packages/pythran"
    }
  },
  "Machine Information": {
    "host": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "linux"
    },
    "build": {
      "cpu": "x86_64",
      "family": "x86_64",
      "endian": "little",
      "system": "linux"
    },
    "cross-compiled": false
  },
  "Build Dependencies": {
    "blas": {
      "name": "openblas",
      "found": true,
      "version": "0.3.21",
      "detection method": "pkgconfig",
      "include directory": "/home/jakeb/miniconda3/envs/scipy-dev/include",
      "lib directory": "/home/jakeb/miniconda3/envs/scipy-dev/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP=0 PRESCOTT MAX_THREADS=2",
      "pc file directory": "/home/jakeb/miniconda3/envs/scipy-dev/lib/pkgconfig"
    },
    "lapack": {
      "name": "openblas",
      "found": true,
      "version": "0.3.21",
      "detection method": "pkgconfig",
      "include directory": "/home/jakeb/miniconda3/envs/scipy-dev/include",
      "lib directory": "/home/jakeb/miniconda3/envs/scipy-dev/lib",
      "openblas configuration": "USE_64BITINT= DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK= NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP=0 PRESCOTT MAX_THREADS=2",
      "pc file directory": "/home/jakeb/miniconda3/envs/scipy-dev/lib/pkgconfig"
    }
  },
  "Python Information": {
    "path": "/home/jakeb/miniconda3/envs/scipy-dev/bin/python3.10",
    "version": "3.10"
  }
}

dimpase · 2023-03-29T22:36:34Z

these benchmarks are obviously very sensitive w.r.t. number of threads, and MAX_THREADS=2 is serioisly low. How about at least 8? And it's not at all clear what's used for this PR.

ilayn · 2023-03-29T22:54:12Z

multithreading does (should) not start until some value later in the plot. That kinda common jump is often L2 cache limit. But one can never be sure. "This PR" is the one I'm about to prepare which is the mydet function. It's just a label, see the code above.

ilayn · 2023-04-01T11:19:55Z

All feedback is welcome in the linked PR. I tried to address the shortcomings I could gather from various comments and issues.

ilayn · 2023-04-09T19:26:59Z

linalg.det function is also worked out and more tests are added accordingly. Since we identified this problem as an upstream issue, we can close it. But thanks everyone for the heads up and contributions.

mkoeppe added the defect A clear bug or issue that prevents SciPy from being installed or used as expected label Mar 28, 2023

ilayn added the scipy.linalg label Mar 28, 2023

mkoeppe mentioned this issue Mar 28, 2023

build/pkgs/openblas/spkg-configure.m4: Reject version 0.3.22 sagemath/sage#35377

Merged

5 tasks

mmuetzel mentioned this issue Mar 29, 2023

openblas: update to 0.3.23 msys2/MINGW-packages#16467

Merged

dimpase mentioned this issue Mar 29, 2023

0.3.22 breaks functions downstream in GNU Octave OpenMathLib/OpenBLAS#3976

Closed

ilayn mentioned this issue Apr 1, 2023

ENH:MAINT:linalg det in Cython and with nDarray support #18225

Merged

ilayn closed this as completed Apr 9, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BUG: complex determinants broken with openblas 0.3.22 #18208

BUG: complex determinants broken with openblas 0.3.22 #18208

mkoeppe commented Mar 28, 2023

ilayn commented Mar 28, 2023

ilayn commented Mar 28, 2023

mkoeppe commented Mar 28, 2023

isuruf commented Mar 28, 2023

andyfaff commented Mar 28, 2023

mkoeppe commented Mar 28, 2023

dimpase commented Mar 28, 2023 •

edited

Loading

andyfaff commented Mar 29, 2023

ilayn commented Mar 29, 2023 •

edited

Loading

ilayn commented Mar 29, 2023 •

edited

Loading

andyfaff commented Mar 29, 2023 •

edited

Loading

j-bowhay commented Mar 29, 2023

dimpase commented Mar 29, 2023

ilayn commented Mar 29, 2023

ilayn commented Apr 1, 2023

ilayn commented Apr 9, 2023

BUG: complex determinants broken with openblas 0.3.22 #18208

BUG: complex determinants broken with openblas 0.3.22 #18208

Comments

mkoeppe commented Mar 28, 2023

Describe your issue.

Reproducing Code Example

Error message

SciPy/NumPy/Python version and system information

ilayn commented Mar 28, 2023

ilayn commented Mar 28, 2023

mkoeppe commented Mar 28, 2023

isuruf commented Mar 28, 2023

andyfaff commented Mar 28, 2023

mkoeppe commented Mar 28, 2023

dimpase commented Mar 28, 2023 • edited Loading

andyfaff commented Mar 29, 2023

ilayn commented Mar 29, 2023 • edited Loading

ilayn commented Mar 29, 2023 • edited Loading

andyfaff commented Mar 29, 2023 • edited Loading

j-bowhay commented Mar 29, 2023

dimpase commented Mar 29, 2023

ilayn commented Mar 29, 2023

ilayn commented Apr 1, 2023

ilayn commented Apr 9, 2023

dimpase commented Mar 28, 2023 •

edited

Loading

ilayn commented Mar 29, 2023 •

edited

Loading

ilayn commented Mar 29, 2023 •

edited

Loading

andyfaff commented Mar 29, 2023 •

edited

Loading