Repository: scipy/scipy
Pull request: Do not segfault in svd(a) with VT.size > INT_MAX
#### Reference issue
<!--Example: Closes gh-WXYZ.-->

closes https://github.com/scipy/scipy/issues/14001

#### What does this implement/fix?
<!--Please explain your changes.-->

`linalg.svd` with too large matrices may segfault if `m*n > int_max` on large-memory machines (on smaller memory machines it may fail with a MemoryError instead). The root cause is an integer overflow in indexing 2D arrays, deep in the LAPACK code.

Thus, detect a possible error condition in the f2py wrapper and error out early.

#### Additional information
<!--Any additional information you think is important.-->

Suggested by @pearu in https://github.com/scipy/scipy/issues/14001#issuecomment-2015872894

Adding this kind of check to f2py may still make sense IMO. No reason to not guard against this specific segfault in the meantime though.

Problem statement:
Pycharm scipy SVD returning error code without message
I am trying to compute the SVD of a large matrix (4799 x 53130) using Pycharm as my IDE. The following returns the exit code -1073741819 (0xC0000005), but no error message is displayed. I have looked up this error code on here and people have found a mixture of things causing it. I have re-installed Pycharm and all dependencies, yet this issue persists. I have tried also using the numpy version, which returns the same error. I'm confident this problem is related to the lapack package, since both use it. Any direction would be greatly appreciated. 

    import numpy as np
    import pandas as pd
    from scipy.linalg import svd
    df = pd.read_pickle(r'filepath\data.pkl')
    df = df.to_numpy()
    svd(df) 

I tried using df=np.ones((4799,53130)) too in case it was my data that was the problem, and this yields the same issue. It isn't the amount of RAM, my machine has 120GB and the process doesn't go near the limit.

Hints from issue discussion:
You probably run out of 32bit indexing. Could you try building both using 64bit Lapack known as ilp64? Here is the first hit after a search about mkl and openblas 

https://github.com/pytorch/pytorch/issues/51720
> You probably run out of 32bit indexing. Could you try building both using 64bit Lapack known as ilp64? Here is the first hit after a search about mkl and openblas
> 
> [pytorch/pytorch#51720](https://github.com/pytorch/pytorch/issues/51720)

Hi Ilayn, what do you mean by building both using 64bit? Sorry, this is above my comprehension.
The default underlying lapack library uses 32 bit integers internally for indexing arrays and for large arrays they might overflow. There is also a 64bit integer version of lapack that can be used to build numpy and scipy. That's what is being mentioned in the linked issue. 
> 64bit integer version of lapack

I looked in the linalg folder of scipy, and found the lapack code that is called, which has the following:
![image](https://user-images.githubusercontent.com/60707891/117434555-b973f580-af24-11eb-9401-428adbed3a85.png)
It seems to already have the 64bit version? 

EDIT: I think from the linked thread it's not the Lapack function that's the problem, it's the scipy package using only the 32-bit version and the OP in that linked thread is requesting 64-bit support from scipy, is this correct? If so I think I will have to wait it out and see if some smart folks figure it out. Do you know if there is anyway to catch the error code and not terminate the program? I tried using try and except but that doesn't work, I think exit or quit is called in one of the packages. 
I can reproduce a scipy.linalg.svd crash with the following gdb traceback:
```
>>> import scipy.linalg
>>> import numpy as np
>>> df=np.ones((4799,53130), dtype=np.float64)
>>> r=scipy.linalg.svd(df)
```

<details>

```
>>> r=scipy.linalg.svd(df)

Thread 1 "python" received signal SIGSEGV, Segmentation fault.
0x00007ffff5d82340 in dlarf_ () from /home/pearu/miniconda3/envs/pytorch-cuda-dev/lib/python3.9/site-packages/numpy/core/../../../../libcblas.so.3
(gdb) bt
#0  0x00007ffff5d82340 in dlarf_ () from /home/pearu/miniconda3/envs/pytorch-cuda-dev/lib/python3.9/site-packages/numpy/core/../../../../libcblas.so.3
#1  0x00007ffff5d9a806 in dorgl2_ () from /home/pearu/miniconda3/envs/pytorch-cuda-dev/lib/python3.9/site-packages/numpy/core/../../../../libcblas.so.3
#2  0x00007ffff5d9aade in dorglq_ () from /home/pearu/miniconda3/envs/pytorch-cuda-dev/lib/python3.9/site-packages/numpy/core/../../../../libcblas.so.3
#3  0x00007ffff5d30761 in dgesdd_ () from /home/pearu/miniconda3/envs/pytorch-cuda-dev/lib/python3.9/site-packages/numpy/core/../../../../libcblas.so.3
#4  0x00007ffec9d2d7e3 in f2py_rout.flapack_dgesdd () from /home/pearu/miniconda3/envs/pytorch-cuda-dev/lib/python3.9/site-packages/scipy/linalg/_flapack.cpython-39-x86_64-linux-gnu.so
#5  0x0000555555699456 in _PyObject_MakeTpCall () at /home/conda/feedstock_root/build_artifacts/python-split_1610246514195/work/Include/internal/pycore_pyerrors.h:14
<snip>
```

</details>

Note that lapack routines use `integer` type for indices that corresponds to 4-byte int type. The size of dgesdd work array is
`LWORK = 4*mn*mn + 6*mn + mx = 92203528` which is less than `2**27`, and `4799*53130 < 2**28`, so, in theory, there should not be integer overflows when the inputs have the given size.

Btw, the following works (while using numpy here, it has the same crash message as scipy when using `compute_uv=True` which is default):
```
>>> df=np.ones((4799, 53130))
>>> r=np.linalg.svd(df, compute_uv=False)
```
that is, if only singular values are required, disable computing the U and V matrices.
If you know that the input matrix has low rank, try using lowrank specific svd algorithms (such as torch.svd_lowrank).
Thanks for the response Pearu, unfortunately I need the full decomposition. The bigger problem is that I cannot use try and except statements, the call of the SVD function actually causes my program to exit. It would be better if an exception was raised with some form of message. 

Anyway, thanks for all your help everyone,

Joel. 
@pearu noted the absolute minimums about the workspace but it might happen that gesdd internally might ask for more. Could you please try the same with converting your array to float32 data type just to make sure? That would mean considerable less memory but still the same indexing constraint.
> @pearu noted the absolute minimums about the workspace but it might happen that gesdd internally might ask for more. Could you please try the same with converting your array to float32 data type just to make sure? That would mean considerable less memory but still the same indexing constraint.

The same behaviour exists when using float32 inputs. Re "gesdd internally might ask for more", AFAIK, gesdd does not manage memory.

Another cause of this issue could be an uncatched integer overflow that might lead to wrong memory access conditions but I cannot confirm this atm. One would need to run the test for an increasing number of sizes and see at which point the crash will occur and debug this parameter case in detail.
what I mean is that we are querying for lwork and that number returned can go up from that minimum see the [Notes section and the inside code comments](http://www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_gad8e0f1c83a78d3d4858eaaa88a1c5ab1.html)
Running Pearu's example on a machine with 16Gb of memory generates a MemoryError:

```
$ cat svd_oops.py 
import scipy.linalg
import numpy as np
df=np.ones((4799,53130), dtype=np.float64)
breakpoint()
scipy.linalg.svd(df)
br@gonzales:~/sweethome/temp/chol$ python svd_oops.py 
> /home/br/sweethome/temp/chol/svd_oops.py(5)<module>()
-> scipy.linalg.svd(df)
(Pdb) c
Traceback (most recent call last):
  File "/home/br/sweethome/temp/chol/svd_oops.py", line 5, in <module>
    scipy.linalg.svd(df)
  File "/home/br/repos/scipy/scipy/build-install/lib/python3.10/site-packages/scipy/linalg/_decomp_svd.py", line 127, in svd
    u, s, v, info = gesXd(a1, compute_uv=compute_uv, lwork=lwork,
numpy.core._exceptions._ArrayMemoryError: Unable to allocate 21.0 GiB for an array with shape (53130, 53130) and data type float64
```

dropping into pdb shows that lwork is OK, but the the unitary matrices in the decomposition are too large:

```
(Pdb) lwork
92155197
(Pdb) 53130**2
2822796900
(Pdb) np.iinfo(np.int32).max
2147483647
```

On a larger-memory machine (not sure how much, definitely more than 128Gb), the snippet above does segfault.

To double-check that that it's indeed `U` and `VT` matrices, changing `scipy.linalg.svd(df)` into `scipy.linalg.svd(df, compute_uv=False)` does not segfault and returns a reasonable result:

```
ev-br@qgpu3:~/temp$ cat svd_oops.py 
import scipy.linalg
import numpy as np
df=np.ones((4799,53130), dtype=np.float64)
u, s, vt  = scipy.linalg.svd(df, full_matrices=False)
print( (u @ np.diag(s) @ vt - df).max())
ev-br@qgpu3:~/temp$ python svd_oops.py 
1.1737300020797647e-10
```

Since `U` and `VT` matrices are generated internally by f2py, it looks like the issue is somewhere in f2py.

Commit message: BUG: linalg: avoid a segfault in ?gesdd

linalg.svd with too large matrices may segfault if
max(m, n)*max(m, n) > int_max on large-memory machines (on smaller memory machines
it may fail with a MemoryError instead). The root cause is an integer overflow
in indexing 2D arrays, deep in the LAPACK code.
Thus, detect a possible error condition and bail out early.

Patch to review:
diff --git a/scipy/linalg/_decomp_svd.py b/scipy/linalg/_decomp_svd.py
index 24be62752e62..5f8bdaac8a37 100644
--- a/scipy/linalg/_decomp_svd.py
+++ b/scipy/linalg/_decomp_svd.py
@@ -116,6 +116,20 @@ def svd(a, full_matrices=True, compute_uv=True, overwrite_a=False,
     if lapack_driver not in ('gesdd', 'gesvd'):
         message = f'lapack_driver must be "gesdd" or "gesvd", not "{lapack_driver}"'
         raise ValueError(message)
+
+    if lapack_driver == 'gesdd' and compute_uv:
+        # XXX: revisit int32 when ILP64 lapack becomes a thing
+        max_mn, min_mn = (m, n) if m > n else (n, m)
+        if full_matrices:
+            if max_mn*max_mn > numpy.iinfo(numpy.int32).max:
+                raise ValueError(f"Indexing a matrix size {max_mn} x {max_mn} "
+                                 f" would incur integer overflow in LAPACK.")
+        else:
+            sz = max(m * min_mn, n * min_mn)
+            if max(m * min_mn, n * min_mn) > numpy.iinfo(numpy.int32).max:
+                raise ValueError(f"Indexing a matrix of {sz} elements would "
+                                  "incur an in integer overflow in LAPACK.")
+
     funcs = (lapack_driver, lapack_driver + '_lwork')
     gesXd, gesXd_lwork = get_lapack_funcs(funcs, (a1,), ilp64='preferred')
 
diff --git a/scipy/linalg/tests/test_decomp.py b/scipy/linalg/tests/test_decomp.py
index bcd5b6611d80..f9e347050627 100644
--- a/scipy/linalg/tests/test_decomp.py
+++ b/scipy/linalg/tests/test_decomp.py
@@ -1089,6 +1089,14 @@ class TestSVD_GESVD(TestSVD_GESDD):
     lapack_driver = 'gesvd'
 
 
+def test_svd_gesdd_nofegfault():
+    # svd(a) with {U,VT}.size > INT_MAX does not segfault
+    # cf https://github.com/scipy/scipy/issues/14001
+    df=np.ones((4799, 53130), dtype=np.float64)
+    with assert_raises(ValueError):
+        svd(df)
+
+
 class TestSVDVals:
 
     def test_empty(self):


Code review comments:
diff_hunk: @@ -116,6 +116,20 @@ def svd(a, full_matrices=True, compute_uv=True, overwrite_a=False,
     if lapack_driver not in ('gesdd', 'gesvd'):
         message = f'lapack_driver must be "gesdd" or "gesvd", not "{lapack_driver}"'
         raise ValueError(message)
+
+    if lapack_driver == 'gesdd' and compute_uv:
+        # XXX: revisit int32 when ILP64 lapack becomes a thing
+        max_mn, min_mn = (m, n) if m > n else (n, m)
+        if full_matrices:
+            if max_mn*max_mn > numpy.iinfo(numpy.int32).max:
+                raise ValueError(f"Indexing a matrix size {max_mn} x {max_mn} "
+                                 f" would incur integer overflow in LAPACK.")
path: scipy/linalg/_decomp_svd.py
text: @user1:
Nit:
```suggestion
                                  " would incur integer overflow in LAPACK.")
```