Repository: scipy/scipy
Pull request: MAINT:linalg: Make lu_factor accept rectangular arrays
#### Reference issue
closes #14243 

#### What does this implement/fix?
Removes check that matrix is square for `linalg.lu_factor`. This function is valid for rectangular matrices as well.
The check was probably there because `lu_solve` only works for square matrices, but `lu_factor` can be used on its own. (My personal reason for doing so is the fact that the permutation matrix is returned as a dense 2D array in `lu`, whereas `lu_factor` returns a 1D integer array of pivot indices containing the same information)

The docs for `lu_factor` are also improved, better highlighting the differences to `lu`.

#### Additional information
I took the same test cases for `lu_factor` as for `lu`, except for 4 'simple' matrices. These are singular, and triggers the warning in `lu_factor` defined here:
```py
    if info > 0:
        warn("Diagonal number %d is exactly zero. Singular matrix." % info,
             LinAlgWarning, stacklevel=2)
```
I think in this particular case the warning can be ignored, but instead I opted to not run the test on these particular matrices. 

Problem statement:
Sparse pivot for linalg.lu
Currently `p, l, u = scipy.linalg.lu(a)` returns a dense permutation matrix `p` for the pivots of shape `(M, M)` for `a` of shape `(M, N)`. If the input matrix is tall and skinny (`M` >> `N`), then the storage costs of `p` are unacceptable. 

A workaround is to make LAPACK calls directly, e.g. use `linalg.lapack.dgetrf(a)`, since the LAPACK implementation of the LU factorization returns a 1-d integer array for the pivots. 

It seems that `linalg.lu` uses `linalg/src/lu.f` rather than the LAPACK implementation. Is this faster than the LAPACK implementation? Either way, I suggest to add an option to return an 1-d integer array for the pivots. If this option is invoked, then the appropriate LAPACK routine is called. This should be very easy to implement, shall I make a PR?

Hints from issue discussion:
If you only need the permutation you can use `scipy.linalg.lu_factor` which is a bit more labourous to use the L, U matrices but it is meant to be used with `scipy.linalg.lu_solve`
`scipy.linalg.lu_factor` is only for square matrices. My use case is precisely tall skinny matrices.
Also, it seems `dgetrf` is faster than `scipy.linalg.lu`, especially for tall matrices. For square matrices (in simple benchmarks) it is a bit faster, but not by much. 
That is a bug if `lu_factor` only works with square arrays and we have to fix it. dgetrs should require squareness not the trf
Right, `lu_factor` calls `getrf`, but it explicitly checks for a square matrix,
I would suggest removing the the check, and updating the documentation. It should also be a little clearer in the documentation what the difference is between `lu_factor` and `lu`. 
For `lu_factor` one can say something along the likes of
"Unlike `scipy.linalg.lu`, this function returns L and U stored in the same matrix LU, with U the upper triangular part, and L the lower triangular part with ones on the diagonal. Furthermore the permutation matrix P is returned as a 1D integer array instead of a  dense matrix"

Shall I make a PR?
That would much be appreciated if you can indeed. Thanks @RikVoorhaar 
As mentioned in my original post
> It seems that `linalg.lu` uses `linalg/src/lu.f` rather than the LAPACK implementation. Is this faster than the LAPACK implementation?

It seems that it is in fact slower than the LAPACK implementation in my own benchmarks. For square matrices the methods are comparable with the Fortran `lu.f` being slightly faster, but for matrices of shape (M,N) with M>>N the LAPACK implementation is faster. 
In extreme cases of a (10,10000) matrix the difference between the two methods is a factor 500.

I can past a more detailed benchmark, (is there a specific format I should adhere to?).

Based on detailed benchmarks we can make `scipy.linalg.lu` call either the Fortran implementation or call LAPACK depending on the shape of the matrix.
Wow I've been looking at this codebase for quite some time now and still get surprises. I just realized that this `lu.f` file is still in use. This probably predates the LAPACK wrappers and hence the custom roll of 4 flavors. As you said we should definitely switch back to regular wrappers unconditionally because there is no need for it. The overhead probably would be negligible depending on the matrix size


Shall I make a PR for a version using a LAPACK wrapper?

Commit message: Merge remote-tracking branch 'upstream/main'

Patch to review:
diff --git a/scipy/linalg/_decomp_lu.py b/scipy/linalg/_decomp_lu.py
index 3375409c30d0..ac1a37197401 100644
--- a/scipy/linalg/_decomp_lu.py
+++ b/scipy/linalg/_decomp_lu.py
@@ -25,7 +25,7 @@ def lu_factor(a, overwrite_a=False, check_finite=True):
 
     Parameters
     ----------
-    a : (M, M) array_like
+    a : (M, N) array_like
         Matrix to decompose
     overwrite_a : bool, optional
         Whether to overwrite data in A (may increase performance)
@@ -36,7 +36,7 @@ def lu_factor(a, overwrite_a=False, check_finite=True):
 
     Returns
     -------
-    lu : (N, N) ndarray
+    lu : (M, N) ndarray
         Matrix containing U in its upper triangle, and L in its lower triangle.
         The unit diagonal elements of L are not stored.
     piv : (N,) ndarray
@@ -45,11 +45,14 @@ def lu_factor(a, overwrite_a=False, check_finite=True):
 
     See Also
     --------
+    lu : gives lu factorization in more user-friendly format
     lu_solve : solve an equation system using the LU factorization of a matrix
 
     Notes
     -----
-    This is a wrapper to the ``*GETRF`` routines from LAPACK.
+    This is a wrapper to the ``*GETRF`` routines from LAPACK. Unlike
+`linalg.lu`, It outputs the L and U factors into a single array,
+and returns pivot indices instead of a permutation matrix.
 
     Examples
     --------
@@ -71,8 +74,6 @@ def lu_factor(a, overwrite_a=False, check_finite=True):
         a1 = asarray_chkfinite(a)
     else:
         a1 = asarray(a)
-    if len(a1.shape) != 2 or (a1.shape[0] != a1.shape[1]):
-        raise ValueError('expected square matrix')
     overwrite_a = overwrite_a or (_datacopied(a1, a))
     getrf, = get_lapack_funcs(('getrf',), (a1,))
     lu, piv, info = getrf(a1, overwrite_a=overwrite_a)
diff --git a/scipy/linalg/tests/test_decomp.py b/scipy/linalg/tests/test_decomp.py
index 77a24c0af5dd..ec5c5c97e07f 100644
--- a/scipy/linalg/tests/test_decomp.py
+++ b/scipy/linalg/tests/test_decomp.py
@@ -10,8 +10,6 @@
 
 import itertools
 import platform
-import sys
-
 import numpy as np
 from numpy.testing import (assert_equal, assert_almost_equal,
                            assert_array_almost_equal, assert_array_equal,
@@ -28,7 +26,8 @@
                           eigh_tridiagonal, null_space, cdf2rdf, LinAlgError)
 
 from scipy.linalg.lapack import (dgbtrf, dgbtrs, zgbtrf, zgbtrs, dsbev,
-                                 dsbevd, dsbevx, zhbevd, zhbevx)
+                                 dsbevd, dsbevx, zhbevd, zhbevx,
+                                 get_lapack_funcs)
 
 from scipy.linalg._misc import norm
 from scipy.linalg._decomp_qz import _select_function
@@ -941,7 +940,15 @@ def _test_common(self, data):
         pl, u = lu(data, permute_l=1)
         assert_array_almost_equal(pl @ u, data)
 
-    # Simple tests
+    def _test_common_lu_factor(self, data):
+        l_and_u1, piv1 = lu_factor(data)
+        (getrf,) = get_lapack_funcs(("getrf",), (data,))
+        l_and_u2, piv2, _ = getrf(data, overwrite_a=False)
+        assert_array_equal(l_and_u1, l_and_u2)
+        assert_array_equal(piv1, piv2)
+
+    # Simple tests.
+    # For lu_factor gives a LinAlgWarning because these matrices are singular
     def test_simple(self):
         self._test_common(self.a)
 
@@ -957,24 +964,30 @@ def test_simple2_complex(self):
     # rectangular matrices tests
     def test_hrectangular(self):
         self._test_common(self.hrect)
+        self._test_common_lu_factor(self.hrect)
 
     def test_vrectangular(self):
         self._test_common(self.vrect)
+        self._test_common_lu_factor(self.vrect)
 
     def test_hrectangular_complex(self):
         self._test_common(self.chrect)
+        self._test_common_lu_factor(self.chrect)
 
     def test_vrectangular_complex(self):
         self._test_common(self.cvrect)
+        self._test_common_lu_factor(self.cvrect)
 
     # Bigger matrices
     def test_medium1(self):
         """Check lu decomposition on medium size, rectangular matrix."""
         self._test_common(self.med)
+        self._test_common_lu_factor(self.med)
 
     def test_medium1_complex(self):
         """Check lu decomposition on medium size, rectangular matrix."""
         self._test_common(self.cmed)
+        self._test_common_lu_factor(self.cmed)
 
     def test_check_finite(self):
         p, l, u = lu(self.a, check_finite=False)
@@ -2033,8 +2046,6 @@ class TestQZ:
     def setup_method(self):
         seed(12345)
 
-    @pytest.mark.xfail(sys.platform == 'darwin',
-        reason="gges[float32] broken for OpenBLAS on macOS, see gh-16949")
     def test_qz_single(self):
         n = 5
         A = random([n, n]).astype(float32)


Code review comments:
diff_hunk: @@ -45,11 +45,14 @@ def lu_factor(a, overwrite_a=False, check_finite=True):
 
     See Also
     --------
+    lu : gives lu factorization in more user-friendly format
     lu_solve : solve an equation system using the LU factorization of a matrix
 
     Notes
     -----
-    This is a wrapper to the ``*GETRF`` routines from LAPACK.
+    This is a wrapper to the ``*GETRF`` routines from LAPACK. Unlike
+`linalg.lu`, It outputs the L and U factors into a single array,
+and returns pivot indices instead of a permutation matrix.
path: scipy/linalg/_decomp_lu.py
text: @user1:
```suggestion
    `linalg.lu`, It outputs the L and U factors into a single array,
    and returns pivot indices instead of a permutation matrix.
```