3 files changed, 183 insertions, 130 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 3d4a209..de88dd9 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -27,6 +27,7 @@ module Numeric.LinearAlgebra.Algorithms (
 -- * Linear Systems
     linearSolve,
     luSolve,
+    linearSolveLS,
     linearSolveSVD,
     inv, pinv,
     det, rank, rcond,
@@ -36,13 +37,13 @@ module Numeric.LinearAlgebra.Algorithms (
     fullSVD,
     thinSVD,
     compactSVD,
-    sv,
+    singularValues,
     leftSV, rightSV,
 -- ** Eigensystems
     eig, eigSH, eigSH',
     eigenvalues, eigenvaluesSH, eigenvaluesSH',
 -- ** QR
-    qr,
+    qr, rq,
 -- ** Cholesky
     chol, cholSH,
 -- ** Hessenberg
@@ -91,6 +92,7 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
     luSolve'     :: (Matrix t, [Int]) -> Matrix t -> Matrix t
     linearSolve' :: Matrix t -> Matrix t -> Matrix t
     linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
+    linearSolveLS'  :: Matrix t -> Matrix t -> Matrix t
     eig'         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
     eigSH''      :: Matrix t -> (Vector Double, Matrix t)
     eigOnly      :: Matrix t -> Vector (Complex Double)
@@ -103,35 +105,134 @@ class (Normed (Matrix t), Linear Vector t, Linear Matrix t) => Field t where
     multiply'    :: Matrix t -> Matrix t -> Matrix t
 
 
--- | Full singular value decomposition using lapack's dgesdd or zgesdd.
+instance Field Double where
+    svd' = svdRd
+    thinSVD' = thinSVDRd
+    sv' = svR
+    luPacked' = luR
+    luSolve' (l_u,perm) = lusR l_u perm
+    linearSolve' = linearSolveR                 -- (luSolve . luPacked) ??
+    linearSolveLS' = linearSolveLSR
+    linearSolveSVD' = linearSolveSVDR Nothing
+    ctrans' = trans
+    eig' = eigR
+    eigSH'' = eigS
+    eigOnly = eigOnlyR
+    eigOnlySH = eigOnlyS
+    cholSH' = cholS
+    qr' = unpackQR . qrR
+    hess' = unpackHess hessR
+    schur' = schurR
+    multiply' = multiplyR
+
+instance Field (Complex Double) where
+    svd' = svdCd
+    thinSVD' = thinSVDCd
+    sv' = svC
+    luPacked' = luC
+    luSolve' (l_u,perm) = lusC l_u perm
+    linearSolve' = linearSolveC
+    linearSolveLS' = linearSolveLSC
+    linearSolveSVD' = linearSolveSVDC Nothing
+    ctrans' = conj . trans
+    eig' = eigC
+    eigOnly = eigOnlyC
+    eigSH'' = eigH
+    eigOnlySH = eigOnlyH
+    cholSH' = cholH
+    qr' = unpackQR . qrC
+    hess' = unpackHess hessC
+    schur' = schurC
+    multiply' = multiplyC
+
+--------------------------------------------------------------
+
+-- | Full singular value decomposition.
 svd :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
 svd = svd'
 
--- | Thin singular value decomposition using lapack's dgesvd or zgesvd.
+-- | A version of 'svd' which returns only the @min (rows m) (cols m)@ singular vectors of @m@.
+--
+-- If @(u,s,v) = thinSVD m@ then @m == u \<> diag s \<> trans v@.
 thinSVD :: Field t => Matrix t -> (Matrix t, Vector Double, Matrix t)
 thinSVD = thinSVD'
 
--- | Singular values using lapack's dgesvd or zgesvd.
-sv :: Field t => Matrix t -> Vector Double
-sv = sv'
+-- | Singular values only.
+singularValues :: Field t => Matrix t -> Vector Double
+singularValues = sv'
+
+-- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
+--
+-- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
+fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
+fullSVD m = (u,d,v) where
+    (u,s,v) = svd m
+    d = diagRect s r c
+    r = rows m
+    c = cols m
+
+-- | Similar to 'thinSVD', returning only the nonzero singular values and the corresponding singular vectors.
+compactSVD :: Field t  => Matrix t -> (Matrix t, Vector Double, Matrix t)
+compactSVD m = (u', subVector 0 d s, v') where
+    (u,s,v) = thinSVD m
+    d = rankSVD (1*eps) m s `max` 1
+    u' = takeColumns d u
+    v' = takeColumns d v
+
+vertical m = rows m >= cols m
+
+-- | Singular values and all right singular vectors.
+rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
+rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
+          | otherwise  = let (_,s,v) = svd m     in (s,v)
+
+-- | Singular values and all right singular vectors.
+leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
+leftSV m  | vertical m = let (u,s,_) = svd m     in (u,s)
+          | otherwise  = let (u,s,_) = thinSVD m in (u,s)
+
+
+{-# DEPRECATED full "use fullSVD instead" #-}
+full svdFun m = (u, d ,v) where
+    (u,s,v) = svdFun m
+    d = diagRect s r c
+    r = rows m
+    c = cols m
+
+{-# DEPRECATED economy "use compactSVD instead" #-}
+economy svdFun m = (u', subVector 0 d s, v') where
+    (u,s,v) = svdFun m
+    d = rankSVD (1*eps) m s `max` 1
+    u' = takeColumns d u
+    v' = takeColumns d v
+
+
+--------------------------------------------------------------
 
 -- | Obtains the LU decomposition of a matrix in a compact data structure suitable for 'luSolve'.
 luPacked    :: Field t => Matrix t -> (Matrix t, [Int])
 luPacked = luPacked'
 
--- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization
---   obtained by 'luPacked'.
+-- | Solution of a linear system (for several right hand sides) from the precomputed LU factorization obtained by 'luPacked'.
 luSolve     :: Field t => (Matrix t, [Int]) -> Matrix t -> Matrix t
 luSolve = luSolve'
 
--- | Solution of a general linear system (for several right-hand sides) using lapacks' dgesv or zgesv.
+-- | Solve a linear system (for square coefficient matrix and several right-hand sides) using the LU decomposition. For underconstrained or overconstrained systems use 'linearSolveLS' or 'linearSolveSVD'.
 -- It is similar to 'luSolve' . 'luPacked', but @linearSolve@ raises an error if called on a singular system.
---  See also other versions of linearSolve in "Numeric.LinearAlgebra.LAPACK".
 linearSolve :: Field t => Matrix t -> Matrix t -> Matrix t
 linearSolve = linearSolve'
+
+-- | Minimum norm solution of a general linear least squares problem Ax=B using the SVD. Admits rank-deficient systems but it is slower than 'linearSolveLS'. The effective rank of A is determined by treating as zero those singular valures which are less than 'eps' times the largest singular value.
 linearSolveSVD :: Field t => Matrix t -> Matrix t -> Matrix t
 linearSolveSVD = linearSolveSVD'
 
+
+-- | Least squared error solution of an overconstrained linear system, or the minimum norm solution of an underconstrained system. For rank-deficient systems use 'linearSolveSVD'.
+linearSolveLS :: Field t => Matrix t -> Matrix t -> Matrix t
+linearSolveLS = linearSolveLS'
+
+--------------------------------------------------------------
+
 -- | Eigenvalues and eigenvectors of a general square matrix.
 --
 -- If @(s,v) = eig m@ then @m \<> v == v \<> diag s@
@@ -150,24 +251,45 @@ eigSH' = eigSH''
 eigenvaluesSH' :: Field t => Matrix t -> Vector Double
 eigenvaluesSH' = eigOnlySH
 
--- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.
-cholSH      :: Field t => Matrix t -> Matrix t
-cholSH = cholSH'
+-- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.
+--
+-- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
+eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
+eigSH m | m `equal` ctrans m = eigSH' m
+        | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
+
+-- | Eigenvalues of a complex hermitian or real symmetric matrix.
+eigenvaluesSH :: Field t => Matrix t -> Vector Double
+eigenvaluesSH m | m `equal` ctrans m = eigenvaluesSH' m
+                | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
+
+--------------------------------------------------------------
 
--- | QR factorization using lapack's dgeqr2 or zgeqr2.
+-- | QR factorization.
 --
 -- If @(q,r) = qr m@ then @m == q \<> r@, where q is unitary and r is upper triangular.
 qr          :: Field t => Matrix t -> (Matrix t, Matrix t)
 qr = qr'
 
--- | Hessenberg factorization using lapack's dgehrd or zgehrd.
+-- | RQ factorization.
+--
+-- If @(r,q) = rq m@ then @m == r \<> q@, where q is unitary and r is upper triangular.
+rq :: Field t => Matrix t -> (Matrix t, Matrix t)
+rq m = (r,q) where
+    (q',r') = qr $ trans $ rev1 m
+    r = rev2 (trans r')
+    q = rev2 (trans q')
+    rev1 = flipud . fliprl
+    rev2 = fliprl . flipud
+
+-- | Hessenberg factorization.
 --
 -- If @(p,h) = hess m@ then @m == p \<> h \<> ctrans p@, where p is unitary
--- and h is in upper Hessenberg form.
+-- and h is in upper Hessenberg form (it has zero entries below the first subdiagonal).
 hess        :: Field t => Matrix t -> (Matrix t, Matrix t)
 hess = hess'
 
--- | Schur factorization using lapack's dgees or zgees.
+-- | Schur factorization.
 --
 -- If @(u,s) = schur m@ then @m == u \<> s \<> ctrans u@, where u is unitary
 -- and s is a Shur matrix. A complex Schur matrix is upper triangular. A real Schur matrix is
@@ -187,137 +309,45 @@ multiply :: Field t => Matrix t -> Matrix t -> Matrix t
 multiply = multiply'
 
 
-instance Field Double where
-    svd' = svdRd
-    thinSVD' = thinSVDRd
-    sv' = svR
-    luPacked' = luR
-    luSolve' (l_u,perm) = lusR l_u perm
-    linearSolve' = linearSolveR                 -- (luSolve . luPacked) ??
-    linearSolveSVD' = linearSolveSVDR Nothing
-    ctrans' = trans
-    eig' = eigR
-    eigSH'' = eigS
-    eigOnly = eigOnlyR
-    eigOnlySH = eigOnlyS
-    cholSH' = cholS
-    qr' = unpackQR . qrR
-    hess' = unpackHess hessR
-    schur' = schurR
-    multiply' = multiplyR
-
-instance Field (Complex Double) where
-    svd' = svdCd
-    thinSVD' = thinSVDCd
-    sv' = svC
-    luPacked' = luC
-    luSolve' (l_u,perm) = lusC l_u perm
-    linearSolve' = linearSolveC
-    linearSolveSVD' = linearSolveSVDC Nothing
-    ctrans' = conj . trans
-    eig' = eigC
-    eigOnly = eigOnlyC
-    eigSH'' = eigH
-    eigOnlySH = eigOnlyH
-    cholSH' = cholH
-    qr' = unpackQR . qrC
-    hess' = unpackHess hessC
-    schur' = schurC
-    multiply' = multiplyC
-
-
--- | Eigenvalues and Eigenvectors of a complex hermitian or real symmetric matrix.
---
--- If @(s,v) = eigSH m@ then @m == v \<> diag s \<> ctrans v@
-eigSH :: Field t => Matrix t -> (Vector Double, Matrix t)
-eigSH m | m `equal` ctrans m = eigSH' m
-        | otherwise = error "eigSH requires complex hermitian or real symmetric matrix"
-
--- | Eigenvalues of a complex hermitian or real symmetric matrix.
-eigenvaluesSH :: Field t => Matrix t -> Vector Double
-eigenvaluesSH m | m `equal` ctrans m = eigenvaluesSH' m
-                | otherwise = error "eigenvaluesSH requires complex hermitian or real symmetric matrix"
+-- | Similar to 'chol' without checking that the input matrix is hermitian or symmetric.
+cholSH      :: Field t => Matrix t -> Matrix t
+cholSH = cholSH'
 
--- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf.
+-- | Cholesky factorization of a positive definite hermitian or symmetric matrix.
 --
 -- If @c = chol m@ then @m == ctrans c \<> c@.
 chol :: Field t => Matrix t ->  Matrix t
 chol m | m `equal` ctrans m = cholSH m
        | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
 
+
+
+
 square m = rows m == cols m
 
--- | determinant of a square matrix, computed from the LU decomposition.
+-- | Determinant of a square matrix.
 det :: Field t => Matrix t -> t
 det m | square m = s * (product $ toList $ takeDiag $ lup)
       | otherwise = error "det of nonsquare matrix"
     where (lup,perm) = luPacked m
           s = signlp (rows m) perm
 
--- | Explicit LU factorization of a general matrix using lapack's dgetrf or zgetrf.
+-- | Explicit LU factorization of a general matrix.
 --
 -- If @(l,u,p,s) = lu m@ then @m == p \<> l \<> u@, where l is lower triangular,
 -- u is upper triangular, p is a permutation matrix and s is the signature of the permutation.
 lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
 lu = luFact . luPacked
 
--- | Inverse of a square matrix using lapacks' dgesv and zgesv.
+-- | Inverse of a square matrix.
 inv :: Field t => Matrix t -> Matrix t
 inv m | square m = m `linearSolve` ident (rows m)
       | otherwise = error "inv of nonsquare matrix"
 
--- | Pseudoinverse of a general matrix using lapack's dgelss or zgelss.
+-- | Pseudoinverse of a general matrix.
 pinv :: Field t => Matrix t -> Matrix t
 pinv m = linearSolveSVD m (ident (rows m))
 
-
-{-# DEPRECATED full "use fullSVD instead" #-}
-full svdFun m = (u, d ,v) where
-    (u,s,v) = svdFun m
-    d = diagRect s r c
-    r = rows m
-    c = cols m
-
--- | A version of 'svd' which returns an appropriate diagonal matrix with the singular values.
---
--- If @(u,d,v) = fullSVD m@ then @m == u \<> d \<> trans v@.
-fullSVD :: Field t => Matrix t -> (Matrix t, Matrix Double, Matrix t)
-fullSVD m = (u,d,v) where
-    (u,s,v) = svd m
-    d = diagRect s r c
-    r = rows m
-    c = cols m
-
-
-{-# DEPRECATED economy "use compactSVD instead" #-}
-economy svdFun m = (u', subVector 0 d s, v') where
-    (u,s,v) = svdFun m
-    d = rankSVD (1*eps) m s `max` 1
-    u' = takeColumns d u
-    v' = takeColumns d v
-
--- | A version of 'svd' which returns only the nonzero singular values and the corresponding rows and columns of the rotations.
---
--- If @(u,s,v) = compactSVD m@ then @m == u \<> diag s \<> trans v@.
-compactSVD :: Field t  => Matrix t -> (Matrix t, Vector Double, Matrix t)
-compactSVD m = (u', subVector 0 d s, v') where
-    (u,s,v) = thinSVD m
-    d = rankSVD (1*eps) m s `max` 1
-    u' = takeColumns d u
-    v' = takeColumns d v
-
-vertical m = rows m >= cols m
-
--- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd.
-rightSV :: Field t => Matrix t -> (Vector Double, Matrix t)
-rightSV m | vertical m = let (_,s,v) = thinSVD m in (s,v)
-          | otherwise  = let (_,s,v) = svd m     in (s,v)
-
--- | Singular values and all right singular vectors, using lapack's dgesvd or zgesvd.
-leftSV :: Field t => Matrix t -> (Matrix t, Vector Double)
-leftSV m  | vertical m = let (u,s,_) = svd m     in (u,s)
-          | otherwise  = let (u,s,_) = thinSVD m in (u,s)
-
 -- | Numeric rank of a matrix from the SVD decomposition.
 rankSVD :: Element t
         => Double   -- ^ numeric zero (e.g. 1*'eps')
@@ -379,12 +409,12 @@ pnormCV PNorm1 = norm1 . mapVector magnitude
 pnormCV Infinity = vectorMax . mapVector magnitude
 --pnormCV _ = error "pnormCV not yet defined"
 
-pnormRM PNorm2 m = sv m @> 0
+pnormRM PNorm2 m = singularValues m @> 0
 pnormRM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (vectorMapR Abs) m
 pnormRM Infinity m = vectorMax $ liftMatrix (vectorMapR Abs) m `mXv` constant 1 (cols m)
 --pnormRM _ _ = error "p norm not yet defined"
 
-pnormCM PNorm2 m = sv m @> 0
+pnormCM PNorm2 m = singularValues m @> 0
 pnormCM PNorm1 m = vectorMax $ constant 1 (rows m) `vXm` liftMatrix (mapVector magnitude) m
 pnormCM Infinity m = vectorMax $ liftMatrix (mapVector magnitude) m `mXv` constant 1 (cols m)
 --pnormCM _ _ = error "p norm not yet defined"
@@ -433,12 +463,12 @@ nullspaceSVD hint a (s,v) = vs where
 
 -- | The nullspace of a matrix. See also 'nullspaceSVD'.
 nullspacePrec :: Field t
-              => Double     -- ^ relative tolerance in 'eps' units
+              => Double     -- ^ relative tolerance in 'eps' units (e.g., use 3 to get 3*'eps')
               -> Matrix t   -- ^ input matrix
               -> [Vector t] -- ^ list of unitary vectors spanning the nullspace
 nullspacePrec t m = nullspaceSVD (Left (t*eps)) m (rightSV m)
 
--- | The nullspace of a matrix, assumed to be one-dimensional, with default tolerance 1*'eps'.
+-- | The nullspace of a matrix, assumed to be one-dimensional, with machine precision.
 nullVector :: Field t => Matrix t -> Vector t
 nullVector = last . nullspacePrec 1
 
@@ -522,11 +552,11 @@ uH (pq, tau) = (p,h)
 -- | Reciprocal of the 2-norm condition number of a matrix, computed from the singular values.
 rcond :: Field t => Matrix t -> Double
 rcond m = last s / head s
-    where s = toList (sv m)
+    where s = toList (singularValues m)
 
 -- | Number of linearly independent rows or columns.
 rank :: Field t => Matrix t -> Int
-rank m = rankSVD eps m (sv m)
+rank m = rankSVD eps m (singularValues m)
 
 {-
 expm' m = case diagonalize (complex m) of
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 7284cd9..8c55486 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -178,6 +178,11 @@ runTests n = do
     putStrLn "------ luSolve"
     test (linearSolveProp (luSolve.luPacked) . rSqWC)
     test (linearSolveProp (luSolve.luPacked) . cSqWC)
+    putStrLn "------ luSolveLS"
+    test (linearSolveProp linearSolveLS . rSqWC)
+    test (linearSolveProp linearSolveLS . cSqWC)
+    test (linearSolveProp2 linearSolveLS . rConsist)
+    test (linearSolveProp2 linearSolveLS . cConsist)
     putStrLn "------ pinv (linearSolveSVD)"
     test (pinvProp  . rM)
     test (pinvProp  . cM)
@@ -220,6 +225,8 @@ runTests n = do
     putStrLn "------ qr"
     test (qrProp     . rM)
     test (qrProp     . cM)
+    test (rqProp     . rM)
+    test (rqProp     . cM)
     putStrLn "------ hess"
     test (hessProp   . rSq)
     test (hessProp   . cSq)
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index 7bb914f..405ce64 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -32,13 +32,13 @@ module Numeric.LinearAlgebra.Tests.Properties (
     svdProp1, svdProp1a, svdProp2, svdProp3, svdProp4,
     svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
     eigProp, eigSHProp, eigProp2, eigSHProp2,
-    qrProp,
+    qrProp, rqProp,
     hessProp,
     schurProp1, schurProp2,
     cholProp,
     expmDiagProp,
     multProp1, multProp2,
-    linearSolveProp
+    linearSolveProp, linearSolveProp2
 ) where
 
 import Numeric.LinearAlgebra
@@ -177,10 +177,11 @@ svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
           (s',v') = rightSVC m
           (u',s'') = leftSVC m
 
-svdProp7 m = s |~| s' && u |~| u' && v |~| v'
+svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
     where (u,s,v) = svd m
           (s',v') = rightSV m
           (u',s'') = leftSV m
+          s''' = singularValues m
 
 ------------------------------------------------------------------
 
@@ -201,6 +202,16 @@ eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m
 qrProp m = q <> r |~| m && unitary q && upperTriang r
     where (q,r) = qr m
 
+rqProp m = r <> q |~| m && unitary q && utr
+    where (r,q) = rq m
+          upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1
+              where t = f-c
+          utr = upptr (rows r) (cols r) * r |~| r
+
+upperTriang' m = rows m == 1 || down |~| z
+    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
 hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
     where (p,h) = hess m
 
@@ -226,3 +237,8 @@ multProp1 (a,b) = a <> b |~| mulH a b
 multProp2 (a,b) = ctrans (a <> b) |~| ctrans b <> ctrans a
 
 linearSolveProp f m = f m m |~| ident (rows m)
+
+linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
+    where q = min (rows a) (cols a)
+          b = a <> x
+          wc = rank a == q