module Linear.System (
backward_substitute,
- forward_substitute )
+ forward_substitute,
+ solve_positive_definite )
where
-import Data.Vector.Fixed ( Arity, N1 )
+import qualified Algebra.Algebraic as Algebraic ( C )
+import Data.Vector.Fixed ( Arity, S )
import NumericPrelude hiding ( (*), abs )
-import qualified NumericPrelude as NP ( (*) )
import qualified Algebra.Field as Field ( C )
-import Linear.Matrix ( Mat(..), (!!!), construct, transpose )
+import Linear.Matrix (
+ Col,
+ Mat(..),
+ cholesky,
+ diagonal,
+ dot,
+ ifoldl2,
+ ifoldr2,
+ is_lower_triangular,
+ is_upper_triangular,
+ row,
+ set_idx,
+ zip2,
+ transpose )
--- | Solve the system m' * x = b', where m' is upper-triangular. Will
+-- | Solve the system m' * x = b', where m' is lower-triangular. Will
-- probably crash if m' is non-singular. The result is the vector x.
--
-- Examples:
--
--- >>> import Linear.Matrix ( Mat2, Mat3, fromList, vec2d, vec3d )
+-- >>> import Linear.Matrix ( Mat2, Mat3, frobenius_norm, fromList )
+-- >>> import Linear.Matrix ( vec2d, vec3d )
+-- >>> import Naturals ( N7 )
--
-- >>> let identity = fromList [[1,0,0],[0,1,0],[0,0,1]] :: Mat3 Double
-- >>> let b = vec3d (1, 2, 3::Double)
-- >>> forward_substitute m b
-- ((0.5),(0.75))
--
-forward_substitute :: forall a m. (Field.C a, Arity m)
- => Mat m m a
- -> Mat m N1 a
- -> Mat m N1 a
-forward_substitute m' b' = x'
- where
- x' = construct lambda
-
- -- Convenient accessor for the elements of b'.
- b :: Int -> a
- b k = b' !!! (k, 0)
-
- -- Convenient accessor for the elements of m'.
- m :: Int -> Int -> a
- m i j = m' !!! (i, j)
+-- >>> let f1 = [0.0418]
+-- >>> let f2 = [0.0805]
+-- >>> let f3 = [0.1007]
+-- >>> let f4 = [-0.0045]
+-- >>> let f5 = [-0.0332]
+-- >>> let f6 = [-0.0054]
+-- >>> let f7 = [-0.0267]
+-- >>> let big_F = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col N7 Double
+-- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
+-- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
+-- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
+-- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
+-- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
+-- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
+-- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
+-- >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
+-- >>> let r = cholesky big_K
+-- >>> let rt = transpose r
+-- >>> let e1 = [0.0170647785413895] :: [Double]
+-- >>> let e2 = [0.0338] :: [Double]
+-- >>> let e3 = [0.07408] :: [Double]
+-- >>> let e4 = [-0.00183711730708738] :: [Double]
+-- >>> let e5 = [-0.0135538432434003] :: [Double]
+-- >>> let e6 = [-0.00220454076850486] :: [Double]
+-- >>> let e7 = [-0.00689391035624920] :: [Double]
+-- >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Col N7 Double
+-- >>> let actual = forward_substitute rt big_F
+-- >>> frobenius_norm (actual - expected) < 1e-10
+-- True
+--
+forward_substitute :: forall a m. (Eq a, Field.C a, Arity m)
+ => Mat (S m) (S m) a
+ -> Col (S m) a
+ -> Col (S m) a
+forward_substitute matrix b
+ | not (is_lower_triangular matrix) =
+ error "forward substitution on non-lower-triangular matrix"
+ | otherwise = ifoldl2 f zero pairs
+ where
+ -- Pairs (m_ii, b_i) needed at each step.
+ pairs :: Col (S m) (a,a)
+ pairs = zip2 (diagonal matrix) b
- -- Convenient accessor for the elements of x'.
- x :: Int -> a
- x k = x' !!! (k, 0)
+ f :: Int -> Int -> Col (S m) a -> (a, a) -> Col (S m) a
+ f i _ x (mii, bi) = set_idx x (i,0) newval
+ where
+ newval = (bi - (x `dot` (transpose $ row matrix i))) / mii
- -- The second argument to lambda should always be zero here, so we
- -- ignore it.
- lambda :: Int -> Int -> a
- lambda 0 _ = (b 0) / (m 0 0)
- lambda k _ = ((b k) - sum [ (m k j) NP.* (x j) |
- j <- [0..k-1] ]) / (m k k)
--- | Solve the system m*x = b, where m is lower-triangular. Will
--- probably crash if m is non-singular. The result is the vector x.
+-- | Solve the system m*x = b, where m is upper-triangular. Will
+-- probably crash if m is non-singular. The result is the vector
+-- x. The implementation is identical to 'forward_substitute' except
+-- with a right-fold.
--
-- Examples:
--
-- >>> (backward_substitute identity b) == b
-- True
--
-backward_substitute :: (Field.C a, Arity m)
- => Mat m m a
- -> Mat m N1 a
- -> Mat m N1 a
-backward_substitute m =
- forward_substitute (transpose m)
+-- >>> let m1 = fromList [[1,1,1], [0,1,1], [0,0,1]] :: Mat3 Double
+-- >>> let b = vec3d (1,1,1::Double)
+-- >>> backward_substitute m1 b
+-- ((0.0),(0.0),(1.0))
+--
+backward_substitute :: forall m a. (Eq a, Field.C a, Arity m)
+ => Mat (S m) (S m) a
+ -> Col (S m) a
+ -> Col (S m) a
+backward_substitute matrix b
+ | not (is_upper_triangular matrix) =
+ error "backward substitution on non-upper-triangular matrix"
+ | otherwise = ifoldr2 f zero pairs
+ where
+ -- Pairs (m_ii, b_i) needed at each step.
+ pairs :: Col (S m) (a,a)
+ pairs = zip2 (diagonal matrix) b
+
+ f :: Int -> Int -> Col (S m) a -> (a, a) -> Col (S m) a
+ f i _ x (mii, bi) = set_idx x (i,0) newval
+ where
+ newval = (bi - (x `dot` (transpose $ row matrix i))) / mii
-- | Solve the linear system m*x = b where m is positive definite.
-{-
-solve_positive_definite :: Mat v w a -> Mat w z a
+--
+-- Examples:
+--
+-- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
+-- >>> import Naturals ( N7 )
+--
+-- >>> let f1 = [0.0418]
+-- >>> let f2 = [0.0805]
+-- >>> let f3 = [0.1007]
+-- >>> let f4 = [-0.0045]
+-- >>> let f5 = [-0.0332]
+-- >>> let f6 = [-0.0054]
+-- >>> let f7 = [-0.0267]
+-- >>> let big_F = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col N7 Double
+--
+-- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
+-- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
+-- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
+-- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
+-- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
+-- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
+-- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
+-- >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
+--
+-- >>> let e1 = [1871/75000] :: [Double]
+-- >>> let e2 = [899/25000] :: [Double]
+-- >>> let e3 = [463/15625] :: [Double]
+-- >>> let e4 = [-3/4000] :: [Double]
+-- >>> let e5 = [-83/15000] :: [Double]
+-- >>> let e6 = [-9/10000] :: [Double]
+-- >>> let e7 = [-89/50000] :: [Double]
+-- >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Col N7 Double
+-- >>> let actual = solve_positive_definite big_K big_F
+-- >>> frobenius_norm (actual - expected) < 1e-12
+-- True
+--
+solve_positive_definite :: (Arity m, Algebraic.C a, Eq a, Field.C a)
+ => Mat (S m) (S m) a
+ -> Col (S m) a
+ -> Col (S m) a
solve_positive_definite m b = x
where
r = cholesky m
- -- First we solve r^T * y == b for y. Then let y=r*x
- rx = forward_substitute (transpose r) b
- -- Now solve r*x == b.
--}
+ -- Now, r^T*r*x = b. Let r*x = y, so the system looks like
+ -- r^T * y = b. We can solve this for y.
+ y = forward_substitute (transpose r) b
+ -- Now solve r*x = y to find the value of x.
+ x = backward_substitute r y