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

[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?:

numpy/numpy/array_api/linalg.py

Line 56 in 2ce972e

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.

[lucascolley]

cc @asmeurer

[izaid]

Also noting the docs for the array API need to be updated. It's correct for the lower decomposition, but for the upper decomposition only it should read x = U^H U not x = U U^H.

[ntessore]

Cf. LAPACK:

 ZPOTRF computes the Cholesky factorization of a complex Hermitian
 positive definite matrix A.

 The factorization has the form
    A = U**H * U,  if UPLO = 'U', or
    A = L  * L**H,  if UPLO = 'L',
 where U is an upper triangular matrix and L is lower triangular.