Repository: scipy/scipy
Pull request: BUG: linalg: fix sqrtm precision preservation for complex input.
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Fix #18250 

#### What does this implement/fix?
<!--Please explain your changes.-->
I fixed the `sqrtm` precision preservation for complex input.
As discussed in #18250, when input is complex, output precision should be the same as the input data size, not twice.

#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
BUG: `complex256` type created by entering `complex128` into `linalg.sqrtm` causes errors.
### Describe your issue.

When `scipy.linalg` operates on a complex array/matrix it returns a complex array with double the precision of the input array. As a result if given a `complex128` type array it will return a `complex256` array.

This type of array is not supported by `numpy.linalg` functions and will result in errors. And  `scipi.linalg` functions it will automatically downcast any `complex256` type to a `complex128` type.

This causes many problems in large code bases with combine `scipy` and `numpy`. 

I have attached a sample code which shows the behavior of both scipy and numpy when acting on this datatype.

### Reproducing Code Example

```python
import numpy as np
import scipy as sp

# create a numpy array of type complex128
a = np.array([[1.0, 3.0], [1.0, 4.0]],dtype=np.complex128)
print(type(np.square(a[0][0])))

# sqrtm returns an array of complex256
b = sp.linalg.sqrtm(a)
print(type(np.square(b[0][0])))

# arrays of type complex128 will run fine for both
# scipy and numpy functions
print(type(np.linalg.eig(a)[0][0]))
print(type(sp.linalg.eig(a)[0][0]))

# arrays of type complex256 will run for
# scipy functions, but not for the exact same numpy function
# it would seem that scipy functions reduce the precision of
# the complex256 inputs while numpy functions throw a TypeError
print(type(sp.linalg.eig(b)[0][0]))
print(type(np.linalg.eig(b)[0][0]))
```


### Error message

```shell
<class 'numpy.complex128'>
<class 'numpy.complex256'>
<class 'numpy.complex128'>
<class 'numpy.complex128'>
<class 'numpy.ndarray'>
Traceback (most recent call last):
  File "/home/cphys/projects/scipy_numpy_bug/scipy_test.py", line 20, in <module>
    print(type(np.linalg.qr(b)[0][0]))
  File "<__array_function__ internals>", line 200, in qr
  File "/home/cphys/.local/lib/python3.10/site-packages/numpy/linalg/linalg.py", line 915, in qr
    t, result_t = _commonType(a)
  File "/home/cphys/.local/lib/python3.10/site-packages/numpy/linalg/linalg.py", line 147, in _commonType
    raise TypeError("array type %s is unsupported in linalg" %
TypeError: array type complex256 is unsupported in linalg
```


### SciPy/NumPy/Python version and system information

```shell
1.10.0 1.24.2 sys.version_info(major=3, minor=10, micro=6, releaselevel='final', serial=0)
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_KNL,AVX512_KNM,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX1.10.0 1.24.2 sys.version_info(major=3, minor=10, micro=6, releaselevel='final', serial=0)
blas_mkl_info:
  NOT AVAILABLE
blis_info:
  NOT AVAILABLE
openblas_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
blas_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
lapack_mkl_info:
  NOT AVAILABLE
openblas_lapack_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
lapack_opt_info:
    libraries = ['openblas', 'openblas']
    library_dirs = ['/usr/local/lib']
    language = c
    define_macros = [('HAVE_CBLAS', None)]
    runtime_library_dirs = ['/usr/local/lib']
Supported SIMD extensions in this NumPy install:
    baseline = SSE,SSE2,SSE3
    found = SSSE3,SSE41,POPCNT,SSE42,AVX,F16C,FMA3,AVX2
    not found = AVX512F,AVX512CD,AVX512_KNL,AVX512_KNM,AVX512_SKX,AVX512_CLX,AVX512_CNL,AVX512_ICL
512_ICL
```

Hints from issue discussion:
That's probably an oversight that came from a contributor with a Windows machine (maybe even from myself). Because on Win, complex256 is complex128 anyways. We don't use complex256 anywhere. I'll check. 
For reference, this was added intentionally in gh-16294 in response to gh-14853, but as noted in the last comment there, this was probably not the correct thing to do. @AtsushiSakai, can you explain why the precision was doubled in that case? I concur with @WarrenWeckesser's comment. It doesn't seem like the right thing to do to me.
I apologize if my PR has caused any issues.

I thought if the input dtype of the `sqrtm` function is real and the output dtype is complex, the precision had to be doubled to maintain accuracy (e.g., if the input was int32, the output would be complex64). 

However, it is true that there may not have been a need to double the precision if both the input and output were complex. In the current implementation, the decision is made based only on the output dtype of the calculation result, so using the input dtype as well might allow for maintaining precision when both input and output are complex.

 I will try to implement it that way and see if it works.
Yes, when the input is real floating, the complex output will be twice the `itemsize` of the input, but that should not be the case when the input is complex. The real part of the input needs to be the same as the real part of the output.

One more thing: in the case of integer input, you need to check the size of the input _after_ it gets converted to a floating point type by `_asarray_validated(A, check_finite=True, as_inexact=True)`. That's the relevant dtype to be looking at, not the integer's `itemsize`.

Commit message: BUG: linalg: fix sqrtm precision preservation based on review comment.

Patch to review:
diff --git a/scipy/linalg/_matfuncs_sqrtm.py b/scipy/linalg/_matfuncs_sqrtm.py
index 867d67cc263f..1e59fae5c6df 100644
--- a/scipy/linalg/_matfuncs_sqrtm.py
+++ b/scipy/linalg/_matfuncs_sqrtm.py
@@ -134,9 +134,7 @@ def sqrtm(A, disp=True, blocksize=64):
     sqrtm : (N, N) ndarray
         Value of the sqrt function at `A`. The dtype is float or complex.
         The precision (data size) is determined based on the precision of
-        input `A`. When the dtype is float, the precision is the same as `A`.
-        When the dtype is complex, the precision is double that of `A`. The
-        precision might be clipped by each dtype precision range.
+        input `A`.
 
     errest : float
         (if disp == False)
@@ -163,7 +161,6 @@ def sqrtm(A, disp=True, blocksize=64):
            [ 1.,  4.]])
 
     """
-    byte_size = np.asarray(A).dtype.itemsize
     A = _asarray_validated(A, check_finite=True, as_inexact=True)
     if len(A.shape) != 2:
         raise ValueError("Non-matrix input to matrix function.")
@@ -181,16 +178,12 @@ def sqrtm(A, disp=True, blocksize=64):
         R = _sqrtm_triu(T, blocksize=blocksize)
         ZH = np.conjugate(Z).T
         X = Z.dot(R).dot(ZH)
-        if not np.iscomplexobj(X):
-            # float byte size range: f2 ~ f16
-            X = X.astype(f"f{np.clip(byte_size, 2, 16)}", copy=False)
-        else:
-            # complex byte size range: c8 ~ c32.
-            # c32(complex256) might not be supported in some environments.
-            if hasattr(np, 'complex256'):
-                X = X.astype(f"c{np.clip(byte_size*2, 8, 32)}", copy=False)
-            else:
-                X = X.astype(f"c{np.clip(byte_size*2, 8, 16)}", copy=False)
+
+        if A.dtype == np.float32 and not np.iscomplexobj(X):
+            X = X.astype(np.float32, copy=False)
+        elif (A.dtype == np.float32 and np.iscomplexobj(X)) or \
+                (A.dtype == np.complex64):
+            X = X.astype(np.complex64, copy=False)
     except SqrtmError:
         failflag = True
         X = np.empty_like(A)
diff --git a/scipy/linalg/tests/test_matfuncs.py b/scipy/linalg/tests/test_matfuncs.py
index 70c55a957b16..ac2e953209da 100644
--- a/scipy/linalg/tests/test_matfuncs.py
+++ b/scipy/linalg/tests/test_matfuncs.py
@@ -412,69 +412,64 @@ def test_gh7839(self):
 
     def test_data_size_preservation_uint_in_float_out(self):
         M = np.zeros((10, 10), dtype=np.uint8)
-        # input bit size is 8, but minimum float bit size is 16
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.uint16)
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.uint32)
-        assert sqrtm(M).dtype == np.float32
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.uint64)
         assert sqrtm(M).dtype == np.float64
 
     def test_data_size_preservation_int_in_float_out(self):
         M = np.zeros((10, 10), dtype=np.int8)
-        # input bit size is 8, but minimum float bit size is 16
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.int16)
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.int32)
-        assert sqrtm(M).dtype == np.float32
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.int64)
         assert sqrtm(M).dtype == np.float64
 
     def test_data_size_preservation_int_in_comp_out(self):
         M = np.array([[2, 4], [0, -2]], dtype=np.int8)
-        # input bit size is 8, but minimum complex bit size is 64
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.int16)
-        # input bit size is 16, but minimum complex bit size is 64
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.int32)
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.int64)
         assert sqrtm(M).dtype == np.complex128
 
     def test_data_size_preservation_float_in_float_out(self):
         M = np.zeros((10, 10), dtype=np.float16)
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.float32)
         assert sqrtm(M).dtype == np.float32
         M = np.zeros((10, 10), dtype=np.float64)
         assert sqrtm(M).dtype == np.float64
         if hasattr(np, 'float128'):
             M = np.zeros((10, 10), dtype=np.float128)
-            assert sqrtm(M).dtype == np.float128
+            assert sqrtm(M).dtype == np.float64
 
     def test_data_size_preservation_float_in_comp_out(self):
         M = np.array([[2, 4], [0, -2]], dtype=np.float16)
-        # input bit size is 16, but minimum complex bit size is 64
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.float32)
         assert sqrtm(M).dtype == np.complex64
         M = np.array([[2, 4], [0, -2]], dtype=np.float64)
         assert sqrtm(M).dtype == np.complex128
         if hasattr(np, 'float128') and hasattr(np, 'complex256'):
             M = np.array([[2, 4], [0, -2]], dtype=np.float128)
-            assert sqrtm(M).dtype == np.complex256
+            assert sqrtm(M).dtype == np.complex128
 
     def test_data_size_preservation_comp_in_comp_out(self):
         M = np.array([[2j, 4], [0, -2j]], dtype=np.complex64)
+        assert sqrtm(M).dtype == np.complex64
+        M = np.array([[2j, 4], [0, -2j]], dtype=np.complex128)
         assert sqrtm(M).dtype == np.complex128
         if hasattr(np, 'complex256'):
-            M = np.array([[2j, 4], [0, -2j]], dtype=np.complex128)
-            assert sqrtm(M).dtype == np.complex256
             M = np.array([[2j, 4], [0, -2j]], dtype=np.complex256)
-            assert sqrtm(M).dtype == np.complex256
+            assert sqrtm(M).dtype == np.complex128
 
 
 class TestFractionalMatrixPower:


Code review comments:
diff_hunk: @@ -412,69 +412,64 @@ def test_gh7839(self):
 
     def test_data_size_preservation_uint_in_float_out(self):
         M = np.zeros((10, 10), dtype=np.uint8)
-        # input bit size is 8, but minimum float bit size is 16
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.uint16)
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.uint32)
-        assert sqrtm(M).dtype == np.float32
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.uint64)
         assert sqrtm(M).dtype == np.float64
 
     def test_data_size_preservation_int_in_float_out(self):
         M = np.zeros((10, 10), dtype=np.int8)
-        # input bit size is 8, but minimum float bit size is 16
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.int16)
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.int32)
-        assert sqrtm(M).dtype == np.float32
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.int64)
         assert sqrtm(M).dtype == np.float64
 
     def test_data_size_preservation_int_in_comp_out(self):
         M = np.array([[2, 4], [0, -2]], dtype=np.int8)
-        # input bit size is 8, but minimum complex bit size is 64
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.int16)
-        # input bit size is 16, but minimum complex bit size is 64
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.int32)
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.int64)
         assert sqrtm(M).dtype == np.complex128
 
     def test_data_size_preservation_float_in_float_out(self):
         M = np.zeros((10, 10), dtype=np.float16)
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
path: scipy/linalg/tests/test_matfuncs.py
text: @author:
When the input is `np.float16`, the output is `np.float64`. Is it OK?

@user1:
All float16 should be upcast to float32.

diff_hunk: @@ -412,69 +412,64 @@ def test_gh7839(self):
 
     def test_data_size_preservation_uint_in_float_out(self):
         M = np.zeros((10, 10), dtype=np.uint8)
-        # input bit size is 8, but minimum float bit size is 16
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.uint16)
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.uint32)
-        assert sqrtm(M).dtype == np.float32
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.uint64)
         assert sqrtm(M).dtype == np.float64
 
     def test_data_size_preservation_int_in_float_out(self):
         M = np.zeros((10, 10), dtype=np.int8)
-        # input bit size is 8, but minimum float bit size is 16
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.int16)
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.int32)
-        assert sqrtm(M).dtype == np.float32
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.int64)
         assert sqrtm(M).dtype == np.float64
 
     def test_data_size_preservation_int_in_comp_out(self):
         M = np.array([[2, 4], [0, -2]], dtype=np.int8)
-        # input bit size is 8, but minimum complex bit size is 64
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.int16)
-        # input bit size is 16, but minimum complex bit size is 64
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.int32)
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
         M = np.array([[2, 4], [0, -2]], dtype=np.int64)
         assert sqrtm(M).dtype == np.complex128
 
     def test_data_size_preservation_float_in_float_out(self):
         M = np.zeros((10, 10), dtype=np.float16)
-        assert sqrtm(M).dtype == np.float16
+        assert sqrtm(M).dtype == np.float64
         M = np.zeros((10, 10), dtype=np.float32)
         assert sqrtm(M).dtype == np.float32
         M = np.zeros((10, 10), dtype=np.float64)
         assert sqrtm(M).dtype == np.float64
         if hasattr(np, 'float128'):
             M = np.zeros((10, 10), dtype=np.float128)
-            assert sqrtm(M).dtype == np.float128
+            assert sqrtm(M).dtype == np.float64
 
     def test_data_size_preservation_float_in_comp_out(self):
         M = np.array([[2, 4], [0, -2]], dtype=np.float16)
-        # input bit size is 16, but minimum complex bit size is 64
-        assert sqrtm(M).dtype == np.complex64
+        assert sqrtm(M).dtype == np.complex128
path: scipy/linalg/tests/test_matfuncs.py
text: @author:
When the input is `np.float16`, the output is `np.complex128`. Is it OK?

@user1:
All float16 should be upcast to float32.