Repository: cupy/cupy
Pull request: Update doctest of `cupyx.scipy.linalg.lu_factor`
Close #4259.

It looks that `cupyx.scipy.linalg.lu_factor` returns the same result to `scipy.linalg.lu_factor`. I remove the note on the document and add a test given a singular matrix.

Problem statement:
Doctest of `lu_factor` is not valid with newer CUDA
### Conditions

<details><summary>the output of `python -c 'import cupy; cupy.show_config()'`</summary>

```
OS                           : Linux-4.15.0-58-generic-x86_64-with-glibc2.27
CuPy Version                 : 9.0.0a1
NumPy Version                : 1.19.2
SciPy Version                : 1.5.2
Cython Build Version         : 0.29.21
Cython Runtime Version       : 0.29.21
CUDA Root                    : /usr/local/cuda
CUDA Build Version           : 10010
CUDA Driver Version          : 11010
CUDA Runtime Version         : 10010
cuBLAS Version               : 10201
cuFFT Version                : 10101
cuRAND Version               : 10101
cuSOLVER Version             : (10, 2, 0)
cuSPARSE Version             : 10300
NVRTC Version                : (10, 1)
Thrust Version               : 100906
CUB Build Version            : 100800
cuDNN Build Version          : 7500
cuDNN Version                : 7500
NCCL Build Version           : 2703
NCCL Runtime Version         : 2703
cuTENSOR Version             : 10200
Device 0 Name                : Tesla V100-SXM2-16GB
Device 0 Compute Capability  : 70
Device 1 Name                : Tesla V100-SXM2-16GB
Device 1 Compute Capability  : 70
Device 2 Name                : Tesla V100-SXM2-16GB
Device 2 Compute Capability  : 70
Device 3 Name                : Tesla V100-SXM2-16GB
Device 3 Compute Capability  : 70
Device 4 Name                : Tesla V100-SXM2-16GB
Device 4 Compute Capability  : 70
Device 5 Name                : Tesla V100-SXM2-16GB
Device 5 Compute Capability  : 70
Device 6 Name                : Tesla V100-SXM2-16GB
Device 6 Compute Capability  : 70
Device 7 Name                : Tesla V100-SXM2-16GB
Device 7 Compute Capability  : 70
```

</details>

### Code to reproduce

`cd docs; make doctest`

### Error messages, stack traces, or logs

```
Document: reference/generated/cupyx.scipy.linalg.lu_factor
----------------------------------------------------------
<doctest default[2]>:1: LinAlgWarning: Diagonal number 1 is exactly zero. Singular matrix.
  scipy.linalg.lu_factor(np.array([[0, 1], [0, 0]], dtype=np.float32))
/mnt/vol21/kataoka/github/cupy/cupyx/scipy/linalg/decomp_lu.py:59: RuntimeWarning: Diagonal number 1 is exactly zero. Singular matrix.
  return _lu_factor(a, overwrite_a, check_finite)
**********************************************************************
File "reference/generated/cupyx.scipy.linalg.lu_factor.rst", line ?, in default
Failed example:
    cupyx.scipy.linalg.lu_factor(cp.array([[0, 1], [0, 0]], dtype=cp.float32))
Expected:
    (array([[ 0.,  1.],
           [nan, nan]], dtype=float32), array([0, 1], dtype=int32))
Got:
    (array([[0., 1.],
           [0., 0.]], dtype=float32), array([ 0, -1], dtype=int32))
**********************************************************************
1 items had failures:
   1 of   6 in default
6 tests in 1 items.
5 passed and 1 failed.
***Test Failed*** 1 failures.
```

Hints from issue discussion:

Commit message: Add test for lu_factor given singular matrix

Patch to review:
diff --git a/cupyx/scipy/linalg/decomp_lu.py b/cupyx/scipy/linalg/decomp_lu.py
index 537b70408e5..96a1700e80d 100644
--- a/cupyx/scipy/linalg/decomp_lu.py
+++ b/cupyx/scipy/linalg/decomp_lu.py
@@ -35,26 +35,6 @@ def lu_factor(a, overwrite_a=False, check_finite=True):
             ``i`` of the matrix was interchanged with row ``piv[i]``
 
     .. seealso:: :func:`scipy.linalg.lu_factor`
-
-    .. note::
-
-        Current implementation returns result different from SciPy when the
-        matrix singular. SciPy returns an array containing ``0.`` while the
-        current implementation returns an array containing ``nan``.
-
-        >>> import numpy as np
-        >>> import scipy.linalg
-        >>> scipy.linalg.lu_factor(np.array([[0, 1], [0, 0]], \
-dtype=np.float32))
-        (array([[0., 1.],
-               [0., 0.]], dtype=float32), array([0, 1], dtype=int32))
-
-        >>> import cupy as cp
-        >>> import cupyx.scipy.linalg
-        >>> cupyx.scipy.linalg.lu_factor(cp.array([[0, 1], [0, 0]], \
-dtype=cp.float32))
-        (array([[ 0.,  1.],
-               [nan, nan]], dtype=float32), array([0, 1], dtype=int32))
     """
     return _lu_factor(a, overwrite_a, check_finite)
 
diff --git a/tests/cupyx_tests/scipy_tests/linalg_tests/test_decomp_lu.py b/tests/cupyx_tests/scipy_tests/linalg_tests/test_decomp_lu.py
index 6d09d1760b3..db3c7645eaa 100644
--- a/tests/cupyx_tests/scipy_tests/linalg_tests/test_decomp_lu.py
+++ b/tests/cupyx_tests/scipy_tests/linalg_tests/test_decomp_lu.py
@@ -59,6 +59,17 @@ def test_lu_factor_reconstruction(self, dtype):
         cupy.testing.assert_allclose(LU, PA, atol=1e-5)
 
 
+@testing.gpu
+@testing.with_requires('scipy')
+class TestLUFactorSingular:
+
+    @testing.for_dtypes('fdFD')
+    @testing.numpy_cupy_allclose(scipy_name='scp')
+    def test_lu_factor_singular(self, xp, scp, dtype):
+        a = xp.array([[0, 1], [0, 0]], dtype=dtype)
+        return scp.linalg.lu_factor(a)
+
+
 @testing.gpu
 @testing.parameterize(*testing.product({
     'shape': [(1, 1), (2, 2), (3, 3), (5, 5), (1, 5), (5, 1), (2, 5), (5, 2)],


Code review comments:
diff_hunk: @@ -59,6 +59,17 @@ def test_lu_factor_reconstruction(self, dtype):
         cupy.testing.assert_allclose(LU, PA, atol=1e-5)
 
 
+@testing.gpu
+@testing.with_requires('scipy')
+class TestLUFactorSingular:
+
+    @testing.for_dtypes('fdFD')
+    @testing.numpy_cupy_allclose(scipy_name='scp')
path: tests/cupyx_tests/scipy_tests/linalg_tests/test_decomp_lu.py
text: @user1:
`numpy_cupy_allclose` is OK just for the matrix `[[0, 1], [0, 0]]` (rank=n-1 and no numerical error).

```
a = rng.normal(size=(100, 50)) @ rng.normal(size=(50, 100))  # 100x100, rank-50
>>> scipy.linalg.lu_factor(a)[1]
array([26, 34,  9, 37, 29, 73, 35, 99, 14, 46, 15, 13, 89, 56, 70, 46, 98,
       63, 88, 84, 71, 21, 76, 68, 41, 57, 29, 39, 80, 39, 60, 66, 60, 49,
       78, 69, 65, 72, 50, 56, 78, 81, 90, 73, 92, 94, 85, 82, 77, 62, 87,
       93, 65, 97, 57, 86, 95, 71, 96, 79, 61, 86, 67, 73, 75, 74, 85, 79,
       85, 97, 73, 93, 89, 82, 82, 82, 89, 86, 84, 94, 89, 98, 86, 94, 99,
       90, 95, 95, 96, 98, 93, 95, 92, 98, 98, 99, 96, 99, 98, 99],
      dtype=int32)
>>> cupyx.scipy.linalg.lu_factor(cupy.array(a))[1]
array([26, 34,  9, 37, 29, 73, 35, 99, 14, 46, 15, 13, 89, 56, 70, 46, 98,
       63, 88, 84, 71, 21, 76, 68, 41, 57, 29, 39, 80, 39, 60, 66, 60, 49,
       78, 69, 65, 72, 50, 56, 78, 81, 90, 73, 92, 94, 85, 82, 77, 62, 64,
       97, 59, 71, 94, 67, 94, 61, 96, 82, 89, 91, 71, 96, 65, 85, 75, 90,
       77, 95, 76, 84, 90, 76, 91, 99, 82, 78, 87, 85, 81, 86, 93, 97, 93,
       96, 98, 96, 96, 95, 93, 94, 97, 94, 99, 98, 99, 99, 99, 99],
      dtype=int32)
```

Comment why the test is OK or test only the reconstructed matrix like `test_lu_factor_reconstruction` above.

To generate a singular matrix, use `cupy.testing.helper.generate_matrix` or follow https://github.com/cupy/cupy/blob/ba8c3fff19b322d5392311312d944b66c84cddf5/tests/cupy_tests/linalg_tests/test_solve.py#L194-L197

@author:
Thanks! I'll fix the test as you show the right way.

@author:
PTAL, @user1! BTW, can I have why the operation above brings a matrix singular?

@user1:
Multiplying the `ones` vector (from left or from right) results zero.  Therefore, `matrix_rank(a) <= min(m-1, n-1) < min(m, n)`.