Skip to content

BUG: array_api: linalg.cholesky returns the conjugate of the expected upper decomposition #24451

Description

@lucascolley

Describe the issue:

I expect numpy.array_api.linalg.cholesky(a, upper=True) to return the same array as scipy.linalg.cholesky(a, lower=False), the upper Cholesky decomposition of a. Instead, it returns the conjugate.

There is a possibility that this is a bug on scipy:main, but I assume that it is here instead. I can't check with np.linalg.cholesky since that only provides lower Cholesky decompositions.

xp.linalg.cholesky: https://data-apis.org/array-api/latest/extensions/generated/array_api.linalg.cholesky.html#array_api.linalg.cholesky

scipy.linalg.cholesky: https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.cholesky.html#scipy.linalg.cholesky

Maybe a conjugation should be included on this line?:

return Array._new(L).mT

Reproduce the code example:

In [1]: import numpy as np

In [2]: import numpy.array_api as npx
<ipython-input-2-f2c2051749b9>:1: UserWarning: The numpy.array_api submodule is still experimental. See NEP 47.
  import numpy.array_api as npx

In [3]: from scipy.linalg import cholesky

In [4]: a = np.asarray([[ 10. +0.j,  13. +9.j,  15. -5.j],
   ...:        [ 13. -9.j,  33. +0.j,  33.-10.j],
   ...:        [ 15. +5.j,  33.+10.j, 143. +0.j]])

In [5]: cholesky(a, lower=False)
Out[5]: 
array([[3.16227766+0.j        , 4.11096096+2.84604989j,
        4.74341649-1.58113883j],
       [0.        +0.j        , 2.82842712+0.j        ,
        6.36396103+3.53553391j],
       [0.        +0.j        , 0.        +0.j        ,
        8.06225775+0.j        ]])

In [6]: b = npx.asarray([[ 10. +0.j,  13. +9.j,  15. -5.j],
   ...:        [ 13. -9.j,  33. +0.j,  33.-10.j],
   ...:        [ 15. +5.j,  33.+10.j, 143. +0.j]])

In [7]: npx.linalg.cholesky(b, upper=True)
Out[7]: 
Array([[3.16227766+0.j        , 4.11096096-2.84604989j,
        4.74341649+1.58113883j],
       [0.        +0.j        , 2.82842712+0.j        ,
        6.36396103-3.53553391j],
       [0.        +0.j        , 0.        +0.j        ,
        8.06225775+0.j        ]], dtype=complex128)

Error message:

No response

Runtime information:

import sys, numpy; print(numpy.__version__); print(sys.version):

1.25.0
3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:40:32) [GCC 12.3.0]

print(numpy.show_runtime()):

[{'numpy_version': '1.25.0',
  'python': '3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:40:32) '
            '[GCC 12.3.0]',
  'uname': uname_result(system='Linux', node='lucas-pc', release='6.2.0-26-generic', version='#26~22.04.1-Ubuntu SMP PREEMPT_DYNAMIC Thu Jul 13 16:27:29 UTC 2', machine='x86_64')},
 {'simd_extensions': {'baseline': ['SSE', 'SSE2', 'SSE3'],
                      'found': ['SSSE3',
                                'SSE41',
                                'POPCNT',
                                'SSE42',
                                'AVX',
                                'F16C',
                                'FMA3',
                                'AVX2'],
                      'not_found': ['AVX512F',
                                    'AVX512CD',
                                    'AVX512_SKX',
                                    'AVX512_CLX',
                                    'AVX512_CNL',
                                    'AVX512_ICL',
                                    'AVX512_SPR']}},
 {'filepath': '/home/lucas/mambaforge/envs/scipy-dev/lib/libmkl_rt.so.2',
  'internal_api': 'mkl',
  'num_threads': 6,
  'prefix': 'libmkl_rt',
  'threading_layer': 'intel',
  'user_api': 'blas',
  'version': '2022.2-Product'},
 {'filepath': '/home/lucas/mambaforge/envs/scipy-dev/lib/libomp.so',
  'internal_api': 'openmp',
  'num_threads': 12,
  'prefix': 'libomp',
  'user_api': 'openmp',
  'version': None},
 {'architecture': 'Zen',
  'filepath': '/home/lucas/mambaforge/envs/scipy-dev/lib/libopenblasp-r0.3.23.so',
  'internal_api': 'openblas',
  'num_threads': 12,
  'prefix': 'libopenblas',
  'threading_layer': 'pthreads',
  'user_api': 'blas',
  'version': '0.3.23'}]
None

Context for the issue:

I am working on array API support in scipy.linalg. This is failing tests over there for numpy.array_api arrays.

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions