1 {-# LANGUAGE RebindableSyntax #-}
2 {-# LANGUAGE ScopedTypeVariables #-}
3 {-# LANGUAGE TypeFamilies #-}
8 solve_positive_definite )
11 import qualified Algebra.Algebraic as Algebraic ( C )
12 import Data.Vector.Fixed ( Arity, S )
13 import NumericPrelude hiding ( (*), abs )
14 import qualified NumericPrelude as NP ( (*) )
15 import qualified Algebra.Field as Field ( C )
17 import Linear.Matrix (
35 -- | Solve the system m' * x = b', where m' is lower-triangular. Will
36 -- probably crash if m' is non-singular. The result is the vector x.
40 -- >>> import Linear.Matrix ( Mat2, Mat3, frobenius_norm, fromList )
41 -- >>> import Linear.Matrix ( vec2d, vec3d )
42 -- >>> import Naturals ( N7 )
44 -- >>> let identity = fromList [[1,0,0],[0,1,0],[0,0,1]] :: Mat3 Double
45 -- >>> let b = vec3d (1, 2, 3::Double)
46 -- >>> forward_substitute identity b
47 -- ((1.0),(2.0),(3.0))
48 -- >>> (forward_substitute identity b) == b
51 -- >>> let m = fromList [[1,0],[1,1]] :: Mat2 Double
52 -- >>> let b = vec2d (1, 1::Double)
53 -- >>> forward_substitute m b
56 -- >>> let m = fromList [[4,0],[0,2]] :: Mat2 Double
57 -- >>> let b = vec2d (2, 1.5 :: Double)
58 -- >>> forward_substitute m b
61 -- >>> let f1 = [0.0418]
62 -- >>> let f2 = [0.0805]
63 -- >>> let f3 = [0.1007]
64 -- >>> let f4 = [-0.0045]
65 -- >>> let f5 = [-0.0332]
66 -- >>> let f6 = [-0.0054]
67 -- >>> let f7 = [-0.0267]
68 -- >>> let big_F = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col N7 Double
69 -- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
70 -- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
71 -- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
72 -- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
73 -- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
74 -- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
75 -- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
76 -- >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
77 -- >>> let r = cholesky big_K
78 -- >>> let rt = transpose r
79 -- >>> let e1 = [0.0170647785413895] :: [Double]
80 -- >>> let e2 = [0.0338] :: [Double]
81 -- >>> let e3 = [0.07408] :: [Double]
82 -- >>> let e4 = [-0.00183711730708738] :: [Double]
83 -- >>> let e5 = [-0.0135538432434003] :: [Double]
84 -- >>> let e6 = [-0.00220454076850486] :: [Double]
85 -- >>> let e7 = [-0.00689391035624920] :: [Double]
86 -- >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Col N7 Double
87 -- >>> let actual = forward_substitute rt big_F
88 -- >>> frobenius_norm (actual - expected) < 1e-10
91 forward_substitute :: forall a m. (Eq a, Field.C a, Arity m)
95 forward_substitute matrix b
96 | not (is_lower_triangular matrix) =
97 error "forward substitution on non-lower-triangular matrix"
98 | otherwise = ifoldl2 f zero pairs
100 -- Pairs (m_ii, b_i) needed at each step.
101 pairs :: Col (S m) (a,a)
102 pairs = zip2 (diagonal matrix) b
104 f :: Int -> Int -> Col (S m) a -> (a, a) -> Col (S m) a
105 f i _ x (mii, bi) = set_idx x (i,0) newval
107 newval = (bi - (x `dot` (transpose $ row matrix i))) / mii
111 -- | Solve the system m*x = b, where m is upper-triangular. Will
112 -- probably crash if m is non-singular. The result is the vector x.
116 -- >>> import Linear.Matrix ( Mat3, fromList, vec3d )
118 -- >>> let identity = fromList [[1,0,0],[0,1,0],[0,0,1]] :: Mat3 Double
119 -- >>> let b = vec3d (1, 2, 3::Double)
120 -- >>> backward_substitute identity b
121 -- ((1.0),(2.0),(3.0))
122 -- >>> (backward_substitute identity b) == b
125 -- >>> let m1 = fromList [[1,1,1], [0,1,1], [0,0,1]] :: Mat3 Double
126 -- >>> let b = vec3d (1,1,1::Double)
127 -- >>> backward_substitute m1 b
128 -- ((0.0),(0.0),(1.0))
130 backward_substitute :: forall m a. (Eq a, Field.C a, Arity m)
134 backward_substitute m' b'
135 | not (is_upper_triangular m') =
136 error "backward substitution on non-upper-triangular matrix"
139 x' = construct lambda
141 -- Convenient accessor for the elements of b'.
145 -- Convenient accessor for the elements of m'.
147 m i j = m' !!! (i, j)
149 -- Convenient accessor for the elements of x'.
153 -- The second argument to lambda should always be zero here, so we
155 lambda :: Int -> Int -> a
157 | k == n = (b k) / (m k k)
158 | otherwise = ((b k) - sum [ (m k j) NP.* (x j) |
159 j <- [k+1..n] ]) / (m k k)
164 -- | Solve the linear system m*x = b where m is positive definite.
168 -- >>> import Linear.Matrix ( Col4, frobenius_norm, fromList )
169 -- >>> import Naturals ( N7 )
171 -- >>> let f1 = [0.0418]
172 -- >>> let f2 = [0.0805]
173 -- >>> let f3 = [0.1007]
174 -- >>> let f4 = [-0.0045]
175 -- >>> let f5 = [-0.0332]
176 -- >>> let f6 = [-0.0054]
177 -- >>> let f7 = [-0.0267]
178 -- >>> let big_F = fromList [f1,f2,f3,f4,f5,f6,f7] :: Col N7 Double
180 -- >>> let k1 = [6, -3, 0, 0, 0, 0, 0] :: [Double]
181 -- >>> let k2 = [-3, 10.5, -7.5, 0, 0, 0, 0] :: [Double]
182 -- >>> let k3 = [0, -7.5, 12.5, 0, 0, 0, 0] :: [Double]
183 -- >>> let k4 = [0, 0, 0, 6, 0, 0, 0] :: [Double]
184 -- >>> let k5 = [0, 0, 0, 0, 6, 0, 0] :: [Double]
185 -- >>> let k6 = [0, 0, 0, 0, 0, 6, 0] :: [Double]
186 -- >>> let k7 = [0, 0, 0, 0, 0, 0, 15] :: [Double]
187 -- >>> let big_K = fromList [k1,k2,k3,k4,k5,k6,k7] :: Mat N7 N7 Double
189 -- >>> let e1 = [1871/75000] :: [Double]
190 -- >>> let e2 = [899/25000] :: [Double]
191 -- >>> let e3 = [463/15625] :: [Double]
192 -- >>> let e4 = [-3/4000] :: [Double]
193 -- >>> let e5 = [-83/15000] :: [Double]
194 -- >>> let e6 = [-9/10000] :: [Double]
195 -- >>> let e7 = [-89/50000] :: [Double]
196 -- >>> let expected = fromList [e1,e2,e3,e4,e5,e6,e7] :: Col N7 Double
197 -- >>> let actual = solve_positive_definite big_K big_F
198 -- >>> frobenius_norm (actual - expected) < 1e-12
201 solve_positive_definite :: (Arity m, Algebraic.C a, Eq a, Field.C a)
205 solve_positive_definite m b = x
208 -- Now, r^T*r*x = b. Let r*x = y, so the system looks like
209 -- r^T * y = b. We can solve this for y.
210 y = forward_substitute (transpose r) b
211 -- Now solve r*x = y to find the value of x.
212 x = backward_substitute r y