linearSolveLS, rq

1 files changed, 157 insertions, 127 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