1 {-# LANGUAGE NoImplicitPrelude #-}
2 {-# LANGUAGE ScopedTypeVariables #-}
4 -- | QR factorization via Givens rotations.
8 eigenvectors_symmetric,
13 import qualified Algebra.Ring as Ring ( C )
14 import qualified Algebra.Algebraic as Algebraic ( C )
15 import Control.Arrow ( first )
16 import Data.Vector.Fixed ( S, ifoldl )
17 import Data.Vector.Fixed.Cont ( Arity )
18 import NumericPrelude hiding ( (*) )
20 import Linear.Matrix (
32 -- | Construct a givens rotation matrix that will operate on row @i@
33 -- and column @j@. This is done to create zeros in some column of
34 -- the target matrix. You must also supply that column's @i@th and
35 -- @j@th entries as arguments.
37 -- Examples (Watkins, p. 193):
39 -- >>> import Normed ( Normed(..) )
40 -- >>> import Linear.Vector ( Vec2, Vec3 )
41 -- >>> import Linear.Matrix ( Mat2, Mat3, fromList, frobenius_norm )
42 -- >>> import qualified Data.Vector.Fixed as V ( map )
44 -- >>> let m = givens_rotator 0 1 1 1 :: Mat2 Double
45 -- >>> let m2 = fromList [[1, -1],[1, 1]] :: Mat2 Double
46 -- >>> m == (1 / (sqrt 2) :: Double) *> m2
49 -- >>> let m = fromList [[2,3],[5,7]] :: Mat2 Double
50 -- >>> let rot = givens_rotator 0 1 2.0 5.0 :: Mat2 Double
51 -- >>> ((transpose rot) * m) !!! (1,0) < 1e-12
53 -- >>> let (Mat rows) = rot
54 -- >>> let (Mat cols) = transpose rot
55 -- >>> V.map norm rows :: Vec2 Double
57 -- >>> V.map norm cols :: Vec2 Double
60 -- >>> let m = fromList [[12,-51,4],[6,167,-68],[-4,24,-41]] :: Mat3 Double
61 -- >>> let rot = givens_rotator 1 2 6 (-4) :: Mat3 Double
62 -- >>> let ex_rot_r1 = [1,0,0] :: [Double]
63 -- >>> let ex_rot_r2 = [0,0.83205,-0.55470] :: [Double]
64 -- >>> let ex_rot_r3 = [0, 0.55470, 0.83205] :: [Double]
65 -- >>> let ex_rot = fromList [ex_rot_r1, ex_rot_r2, ex_rot_r3] :: Mat3 Double
66 -- >>> frobenius_norm ((transpose rot) - ex_rot) < 1e-4
68 -- >>> ((transpose rot) * m) !!! (2,0) == 0
70 -- >>> let (Mat rows) = rot
71 -- >>> let (Mat cols) = transpose rot
72 -- >>> V.map norm rows :: Vec3 Double
73 -- fromList [1.0,1.0,1.0]
74 -- >>> V.map norm cols :: Vec3 Double
75 -- fromList [1.0,1.0,1.0]
77 givens_rotator :: forall m a. (Eq a, Ring.C a, Algebraic.C a, Arity m)
78 => Int -> Int -> a -> a -> Mat m m a
79 givens_rotator i j xi xj =
82 xnorm = sqrt $ xi^2 + xj^2
83 c = if xnorm == (fromInteger 0) then (fromInteger 1) else xi / xnorm
84 s = if xnorm == (fromInteger 0) then (fromInteger 0) else xj / xnorm
88 | y == i && z == i = c
89 | y == j && z == j = c
90 | y == i && z == j = negate s
91 | y == j && z == i = s
92 | y == z = fromInteger 1
93 | otherwise = fromInteger 0
96 -- | Compute the QR factorization of a matrix (using Givens
97 -- rotations). This is accomplished with two folds: the first one
98 -- traverses the columns of the matrix from left to right, and the
99 -- second traverses the entries of the column from top to
102 -- The state that is passed through the fold is the current (q,r)
103 -- factorization. We keep the pair updated by multiplying @q@ and
104 -- @r@ by the new rotator (or its transpose).
106 -- We do not require that the diagonal elements of R are positive,
107 -- so our factorization is a little less unique than usual.
111 -- >>> import Linear.Matrix
113 -- >>> let m = fromList [[1,2],[1,3]] :: Mat2 Double
114 -- >>> let (q,r) = qr m
115 -- >>> let c = (1 / (sqrt 2 :: Double))
116 -- >>> let ex_q = c *> (fromList [[1,-1],[1,1]] :: Mat2 Double)
117 -- >>> let ex_r = c *> (fromList [[2,5],[0,1]] :: Mat2 Double)
118 -- >>> frobenius_norm (q - ex_q) == 0
120 -- >>> frobenius_norm (r - ex_r) == 0
123 -- >>> frobenius_norm (m - m') < 1e-10
125 -- >>> is_upper_triangular' 1e-10 r
128 -- >>> let m = fromList [[2,3],[5,7]] :: Mat2 Double
129 -- >>> let (q,r) = qr m
130 -- >>> frobenius_norm (m - (q*r)) < 1e-12
132 -- >>> is_upper_triangular' 1e-10 r
135 -- >>> let m = fromList [[12,-51,4],[6,167,-68],[-4,24,-41]] :: Mat3 Double
136 -- >>> let (q,r) = qr m
137 -- >>> frobenius_norm (m - (q*r)) < 1e-12
139 -- >>> is_upper_triangular' 1e-10 r
142 qr :: forall m n a. (Arity m, Arity n, Eq a, Algebraic.C a, Ring.C a)
144 -> (Mat (S m) (S m) a, Mat (S m) (S n) a)
146 ifoldl col_function initial_qr columns
148 Mat columns = transpose matrix
149 initial_qr = (identity_matrix, matrix)
151 -- | Process the column and spit out the current QR
152 -- factorization. In the first column, we want to get rotators
153 -- Q12, Q13, Q14,... In the second column, we want rotators Q23,
155 col_function (q,r) col_idx col =
156 ifoldl (f col_idx) (q,r) col
158 -- | Process the entries in a column, doing basically the same
159 -- thing as col_function does. It updates the QR factorization,
160 -- maybe, and returns the current one.
161 f col_idx (q,r) idx _ -- ignore the current element
162 | idx <= col_idx = (q,r) -- leave it alone
163 | otherwise = (q*rotator, (transpose rotator)*r)
165 y = r !!! (idx, col_idx)
166 rotator :: Mat (S m) (S m) a
167 rotator = givens_rotator col_idx idx (r !!! (col_idx, col_idx)) y
171 -- | Determine the eigenvalues of the given @matrix@ using the
172 -- iterated QR algorithm (see Golub and Van Loan, \"Matrix
175 -- Warning: this may not converge if there are repeated eigenvalues
180 -- >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
181 -- >>> import Linear.Matrix ( frobenius_norm, fromList, identity_matrix )
183 -- >>> let m = fromList [[1,1],[-2,4]] :: Mat2 Double
184 -- >>> let actual = eigenvalues 1000 m
185 -- >>> let expected = fromList [[3],[2]] :: Col2 Double
186 -- >>> frobenius_norm (actual - expected) < 1e-12
189 -- >>> let m = identity_matrix :: Mat2 Double
190 -- >>> let actual = eigenvalues 10 m
191 -- >>> let expected = fromList [[1],[1]] :: Col2 Double
192 -- >>> frobenius_norm (actual - expected) < 1e-12
195 -- >>> let m = fromList [[0,1,0],[0,0,1],[1,-3,3]] :: Mat3 Double
196 -- >>> let actual = eigenvalues 1000 m
197 -- >>> let expected = fromList [[1],[1],[1]] :: Col3 Double
198 -- >>> frobenius_norm (actual - expected) < 1e-2
201 eigenvalues :: forall m a. (Arity m, Algebraic.C a, Eq a)
205 eigenvalues iterations matrix
206 | iterations < 0 = error "negative iterations requested"
207 | iterations == 0 = diagonal matrix
209 diagonal (ut_approximation (iterations - 1))
211 ut_approximation :: Int -> Mat (S m) (S m) a
212 ut_approximation 0 = matrix
213 ut_approximation k = ut_next
215 ut_prev = ut_approximation (k-1)
221 -- | Compute the eigenvalues and eigenvectors of a symmetric matrix
222 -- using an iterative QR algorithm. This is similar to what we do in
223 -- 'eigenvalues' except we also return the product of all \"Q\"
224 -- matrices that we have generated. This turns out to me the matrix
225 -- of eigenvectors when the original matrix is symmetric. For
226 -- references see Goluv and Van Loan, \"Matrix Computations\", or
227 -- \"Calculation of Gauss Quadrature Rules\" by Golub and Welsch.
229 -- Warning: this may not converge if there are repeated eigenvalues
234 -- >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
235 -- >>> import Linear.Matrix ( column, frobenius_norm, fromList )
236 -- >>> import Linear.Matrix ( identity_matrix, vec3d )
237 -- >>> import Normed ( Normed(..) )
239 -- >>> let m = identity_matrix :: Mat3 Double
240 -- >>> let (vals, vecs) = eigenvectors_symmetric 100 m
241 -- >>> let expected_vals = fromList [[1],[1],[1]] :: Col3 Double
242 -- >>> let expected_vecs = m
243 -- >>> vals == expected_vals
245 -- >>> vecs == expected_vecs
248 -- >>> let m = fromList [[3,2,4],[2,0,2],[4,2,3]] :: Mat3 Double
249 -- >>> let (vals, vecs) = eigenvectors_symmetric 1000 m
250 -- >>> let expected_vals = fromList [[8],[-1],[-1]] :: Col3 Double
251 -- >>> let v0' = vec3d (2, 1, 2) :: Col3 Double
252 -- >>> let v0 = (1 / (norm v0') :: Double) *> v0'
253 -- >>> let v1' = vec3d (-1, 2, 0) :: Col3 Double
254 -- >>> let v1 = (1 / (norm v1') :: Double) *> v1'
255 -- >>> let v2' = vec3d (-4, -2, 5) :: Col3 Double
256 -- >>> let v2 = (1 / (norm v2') :: Double) *> v2'
257 -- >>> frobenius_norm ((column vecs 0) - v0) < 1e-12
259 -- >>> frobenius_norm ((column vecs 1) - v1) < 1e-12
261 -- >>> frobenius_norm ((column vecs 2) - v2) < 1e-12
264 eigenvectors_symmetric :: forall m a. (Arity m, Algebraic.C a, Eq a)
267 -> (Col (S m) a, Mat (S m) (S m) a)
268 eigenvectors_symmetric iterations matrix
269 | iterations < 0 = error "negative iterations requested"
270 | iterations == 0 = (diagonal matrix, identity_matrix)
271 | not $ symmetric matrix = error "argument is not symmetric"
275 -- | We think of \"T\" as an approximation to an
276 -- upper-triangular matrix from which we get our
277 -- eigenvalues. The matrix \"P\" is the product of all
278 -- previous \"Q\"s and its columns approximate the
280 tp_pair :: Int -> (Mat (S m) (S m) a, Mat (S m) (S m) a)
281 tp_pair 0 = (matrix, identity_matrix)
284 (t_prev, p_prev) = tp_pair (k-1)
289 (values, vectors) = (first diagonal) (tp_pair iterations)