1 {-# LANGUAGE NoImplicitPrelude #-}
2 {-# LANGUAGE ScopedTypeVariables #-}
4 -- | QR factorization via Givens rotations.
11 import qualified Algebra.Ring as Ring ( C )
12 import qualified Algebra.Algebraic as Algebraic ( C )
13 import Data.Vector.Fixed ( ifoldl )
14 import Data.Vector.Fixed.Cont ( Arity )
15 import NumericPrelude hiding ( (*) )
17 import Linear.Matrix (
26 -- | Construct a givens rotation matrix that will operate on row @i@
27 -- and column @j@. This is done to create zeros in some column of
28 -- the target matrix. You must also supply that column's @i@th and
29 -- @j@th entries as arguments.
31 -- Examples (Watkins, p. 193):
33 -- >>> import Linear.Matrix ( Mat2, fromList )
34 -- >>> let m = givens_rotator 0 1 1 1 :: Mat2 Double
35 -- >>> let m2 = fromList [[1, -1],[1, 1]] :: Mat2 Double
36 -- >>> m == (1 / (sqrt 2) :: Double) *> m2
39 givens_rotator :: forall m a. (Ring.C a, Algebraic.C a, Arity m)
40 => Int -> Int -> a -> a -> Mat m m a
41 givens_rotator i j xi xj =
44 xnorm = sqrt $ xi^2 + xj^2
50 | y == i && z == i = c
51 | y == j && z == j = c
52 | y == i && z == j = negate s
53 | y == j && z == i = s
54 | y == z = fromInteger 1
55 | otherwise = fromInteger 0
58 -- | Compute the QR factorization of a matrix (using Givens
59 -- rotations). This is accomplished with two folds: the first one
60 -- traverses the columns of the matrix from left to right, and the
61 -- second traverses the entries of the column from top to
64 -- The state that is passed through the fold is the current (q,r)
65 -- factorization. We keep the pair updated by multiplying @q@ and
66 -- @r@ by the new rotator (or its transpose).
68 qr :: forall m n a. (Arity m, Arity n, Algebraic.C a, Ring.C a)
69 => Mat m n a -> (Mat m m a, Mat m n a)
71 ifoldl col_function initial_qr columns
73 Mat columns = transpose matrix
74 initial_qr = (identity_matrix, matrix)
76 -- | Process the column and spit out the current QR
77 -- factorization. In the first column, we want to get rotators
78 -- Q12, Q13, Q14,... In the second column, we want rotators Q23,
80 col_function (q,r) col_idx col =
81 ifoldl (f col_idx) (q,r) col
83 -- | Process the entries in a column, doing basically the same
84 -- thing as col_dunction does. It updates the QR factorization,
85 -- maybe, and returns the current one.
87 | idx <= col_idx = (q,r) -- leave it alone.
89 (q*rotator, (transpose rotator)*r)
92 rotator = givens_rotator col_idx idx (r !!! (idx, col_idx)) x