diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -6,32 +6,38 @@ This library provides common numerical analysis functionality, without requiring any external bindings. It is not optimized for performance (yet), but it serves as an experimental platform for scientific computation in a purely functional setting. -Algorithms : +Contents : * Iterative linear solvers - * BiConjugate Gradient (BCG) + * BiConjugate Gradient (`bcg`) - * Conjugate Gradient Squared (CGS) + * Conjugate Gradient Squared (`cgs`) - * BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems) + * BiConjugate Gradient Stabilized (`bicgstab`) (non-Hermitian systems) - * Transpose-Free Quasi-Minimal Residual (TFQMR) + * Transpose-Free Quasi-Minimal Residual (`tfqmr`) +* Direct linear solvers + + * LU-based (`luSolve`) + * Matrix factorization algorithms - * QR + * QR (`qr`) - * LU + * LU (`lu`) - * Cholesky + * Cholesky (`chol`) * Eigenvalue algorithms - * QR algorithm + * Arnoldi iteration (`arnoldi`) - * Rayleigh quotient iteration + * QR (`eigsQR`) + * Rayleigh quotient iteration (`eigRayleigh`) + * Utilities : Vector and matrix norms, matrix condition number, Givens rotation, Householder reflection * Predicates : Matrix orthogonality test (A^T A ~= I) @@ -43,9 +49,9 @@ The module `Numeric.LinearAlgebra.Sparse` contains the user interface. -### Creation and pretty-printing +### Creation of sparse data -To create a sparse matrix from an array of its entries we use `fromListSM` : +The `fromListSM` function creates a sparse matrix from an array of its entries we use : fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a @@ -53,8 +59,14 @@ > amat = fromListSM (3,3) [(0,0,2),(1,0,4),(1,1,3),(1,2,2),(2,2,5)] -And similarly for sparse vectors : `fromListSV :: Int -> [(Int, a)] -> SpVector a`. +and similarly + fromListSV :: Int -> [(Int, a)] -> SpVector a + +can be used to create sparse vectors. + +### Displaying sparse data + Both sparse vectors and matrices can be pretty-printed using `prd`: > prd amat @@ -64,7 +76,7 @@ [4,3,2] [0,0,5] -The zeros are just added at pretty printing time; sparse vectors and matrices should only contain non-zero entries. +The zeros are just added at printing time; sparse vectors and matrices should only contain non-zero entries. ### Matrix operations @@ -79,7 +91,7 @@ [0.0,0.0,5.0] Notice that the result is _dense_, i.e. certain entries are numerically zero but have been inserted into the result along with all the others (thus taking up memory!). -To preserve sparsity, we can use a sparsifying matrix-matrix product `#~#`, which filters out all the elements x for which `|x| <= eps`, where `eps` (defined in `Numeric.Eps`) is fixed at 10^-8. +To preserve sparsity, we can use a sparsifying matrix-matrix product `#~#`, which filters out all the elements x for which `|x| <= eps`, where `eps` (defined in `Numeric.Eps`) depends on the numerical type used (e.g. it is 10^-6 for `Float`s and 10^-12 for `Double`s). > prd $ l #~# u ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 ) @@ -91,7 +103,7 @@ ### Linear systems -Linear systems can be solved with either `linSolve` (which also requires choosing a method) or with `<\>` (which uses BiCGSTAB as default) : +Large sparse linear systems are best solved with iterative methods. `sparse-linear-algebra` provides a selection of these via the `linSolve` function, or alternatively `<\>` (which uses BiCGSTAB as default) : > b = fromListSV 3 [(0,3),(1,2),(2,5)] > x = amat <\> b @@ -106,6 +118,16 @@ ( 3 elements ) , 3 NZ ( sparsity 1.0 ) [2.9999999999999996,1.9999999999999996,4.999999999999999] + +The library also provides a forward-backward substitution solver (`luSolve`) based on a triangular factorization of the system matrix (usually LU). This should be the preferred for solving smaller, dense systems. Using the data defined above we can cross-verify the two solution methods: + + > x' = luSolve l u b + > prd x' + + ( 3 elements ) , 3 NZ ( sparsity 1.0 ) + + [1.5,-2.0,1.0] + diff --git a/sparse-linear-algebra.cabal b/sparse-linear-algebra.cabal --- a/sparse-linear-algebra.cabal +++ b/sparse-linear-algebra.cabal @@ -1,5 +1,5 @@ name: sparse-linear-algebra -version: 0.2.1.0 +version: 0.2.1.1 synopsis: Numerical computation in native Haskell description: Currently the library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details homepage: https://github.com/ocramz/sparse-linear-algebra @@ -34,6 +34,7 @@ , primitive >= 0.6.1.0 , mtl >= 2.2.1 , mwc-random + , vector executable sparse-linear-algebra diff --git a/src/Data/Sparse/Common.hs b/src/Data/Sparse/Common.hs --- a/src/Data/Sparse/Common.hs +++ b/src/Data/Sparse/Common.hs @@ -11,18 +11,20 @@ ( module X, insertRowWith, insertRow, insertColWith, insertCol, diagonalSM, - outerProdSV, (><), toSV, svToSM, + outerProdSV, (><), toSV, svToSM, + lookupRowSM, extractCol, extractRow, extractVectorDenseWith, extractRowDense, extractColDense, extractDiagDense, extractSubRow, extractSubCol, extractSubRow_RK, extractSubCol_RK, matVec, (#>), vecMat, (<#), + fromCols, prd) where import Data.Sparse.Utils as X import Data.Sparse.Types as X -import Data.Sparse.IntMap2.IntMap2 as X +import Data.Sparse.IntMap2.IntMap2 -- as X import Data.Sparse.SpMatrix as X import Data.Sparse.SpVector as X @@ -31,8 +33,20 @@ import qualified Data.IntMap as IM +import Data.Maybe (fromMaybe, maybe) +import qualified Data.Vector as V +-- withBoundsSM m ij e f +-- | isValidIxSM m ij = f m ij +-- | otherwise = error e +-- | modify the size of a SpVector. Do not use directly +resizeSV :: Int -> SpVector a -> SpVector a +resizeSV d2 (SV _ sv) = SV d2 sv + + +mapKeysSV fk (SV d sv) = SV d $ IM.mapKeys fk sv + -- * Insert row/column vector in matrix -- | Insert row , using the provided row index transformation function @@ -107,12 +121,20 @@ | otherwise = error $ "toSV : incompatible matrix dimension " ++ show (m,n) + +-- | Lookup a row in a SpMatrix; returns an SpVector with the row, if this is non-empty +lookupRowSM :: SpMatrix a -> IxRow -> Maybe (SpVector a) +lookupRowSM sm i = SV (ncols sm) <$> IM.lookup i (dat sm) + + -- * Extract a SpVector from an SpMatrix -- ** Sparse extract -- |Extract ith row extractRow :: SpMatrix a -> IxRow -> SpVector a -extractRow m i = toSV $ extractRowSM m i +extractRow m i + | inBounds0 (nrows m) i = fromMaybe (zeroSV (ncols m)) (lookupRowSM m i) + | otherwise = error $ unwords ["extractRow : index",show i,"out of bounds"] -- |Extract jth column extractCol :: SpMatrix a -> IxCol -> SpVector a @@ -143,14 +165,32 @@ --- | extract row interval -extractSubRow :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a -extractSubRow m i (j1, j2) = toSV $ extractSubRowSM m i (j1, j2) +-- | extract row interval (all entries between columns j1 and j2, INCLUDED, are returned) +-- extractSubRow :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a +-- extractSubRow m i (j1, j2) = case lookupRowSM m i of +-- Nothing -> zeroSV (ncols m) +-- Just rv -> ifilterSV (\j _ -> j >= j1 && j <= j2) rv +-- |", returning in Maybe +-- extractSubRow :: SpMatrix a -> IxRow -> (Int, Int) -> Maybe (SpVector a) +-- extractSubRow m i (j1, j2) = +-- resizeSV (j2 - j1) . ifilterSV (\j _ -> j >= j1 && j <= j2) <$> lookupRowSM m i + +-- | Extract an interval of SpVector components, changing accordingly the resulting SpVector size. Keys are _not_ rebalanced, i.e. components are still labeled according with respect to the source matrix. +extractSubRow :: SpMatrix a -> IxRow -> (Int, Int) -> SpVector a +extractSubRow m i (j1, j2) = fromMaybe (zeroSV deltaj) vfilt where + deltaj = j2 - j1 + 1 + vfilt = resizeSV deltaj . + ifilterSV (\j _ -> j >= j1 && j <= j2) <$> lookupRowSM m i + -- | extract row interval, rebalance keys by subtracting lowest one extractSubRow_RK :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a -extractSubRow_RK m i (j1, j2) = toSV $ extractSubRowSM_RK m i (j1, j2) +extractSubRow_RK m i (j1, j2) = mapKeysSV (subtract j1) $ extractSubRow m i (j1, j2) + -- toSV $ extractSubRowSM_RK m i (j1, j2) + + + -- | extract column interval extractSubCol :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpVector a extractSubCol m j (i1, i2) = toSV $ extractSubColSM m j (i1, i2) @@ -190,6 +230,16 @@ (<#) = vecMat + + + + +-- | Pack a V.Vector of SpVectors as columns of an SpMatrix +fromCols :: V.Vector (SpVector a) -> SpMatrix a +fromCols qv = V.ifoldl' ins (zeroSM m n) qv where + n = V.length qv + m = dim $ V.head qv + ins mm i c = insertCol mm c i diff --git a/src/Data/Sparse/IntMap2/IntMap2.hs b/src/Data/Sparse/IntMap2/IntMap2.hs --- a/src/Data/Sparse/IntMap2/IntMap2.hs +++ b/src/Data/Sparse/IntMap2/IntMap2.hs @@ -3,8 +3,9 @@ import Numeric.LinearAlgebra.Class import qualified Data.IntMap.Strict as IM - +import Data.Sparse.Types +import Data.Maybe @@ -60,7 +61,14 @@ lookupIM2 i j imm = IM.lookup i imm >>= IM.lookup j {-# inline lookupIM2 #-} --- |Ppopulate an IM2 from a list of (row index, column index, value) +-- | Lookup with default 0 +lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IxRow, IxCol) -> a +lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j) + + + + +-- |Populate an IM2 from a list of (row index, column index, value) fromListIM2 :: Foldable t => t (IM.Key, IM.Key, a) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) diff --git a/src/Data/Sparse/SpMatrix.hs b/src/Data/Sparse/SpMatrix.hs --- a/src/Data/Sparse/SpMatrix.hs +++ b/src/Data/Sparse/SpMatrix.hs @@ -75,7 +75,12 @@ instance Sparse SpMatrix a where spy = spySM - +instance Num a => SpContainer SpMatrix a where + type ScIx SpMatrix = (Rows, Cols) + scInsert (i,j) = insertSpMatrix i j + scLookup m (i, j) = lookupSM m i j + m @@ d | isValidIxSM m d = m @@! d + | otherwise = error $ "@@ : incompatible indices : matrix size is " ++ show (dim m) ++ ", but user looked up " ++ show d @@ -200,28 +205,29 @@ -- ** Lookup + + + lookupSM :: SpMatrix a -> IxRow -> IxCol -> Maybe a lookupSM (SM _ im) i j = IM.lookup i im >>= IM.lookup j -- | Looks up an element in the matrix with a default (if the element is not found, zero is returned) -lookupWD_SM, (@@!), (@@) :: Num a => SpMatrix a -> (IxRow, IxCol) -> a +lookupWD_SM, (@@!):: Num a => SpMatrix a -> (IxRow, IxCol) -> a lookupWD_SM sm (i,j) = fromMaybe 0 (lookupSM sm i j) -lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IxRow, IxCol) -> a -lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j) + -- | Zero-default lookup, infix form (no bound checking) (@@!) = lookupWD_SM --- | Zero-default lookup, infix form ("safe" : throws exception if lookup is outside matrix bounds) -m @@ d | isValidIxSM m d = m @@! d - | otherwise = error $ "@@ : incompatible indices : matrix size is " ++ show (dim m) ++ ", but user looked up " ++ show d + + -- FIXME : to throw an exception or just ignore the out-of-bound access ? @@ -287,21 +293,11 @@ -- *** Extract i'th row -- | Extract whole row -extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a -extractRowSM sm i = extractSubmatrix sm (i, i) (0, ncols sm - 1) +-- -- moved to Data.Sparse.Common --- | Extract column within a row range -extractSubRowSM :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpMatrix a -extractSubRowSM sm i (j1, j2) = extractSubmatrix sm (i, i) (j1, j2) --- | Extract column within a row range, rebalance keys -extractSubRowSM_RK :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpMatrix a -extractSubRowSM_RK sm i = - extractSubmatrixRebalanceKeys sm (i, i) - - -- *** Extract j'th column -- | Extract whole column extractColSM :: SpMatrix a -> IxCol -> SpMatrix a @@ -346,6 +342,14 @@ d = IM.filterWithKey ff (immSM m) ff irow row = IM.size row == 1 && IM.size (IM.filterWithKey (\j _ -> j == irow) row) == 1 + +-- | Is the matrix lower/upper triangular? +isLowerTriSM, isUpperTriSM :: Eq a => SpMatrix a -> Bool +isLowerTriSM m = m == lm where + lm = ifilterSM (\i j _ -> i >= j) m + +isUpperTriSM m = m == lm where + lm = ifilterSM (\i j _ -> i <= j) m -- |Is the matrix orthogonal? i.e. Q^t ## Q == I isOrthogonalSM :: (Eq a, Epsilon a) => SpMatrix a -> Bool diff --git a/src/Data/Sparse/SpVector.hs b/src/Data/Sparse/SpVector.hs --- a/src/Data/Sparse/SpVector.hs +++ b/src/Data/Sparse/SpVector.hs @@ -22,6 +22,7 @@ import qualified Data.IntMap as IM import qualified Data.Foldable as F +import qualified Data.Vector as V -- * Sparse Vector @@ -74,6 +75,16 @@ spy = spySV + +instance Num a => SpContainer SpVector a where + type ScIx SpVector = Int + scInsert = insertSpVector + scLookup v i = lookupSV i v + v @@ i = lookupDenseSV i v + + + + instance Hilbert SpVector where a `dot` b | dim a == dim b = dot (dat a) (dat b) | otherwise = @@ -133,6 +144,10 @@ + + + + -- | one-hot encoding : `oneHotSV n k` produces a SpVector of length n having 1 at the k-th position oneHotSVU :: Num a => Int -> IxRow -> SpVector a oneHotSVU n k = SV n (IM.singleton k 1) @@ -152,6 +167,28 @@ + + +-- *** Vector-related + +-- | Populate a SpVector with the contents of a Vector. +fromVector :: V.Vector a -> SpVector a +fromVector qv = V.ifoldl' ins (zeroSV n) qv where + n = V.length qv + ins vv i x = insertSpVector i x vv + +-- | Populate a Vector with the entries of a SpVector, discarding the indices (NB: loses sparsity information). +toVector :: SpVector a -> V.Vector a +toVector = V.fromList . snd . unzip . toListSV + +-- | -- | Populate a Vector with the entries of a SpVector, replacing the missing entries with 0 +toVectorDense :: Num a => SpVector a -> V.Vector a +toVectorDense = V.fromList . toDenseListSV + + + + + -- ** Element insertion -- |insert element `x` at index `i` in a preexisting SpVector @@ -220,6 +257,13 @@ -- | Head element headSV :: Num a => SpVector a -> a headSV sv = fromMaybe 0 (IM.lookup 0 (dat sv)) + +-- | Keep the first n components of the SpVector (like `take` for lists) +takeSV, dropSV :: Int -> SpVector a -> SpVector a +takeSV n (SV _ sv) = SV n $ IM.filterWithKey (\i _ -> i < n) sv +-- | Discard the first n components of the SpVector and rebalance the keys (like `drop` for lists) +dropSV n (SV n0 sv) = SV (n0 - n) $ IM.mapKeys (subtract n) $ IM.filterWithKey (\i _ -> i >= n) sv + diff --git a/src/Data/Sparse/Utils.hs b/src/Data/Sparse/Utils.hs --- a/src/Data/Sparse/Utils.hs +++ b/src/Data/Sparse/Utils.hs @@ -1,8 +1,16 @@ module Data.Sparse.Utils where +import qualified Data.Vector as V + -- * Misc. utilities +-- | Wrap a function with a null check, returning in Maybe +harness :: (t -> Bool) -> (t -> a) -> t -> Maybe a +harness q f v | q v = Nothing + | otherwise = Just $ f v + + -- | Componentwise tuple operations -- TODO : use semilattice properties instead maxTup, minTup :: Ord t => (t, t) -> (t, t) -> (t, t) @@ -46,7 +54,7 @@ go _ _ [] = mneutral --- *** Bounds checking +-- ** Bounds checking type LB = Int type UB = Int @@ -63,3 +71,16 @@ inBounds02 :: (UB, UB) -> (Int, Int) -> Bool inBounds02 (bx,by) (i,j) = inBounds0 bx i && inBounds0 by j + + + + +-- ** Safe indexing + + + +head' :: V.Vector a -> Maybe a +head' = harness V.null V.head + +tail' :: V.Vector a -> Maybe (V.Vector a) +tail' = harness V.null V.tail diff --git a/src/Numeric/LinearAlgebra/Class.hs b/src/Numeric/LinearAlgebra/Class.hs --- a/src/Numeric/LinearAlgebra/Class.hs +++ b/src/Numeric/LinearAlgebra/Class.hs @@ -181,9 +181,16 @@ --- class (Set f, Sparse f a) => SparseSet f a --- instance SparseSet SpVector a where + +class Sparse c a => SpContainer c a where + type ScIx c :: * + scInsert :: ScIx c -> a -> c a -> c a + scLookup :: c a -> ScIx c -> Maybe a + -- -- | Lookup with default, infix form ("safe" : should throw an exception if lookup is outside matrix bounds) + (@@) :: c a -> ScIx c -> a + + diff --git a/src/Numeric/LinearAlgebra/Sparse.hs b/src/Numeric/LinearAlgebra/Sparse.hs --- a/src/Numeric/LinearAlgebra/Sparse.hs +++ b/src/Numeric/LinearAlgebra/Sparse.hs @@ -5,19 +5,21 @@ -- * Matrix factorizations qr, lu, chol, - -- * Incomplete LU - ilu0, -- * Condition number conditionNumberSM, -- * Householder reflection hhMat, hhRefl, -- * Givens' rotation givens, + -- * Arnoldi iteration + arnoldi, -- * Eigensolvers eigsQR, eigRayleigh, -- * Linear solvers linSolve, LinSolveMethod, (<\>), - -- ** Methods + -- ** Direct methods + luSolve, + -- ** Iterative methods cgne, tfqmr, bicgstab, cgs, bcg, _xCgne, _xTfq, _xBicgstab, _x, _xBcg, cgsStep, bicgstabStep, @@ -35,7 +37,7 @@ -- * Sparsify data sparsifySV, -- * Iteration combinators - modifyInspectN, runAppendN', + modifyInspectN, runAppendN', untilConverged, diffSqL ) where @@ -65,7 +67,7 @@ -- import qualified Data.List as L import Data.Maybe - +import qualified Data.Vector as V @@ -73,7 +75,7 @@ -- * Sparsify : remove almost-0 elements (|x| < eps) -- | Sparsify an SpVector sparsifySV :: Epsilon a => SpVector a -> SpVector a -sparsifySV (SV d im) = SV d $ IM.filter isNz im +sparsifySV = filterSV isNz @@ -152,8 +154,12 @@ QR.C1 ) To zero out entry A(i, j) we must find row k such that A(k, j) is non-zero but A has zeros in row k for all columns less than j. --} +NB: the current version is quite inefficient in that: +1. the Givens' matrix `G_i` is different from Identity only in 4 entries +2. at each iteration `i` we multiply `G_i` by the previous partial result `M`. Since this corresponds to a rotation, and the `givensCoef` function already computes the value of the resulting non-zero component (output `r`), `G_i ## M` can be simplified by just changing two entries of `M` (i.e. zeroing one out and changing the other into `r`). +-} +{-# inline givens #-} givens :: (Floating a, Epsilon a, Ord a) => SpMatrix a -> IxRow -> IxCol -> SpMatrix a givens mm i j | isValidIxSM mm (i,j) && isSquareSM mm = @@ -165,11 +171,6 @@ a = mm @@ (i', j) b = mm @@ (i, j) -- element to zero out --- |Is the `k`th the first nonzero column in the row? -firstNonZeroColumn :: IM.IntMap a -> IxRow -> Bool -firstNonZeroColumn mm k = isJust (IM.lookup k mm) && - isNothing (IM.lookupLT k mm) - -- |Returns a set of rows {k} that satisfy QR.C1 candidateRows :: IM.IntMap (IM.IntMap a) -> IxRow -> IxCol -> Maybe [IM.Key] candidateRows mm i j | IM.null u = Nothing @@ -177,6 +178,11 @@ u = IM.filterWithKey (\irow row -> irow /= i && firstNonZeroColumn row j) mm +-- |Is the `k`th the first nonzero column in the row? +{-# inline firstNonZeroColumn #-} +firstNonZeroColumn :: IM.IntMap a -> IxRow -> Bool +firstNonZeroColumn mm k = isJust (IM.lookup k mm) && + isNothing (IM.lookupLT k mm) @@ -382,7 +388,7 @@ -- | Apply a function over a range of integer indices, zip the result with it and filter out the almost-zero entries -onRangeSparse :: (Epsilon b, Real b) => (Int -> b) -> [Int] -> [(Int, b)] +onRangeSparse :: Epsilon b => (Int -> b) -> [Int] -> [(Int, b)] onRangeSparse f ixs = filter (isNz . snd) $ zip ixs $ map f ixs @@ -404,14 +410,14 @@ --- Produces the permutation matrix necessary to have a nonzero in position (iref, jref). This is used in the LU factorization -permutAA :: Num b => IxRow -> IxCol -> SpMatrix a -> Maybe (SpMatrix b) -permutAA iref jref (SM (nro,_) mm) - | isJust (lookupIM2 iref jref mm) = Nothing -- eye nro - | otherwise = Just $ permutationSM nro [head u] where - u = IM.keys (ifilterIM2 ff mm) - ff i j _ = i /= iref && - j == jref +-- -- Produces the permutation matrix necessary to have a nonzero in position (iref, jref). This is used in the LU factorization +-- permutAA :: Num b => IxRow -> IxCol -> SpMatrix a -> Maybe (SpMatrix b) +-- permutAA iref jref (SM (nro,_) mm) +-- | isJust (lookupIM2 iref jref mm) = Nothing -- eye nro +-- | otherwise = Just $ permutationSM nro [head u] where +-- u = IM.keys (ifilterIM2 ff mm) +-- ff i j _ = i /= iref && +-- j == jref @@ -433,7 +439,7 @@ (l, u) = lu aa lh = sparsifyLU l aa uh = sparsifyLU u aa - sparsifyLU m m2 = SM (dim m) $ ifilterIM2 f (dat m) where + sparsifyLU m m2 = ifilterSM f m where f i j _ = isJust (lookupSM m2 i j) @@ -441,10 +447,50 @@ +-- * Arnoldi iteration +-- | Given a matrix A and a positive integer `n`, this procedure finds the basis of an order `n` Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q. +-- At the i`th iteration, it finds (i + 1) coefficients (the i`th column of the Hessenberg matrix H) and the (i + 1)`th Krylov vector. +arnoldi :: + (Floating a, Eq a) => SpMatrix a -> Int -> (SpMatrix a, SpMatrix a) +arnoldi aa kn = (fromCols qvfin, hhfin) where + (qvfin, hhfin, _) = execState (modifyUntil tf arnoldiStep) arnInit + tf (_, _, ii) = ii == kn -- termination criterion + (m, n) = dim aa + arnInit = (qv1, hh1, 1) where + q0 = normalize 2 $ onesSV n -- starting basis vector + aq0 = aa #> q0 + h11 = q0 `dot` aq0 + q1nn = (aq0 ^-^ (h11 .* q0)) + hh1 = fromListSM (m + 1, n) [(0, 0, h11), (1, 0, h21)] where + h21 = norm 2 q1nn + q1 = normalize 2 q1nn + qv1 = V.fromList [q0, q1] + arnoldiStep (qv, hh, i) = (qv', hh', i + 1) + where + qi = V.last qv + aqi = aa #> qi + hhcoli = fmap (`dot` aqi) qv -- H_{1, i}, H_{2, i}, .. , H_{m + 1, i} + zv = zeroSV m + qipnn = + aqi ^-^ (V.foldl' (^+^) zv (V.zipWith (.*) hhcoli qv)) -- unnormalized q_{i+1} + qipnorm = singletonSV $ norm 2 qipnn -- normalization factor H_{i+1, i} + qip = normalize 2 qipnn -- q_{i + 1} + hh' = insertCol hh (concatSV (fromVector hhcoli) qipnorm) i -- update H + qv' = V.snoc qv qip -- append q_{i+1} to Krylov basis Q_i + + + + + + + + + + -- * Preconditioning -- | Partition a matrix into strictly subdiagonal, diagonal and strictly superdiagonal parts @@ -472,10 +518,46 @@ --- Linear solver, LU-based +-- * Linear solver, LU-based +-- | Direct solver based on a triangular factorization of the system matrix. +luSolve :: + (Fractional a, Eq a, Epsilon a) => SpMatrix a -> SpMatrix a -> SpVector a -> SpVector a +luSolve ll uu b + | isLowerTriSM ll && isUpperTriSM uu = lubwSolve uu (lufwSolve ll b) + | otherwise = error "luSolve : factors must be triangular matrices" +lufwSolve ll b = sparsifySV v where + (v, _) = execState (modifyUntil q lStep) lInit where + q (_, i) = i == dim b + lStep (ww, i) = (wwi, i + 1) where + lii = ll @@ (i, i) + bi = b @@ i + wi = (bi - r)/lii where + r = extractSubRow ll i (0, i-1) `dot` takeSV i ww + wwi = insertSpVector i wi ww + lInit = (ww0, 1) where + l00 = ll @@ (0, 0) + b0 = b @@ 0 + w0 = b0 / l00 + ww0 = insertSpVector 0 w0 $ zeroSV (dim b) +-- | NB in the computation of `xi` we must rebalance the subrow indices because `dropSV` does that too, in order to take the inner product with consistent index pairs +lubwSolve uu w = sparsifySV x where + (x, _) = execState (modifyUntil q uStep) uInit + q (_, i) = i == (- 1) + uStep (xx, i) = (xxi, i - 1) where + uii = uu @@ (i, i) + wi = w @@ i + xi = (wi - r) / uii where + r = extractSubRow_RK uu i (i + 1, dim w - 1) `dot` dropSV (i + 1) xx + xxi = insertSpVector i xi xx + uInit = (xx0, i - 1) where + i = dim w - 1 + u00 = uu @@ (i, i) + w0 = w @@ i + x0 = w0 / u00 + xx0 = insertSpVector i x0 $ zeroSV (dim w) diff --git a/test/LibSpec.hs b/test/LibSpec.hs --- a/test/LibSpec.hs +++ b/test/LibSpec.hs @@ -20,8 +20,6 @@ import Data.Sparse.Common -import qualified Data.IntMap as IM - import Control.Monad (replicateM) import Control.Monad.State.Strict (execState) @@ -69,12 +67,14 @@ it "permutPairsSM : permutation matrices are orthogonal" $ do let pm0 = permutPairsSM 3 [(0,2), (1,2)] :: SpMatrix Double pm0 ##^ pm0 `shouldBe` eye 3 - pm0 #^# pm0 `shouldBe` eye 3 + pm0 #^# pm0 `shouldBe` eye 3 + it "isLowerTriSM : checks whether matrix is lower triangular" $ + isLowerTriSM tm8' && isUpperTriSM tm8 `shouldBe` True it "modifyInspectN : early termination by iteration count" $ execState (modifyInspectN 2 (nearZero . diffSqL) (/2)) (1 :: Double) `shouldBe` 1/8 it "modifyInspectN : termination by value convergence" $ nearZero (execState (modifyInspectN (2^16) (nearZero . head) (/2)) (1 :: Double)) `shouldBe` True - describe "Numeric.LinearAlgebra.Sparse : Linear solvers" $ do + describe "Numeric.LinearAlgebra.Sparse : Iterative linear solvers" $ do -- it "TFQMR (2 x 2 dense)" $ -- normSq (_xTfq (tfqmr aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True it "BCG (2 x 2 dense)" $ @@ -83,6 +83,9 @@ nearZero (normSq (aa0 <\> b0 ^-^ x0true)) `shouldBe` True it "CGS (2 x 2 dense)" $ nearZero (normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true)) `shouldBe` True + describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers" $ do + it "LU (unoptimized) (4 x 4 sparse)" $ + checkLuSolve aa1 b1 `shouldBe` True describe "Numeric.LinearAlgebra.Sparse : QR decomposition" $ do it "QR (4 x 4 sparse)" $ checkQr tm4 `shouldBe` True @@ -93,11 +96,63 @@ checkLu tm6 `shouldBe` True it "LU (10 x 10 sparse)" $ checkLu tm7 `shouldBe` True - describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices only)" $ do + describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices only)" $ it "chol (5 x 5 sparse)" $ checkChol tm7 `shouldBe` True + describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration" $ do + it "Arnoldi iteration (3 x 3 dense)" $ + checkArnoldi aa2 3 `shouldBe` True + it "Arnoldi iteration (5 x 5 sparse)" $ + checkArnoldi tm7 5 `shouldBe` True + +{- QR-} + + +checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool +checkQr a = c1 && c2 where + (q, r) = qr a + c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a) + c2 = isOrthogonalSM q + + + +{- LU -} + +checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool +checkLu a = lup == a where + (l, u) = lu a + lup = l #~# u + + + +{- Cholesky -} + +checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool +checkChol a = nearZero $ normFrobenius ((l ##^ l) ^-^ a) where + l = chol a + + +{- direct linear solver -} + +checkLuSolve :: (Epsilon a, Real a, Floating a) => SpMatrix a -> SpVector a -> Bool +checkLuSolve amat rhs = nearZero (normSq ( (lmat #> (umat #> xlu)) ^-^ rhs )) + where + (lmat, umat) = lu amat + xlu = luSolve lmat umat rhs + + +{- Arnoldi iteration -} +checkArnoldi :: (Epsilon a, Floating a, Eq a) => SpMatrix a -> Int -> Bool +checkArnoldi aa kn = nearZero $ normFrobenius $ (aa #~# qvprev) ^-^ (qv #~# hh) where + (qv, hh) = arnoldi aa kn + (m, n) = dim qv + qvprev = extractSubmatrix qv (0, m - 1) (0, n - 2) + + + + {- example 0 : 2x2 linear system @@ -112,13 +167,9 @@ -} -aa0 :: SpMatrix Double -aa0 = SM (2,2) im where - im = IM.fromList [(0, aa0r0), (1, aa0r1)] -aa0r0, aa0r1 :: IM.IntMap Double -aa0r0 = IM.fromList [(0,1),(1,2)] -aa0r1 = IM.fromList [(0,3),(1,4)] +aa0 :: SpMatrix Double +aa0 = fromListDenseSM 2 [1,3,2,4] -- b0, x0 : r.h.s and initial solution resp. @@ -161,7 +212,9 @@ b2 = mkSpVectorD 3 [4,-2,4] +aa22 = fromListDenseSM 2 [2,1,1,2] :: SpMatrix Double + -- -- {- @@ -262,41 +315,13 @@ -{- QR-} -checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool -checkQr a = c1 && c2 where - (q, r) = qr a - c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a) - c2 = isOrthogonalSM q -aa22 = fromListDenseSM 2 [2,1,1,2] :: SpMatrix Double - -{- LU -} - -checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool -checkLu a = lup == a where - (l, u) = lu a - lup = l #~# u - - - -{- Cholesky -} - -checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool -checkChol a = nearZero $ normFrobenius ((l ##^ l) ^-^ a) where - l = chol a - - - - - - {- eigenvalues -} @@ -326,8 +351,7 @@ tv0, tv1 :: SpVector Double tv0 = mkSpVectorD 2 [5, 6] - -tv1 = SV 2 $ IM.singleton 0 1 +tv1 = fromListSV 2 [(0,1)] -- wikipedia test matrix for Givens rotation @@ -403,3 +427,10 @@ -- xhatC = head $ runNCGS niter aa b -- -- printDenseSM aa -- return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue)) + + +tm8 :: SpMatrix Double +tm8 = fromListSM (2,2) [(0,0,1), (0,1,1), (1,1,1)] + +tm8' :: SpMatrix Double +tm8' = fromListSM (2,2) [(0,0,1), (1,0,1), (1,1,1)]