Repository: scipy/scipy
Pull request: Add _gels functions
This solve the least square problem using the qr decomposition which is faster than gelsy, gelss and gelsd which are svd based. This closes #7316 .

Problem statement:
dgels is missing
``dgels`` doesn't seem to be exposed in ``sp.linalg.lapack``. Would it make sense to expose it? Also it could be added as a driver to ``lstsq``. All the current drivers are SVD based, so it would be nice to have a method based on the QR decomposition. It should be a bit faster too. I can work on a PR if this sounds reasonable.

Hints from issue discussion:

Commit message: more complete tests

Patch to review:
diff --git a/scipy/linalg/flapack.pyf.src b/scipy/linalg/flapack.pyf.src
index c9c1788814ce..43f74e2018de 100644
--- a/scipy/linalg/flapack.pyf.src
+++ b/scipy/linalg/flapack.pyf.src
@@ -1188,6 +1188,98 @@ interface
 
    end subroutine <prefix2c>gesvd_lwork
 
+   subroutine <prefix2>gels(trans,m,n,nrhs,a,lda,b,ldb,work,lwork,info)
+
+   ! lqr,x,info = gels(a,b,lwork=..,overwrite_a=False,overwrite_b=False)
+   ! Solve Minimize 2-norm(A * X - B).
+
+     fortranname <prefix2>gels
+     callstatement (*f2py_func)(trans,&m,&n,&nrhs,a,&lda,b,&ldb,work,&lwork,&info)
+     callprotoargument char*,int*,int*,int*,<ctype2>*,int*,<ctype2>*,int*,<ctype2>*,int*,int*
+
+     character optional,intent(in),check(*trans=='N'||*trans=='T'):: trans = 'N'
+     integer depend(b),intent(hide):: nrhs = shape(b,1)
+     integer depend(a),intent(hide):: lda = shape(a,0)
+     integer intent(hide),depend(a):: m = shape(a,0)
+     integer intent(hide),depend(a):: n = shape(a,1)
+     integer depend(b),intent(hide):: ldb = shape(b,0)
+
+     <ftype2> dimension(m,nrhs),intent(in,out,copy,out=x) :: b
+     <ftype2> dimension(m,n),intent(in,out,copy,out=lqr) :: a
+     <ftype2> depend(lwork),dimension(MAX(lwork,1)),intent(hide,cache):: work
+     integer intent(out)::info
+
+   end subroutine <prefix2>gels
+
+   subroutine <prefix2>gels_lwork(trans,m,n,nrhs,a,lda,b,ldb,work,lwork,info)
+
+   ! work,info = gels_lwork(m,n,nrhs,trans='N')
+   ! Query for optimal lwork size
+
+     fortranname <prefix2>gels
+     callstatement (*f2py_func)(trans,&m,&n,&nrhs,&a,&lda,&b,&ldb,&work,&lwork,&info)
+     callprotoargument char*,int*,int*,int*,<ctype2>*,int*,<ctype2>*,int*,<ctype2>*,int*,int*
+
+     character optional,intent(in),check(*trans=='N'||*trans=='T'):: trans = 'N'
+     integer intent(in) :: m
+     integer intent(in) :: n
+     integer intent(in) :: nrhs
+     integer depend(m),intent(hide):: lda = m
+     integer depend(m),intent(hide):: ldb = m
+     integer intent(hide) :: lwork = -1
+     <ftype2> intent(hide) :: a
+     <ftype2> intent(hide) :: b
+     <ftype2> intent(out) :: work
+     integer intent(out)::info
+
+   end subroutine <prefix2>gels_lwork
+
+   subroutine <prefix2c>gels(trans,m,n,nrhs,a,lda,b,ldb,work,lwork,info)
+
+   ! lqr,x,info = gels(a,b,lwork=..,overwrite_a=False,overwrite_b=False)
+   ! Solve Minimize 2-norm(A * X - B).
+
+     fortranname <prefix2c>gels
+     callstatement (*f2py_func)(trans,&m,&n,&nrhs,a,&lda,b,&ldb,work,&lwork,&info)
+     callprotoargument char*,int*,int*,int*,<ctype2c>*,int*,<ctype2c>*,int*,<ctype2c>*,int*,int*
+
+     character optional,intent(in),check(*trans=='N'||*trans=='T'):: trans = 'N'
+     integer depend(b),intent(hide):: nrhs = shape(b,1)
+     integer depend(a),intent(hide):: lda = shape(a,0)
+     integer intent(hide),depend(a):: m = shape(a,0)
+     integer intent(hide),depend(a):: n = shape(a,1)
+     integer depend(b),intent(hide):: ldb = shape(b,0)
+
+     <ftype2c> dimension(m,n),intent(in,out,copy,out=lqr) :: a
+     <ftype2c> dimension(m,nrhs),intent(in,out,copy,out=x) :: b
+     <ftype2c> depend(lwork),dimension(MAX(lwork,1)),intent(hide,cache):: work
+     integer intent(out)::info
+
+   end subroutine <prefix2c>gels
+
+   subroutine <prefix2c>gels_lwork(trans,m,n,nrhs,a,lda,b,ldb,work,lwork,info)
+
+   ! work,info = gels_lwork(m,n,nrhs,trans='N')
+   ! Query for optimal lwork size
+
+     fortranname <prefix2c>gels
+     callstatement (*f2py_func)(trans,&m,&n,&nrhs,&a,&lda,&b,&ldb,&work,&lwork,&info)
+     callprotoargument char*,int*,int*,int*,<ctype2c>*,int*,<ctype2c>*,int*,<ctype2c>*,int*,int*
+
+     character optional,intent(in),check(*trans=='N'||*trans=='T'):: trans = 'N'
+     integer intent(in) :: m
+     integer intent(in) :: n
+     integer intent(in) :: nrhs
+     integer depend(m),intent(hide):: lda = m
+     integer depend(m),intent(hide):: ldb = m
+     integer intent(hide) :: lwork = -1
+     <ftype2c> intent(hide) :: a
+     <ftype2c> intent(hide) :: b
+     <ftype2c> intent(out) :: work
+     integer intent(out)::info
+
+   end subroutine <prefix2>gels_lwork
+
    subroutine <prefix2>gelss(m,n,minmn,maxmn,nrhs,a,b,s,cond,r,work,lwork,info)
 
    ! v,x,s,rank,work,info = gelss(a,b,cond=-1.0,overwrite_a=0,overwrite_b=0)
diff --git a/scipy/linalg/tests/test_lapack.py b/scipy/linalg/tests/test_lapack.py
index 71adb8bfbc0f..4db6a62a9f0c 100755
--- a/scipy/linalg/tests/test_lapack.py
+++ b/scipy/linalg/tests/test_lapack.py
@@ -138,6 +138,66 @@ def test_clapack(self):
 
 class TestLeastSquaresSolvers(TestCase):
 
+    def test_gels(self):
+        for dtype in REAL_DTYPES:
+            a1 = np.array([[1.0,2.0],
+                          [4.0,5.0],
+                          [7.0,8.0]], dtype=dtype)
+            b1 = np.array([16.0, 17.0, 20.0], dtype=dtype)
+            gels, gels_lwork = get_lapack_funcs(('gels','gels_lwork'),
+                                                (a1, b1))
+
+            m, n = a1.shape
+            if len(b1.shape) == 2:
+                nrhs = b1.shape[1]
+            else:
+                nrhs = 1
+
+            # Request of sizes
+            work, info = gels_lwork(m,n,nrhs)
+            lwork = int(np.real(work))
+
+            lqr, x, info = gels(a1, b1, lwork)
+            assert_allclose(x[:-1], np.array([-14.333333333333323,
+                                              14.999999999999991], dtype=dtype),
+                            rtol=25*np.finfo(dtype).eps)
+            assert_allclose(lqr, np.array([[-8.1240384 , -9.6011363 ],
+                                           [ 0.43840236,  0.90453403],
+                                           [ 0.76720414,  0.90907633]],
+                                          dtype=dtype),
+                            rtol=25*np.finfo(dtype).eps)
+
+        for dtype in COMPLEX_DTYPES:
+            a1 = np.array([[1.0+4.0j,2.0],
+                          [4.0+0.5j,5.0-3.0j],
+                          [7.0-2.0j,8.0+0.7j]], dtype=dtype)
+            b1 = np.array([16.0, 17.0+2.0j, 20.0-4.0j], dtype=dtype)
+            gels, gels_lwork = get_lapack_funcs(('gels','gels_lwork'),
+                                                (a1, b1))
+
+            m, n = a1.shape
+            if len(b1.shape) == 2:
+                nrhs = b1.shape[1]
+            else:
+                nrhs = 1
+
+            # Request of sizes
+            work, info = gels_lwork(m,n,nrhs)
+            lwork = int(np.real(work))
+
+            lqr, x, info = gels(a1, b1, lwork)
+            assert_allclose(x[:-1],
+                            np.array([1.161753632288328-1.901075709391912j,
+                                      1.735882340522193+1.521240901196909j],
+                            dtype=dtype), rtol=25*np.finfo(dtype).eps)
+            assert_allclose(lqr,
+                            np.array([[-9.28708781+0.j, -8.08649617+0.17228221j],
+                                      [ 0.35418546-0.08911578j, -6.08842332+0.j],
+                                      [ 0.52542621-0.39872361j, 0.18127551+0.68807219j]],
+                                     dtype=dtype),
+                            rtol=25*np.finfo(dtype).eps)
+
+
     def test_gelsd(self):
         for dtype in REAL_DTYPES:
             a1 = np.array([[1.0,2.0],
@@ -343,7 +403,7 @@ def test_gh_2691(self):
                     assert_allclose(np.tril(dpt), np.tril(inv(a)))
                 else:
                     assert_allclose(np.triu(dpt), np.triu(inv(a)))
-                    
+
 class TestDlasd4(TestCase):
     def test_sing_val_update(self):
 
@@ -425,7 +485,7 @@ def test_rot():
         assert_allclose(rot(u, v, c, s, offx=2, incy=2, n=2), [[3,3,5,5],[0,f,0,f]], atol=atol)
         assert_allclose(rot(u, v, c, s, offx=2, incx=2, offy=2, incy=2, n=1), [[3,3,5,3],[f,f,0,f]], atol=atol)
         assert_allclose(rot(u, v, c, s, incx=-2, incy=-2, n=2), [[5,3,5,3],[0,f,0,f]], atol=atol)
-    
+
         a, b = rot(u, v, c, s, overwrite_x=1, overwrite_y=1)
         assert_(a is u)
         assert_(b is v)
@@ -458,7 +518,7 @@ def test_larfg_larf():
         expected = np.zeros_like(a[:,0])
         expected[0] = a[0,0]
         expected[1] = alpha
-        
+
         # assemble householder vector
         v = np.zeros_like(a[1:,0])
         v[0] = 1.0
@@ -469,7 +529,7 @@ def test_larfg_larf():
 
         # apply transform from the right
         a[:,1:] = larf(v, tau, a[:,1:], np.zeros(a.shape[0]), side='R')
-        
+
         assert_allclose(a[:,0], expected, atol=1e-5)
         assert_allclose(a[0,:], expected, atol=1e-5)
 


Code review comments:
diff_hunk: @@ -138,6 +138,66 @@ def test_clapack(self):
 
 class TestLeastSquaresSolvers(TestCase):
 
+    def test_gels(self):
+        for dtype in REAL_DTYPES:
+            a1 = np.array([[1.0,2.0],
+                          [4.0,5.0],
+                          [7.0,8.0]], dtype=dtype)
+            b1 = np.array([16.0, 17.0, 20.0], dtype=dtype)
+            gels, gels_lwork = get_lapack_funcs(('gels','gels_lwork'),
+                                                (a1, b1))
+
+            m, n = a1.shape
+            if len(b1.shape) == 2:
+                nrhs = b1.shape[1]
+            else:
+                nrhs = 1
+
+            # Request of sizes
+            work, info = gels_lwork(m,n,nrhs)
+            lwork = int(np.real(work))
+
+            lqr, x, info = gels(a1, b1, lwork)
+            assert_allclose(x[:-1], np.array([-14.333333333333323,
+                                              14.999999999999991], dtype=dtype),
+                            rtol=25*np.finfo(dtype).eps)
+            assert_allclose(lqr, np.array([[-8.1240384 , -9.6011363 ],
+                                           [ 0.43840236,  0.90453403],
+                                           [ 0.76720414,  0.90907633]],
+                                          dtype=dtype),
+                            rtol=25*np.finfo(dtype).eps)
path: scipy/linalg/tests/test_lapack.py
text: @user1:
https://travis-ci.org/scipy/scipy/jobs/227275929 looks real. On my machine the abs difference is of the order of 1e-9.

diff_hunk: @@ -138,6 +138,66 @@ def test_clapack(self):
 
 class TestLeastSquaresSolvers(TestCase):
 
+    def test_gels(self):
+        for dtype in REAL_DTYPES:
+            a1 = np.array([[1.0,2.0],
+                          [4.0,5.0],
+                          [7.0,8.0]], dtype=dtype)
+            b1 = np.array([16.0, 17.0, 20.0], dtype=dtype)
+            gels, gels_lwork = get_lapack_funcs(('gels','gels_lwork'),
+                                                (a1, b1))
+
+            m, n = a1.shape
+            if len(b1.shape) == 2:
+                nrhs = b1.shape[1]
+            else:
+                nrhs = 1
+
+            # Request of sizes
+            work, info = gels_lwork(m,n,nrhs)
path: scipy/linalg/tests/test_lapack.py
text: @user1:
PEP8 says spaces after commas please. Also several other PEP8 violations, https://travis-ci.org/scipy/scipy/jobs/227275928

@user2:
Minor comment: I'd maybe rewrite this to check for `lwork = _compute_lwork(gels_lwork, m, n, nrhs)` where that routine comes from `from scipy.linalg.lapack import _compute_lwork`