src/Linear/QR.hs

   1 {-# LANGUAGE NoImplicitPrelude #-}
   2 {-# LANGUAGE ScopedTypeVariables #-}
   3
   4 -- | QR factorization via Givens rotations.
   5 --
   6 module Linear.QR (
   7   eigenvalues,
   8   givens_rotator,
   9   qr )
  10 where
  11
  12 import qualified Algebra.Ring as Ring ( C )
  13 import qualified Algebra.Algebraic as Algebraic ( C )
  14 import Data.Vector.Fixed ( N1, S, ifoldl )
  15 import Data.Vector.Fixed.Cont ( Arity )
  16 import NumericPrelude hiding ( (*) )
  17
  18 import Linear.Matrix (
  19   Mat(..),
  20   (*),
  21   (!!!),
  22   construct,
  23   diagonal,
  24   identity_matrix,
  25   transpose )
  26
  27
  28 -- | Construct a givens rotation matrix that will operate on row @i@
  29 --   and column @j@. This is done to create zeros in some column of
  30 --   the target matrix. You must also supply that column's @i@th and
  31 --   @j@th entries as arguments.
  32 --
  33 --   Examples (Watkins, p. 193):
  34 --
  35 --   >>> import Normed ( Normed(..) )
  36 --   >>> import Linear.Vector ( Vec2, Vec3 )
  37 --   >>> import Linear.Matrix ( Mat2, Mat3, fromList, frobenius_norm )
  38 --   >>> import qualified Data.Vector.Fixed as V ( map )
  39 --
  40 --   >>> let m  = givens_rotator 0 1 1 1 :: Mat2 Double
  41 --   >>> let m2 = fromList [[1, -1],[1, 1]] :: Mat2 Double
  42 --   >>> m == (1 / (sqrt 2) :: Double) *> m2
  43 --   True
  44 --
  45 --   >>> let m = fromList [[2,3],[5,7]] :: Mat2 Double
  46 --   >>> let rot =  givens_rotator 0 1 2.0 5.0 :: Mat2 Double
  47 --   >>> ((transpose rot) * m) !!! (1,0) < 1e-12
  48 --   True
  49 --   >>> let (Mat rows) = rot
  50 --   >>> let (Mat cols) = transpose rot
  51 --   >>> V.map norm rows :: Vec2 Double
  52 --   fromList [1.0,1.0]
  53 --   >>> V.map norm cols :: Vec2 Double
  54 --   fromList [1.0,1.0]
  55 --
  56 --   >>> let m = fromList [[12,-51,4],[6,167,-68],[-4,24,-41]] :: Mat3 Double
  57 --   >>> let rot = givens_rotator 1 2 6 (-4) :: Mat3 Double
  58 --   >>> let ex_rot_r1 = [1,0,0] :: [Double]
  59 --   >>> let ex_rot_r2 = [0,0.83205,-0.55470] :: [Double]
  60 --   >>> let ex_rot_r3 = [0, 0.55470, 0.83205] :: [Double]
  61 --   >>> let ex_rot = fromList [ex_rot_r1, ex_rot_r2, ex_rot_r3] :: Mat3 Double
  62 --   >>> frobenius_norm ((transpose rot) - ex_rot) < 1e-4
  63 --   True
  64 --   >>> ((transpose rot) * m) !!! (2,0) == 0
  65 --   True
  66 --   >>> let (Mat rows) = rot
  67 --   >>> let (Mat cols) = transpose rot
  68 --   >>> V.map norm rows :: Vec3 Double
  69 --   fromList [1.0,1.0,1.0]
  70 --   >>> V.map norm cols :: Vec3 Double
  71 --   fromList [1.0,1.0,1.0]
  72 --
  73 givens_rotator :: forall m a. (Eq a, Ring.C a, Algebraic.C a, Arity m)
  74                => Int -> Int -> a -> a -> Mat m m a
  75 givens_rotator i j xi xj =
  76   construct f
  77   where
  78     xnorm = sqrt $ xi^2 + xj^2
  79     c = if xnorm == (fromInteger 0) then (fromInteger 1) else xi / xnorm
  80     s = if xnorm == (fromInteger 0) then (fromInteger 0) else xj / xnorm
  81
  82     f :: Int -> Int -> a
  83     f y z
  84       | y == i && z == i = c
  85       | y == j && z == j = c
  86       | y == i && z == j = negate s
  87       | y == j && z == i = s
  88       | y == z = fromInteger 1
  89       | otherwise = fromInteger 0
  90
  91
  92 -- | Compute the QR factorization of a matrix (using Givens
  93 --   rotations). This is accomplished with two folds: the first one
  94 --   traverses the columns of the matrix from left to right, and the
  95 --   second traverses the entries of the column from top to
  96 --   bottom.
  97 --
  98 --   The state that is passed through the fold is the current (q,r)
  99 --   factorization. We keep the pair updated by multiplying @q@ and
 100 --   @r@ by the new rotator (or its transpose).
 101 --
 102 --   We do not require that the diagonal elements of R are positive,
 103 --   so our factorization is a little less unique than usual.
 104 --
 105 --   Examples:
 106 --
 107 --   >>> import Linear.Matrix
 108 --
 109 --   >>> let m = fromList [[1,2],[1,3]] :: Mat2 Double
 110 --   >>> let (q,r) = qr m
 111 --   >>> let c = (1 / (sqrt 2 :: Double))
 112 --   >>> let ex_q = c *> (fromList [[1,-1],[1,1]] :: Mat2 Double)
 113 --   >>> let ex_r = c *> (fromList [[2,5],[0,1]] :: Mat2 Double)
 114 --   >>> frobenius_norm (q - ex_q) == 0
 115 --   True
 116 --   >>> frobenius_norm (r - ex_r) == 0
 117 --   True
 118 --   >>> let m' = q*r
 119 --   >>> frobenius_norm (m - m') < 1e-10
 120 --   True
 121 --   >>> is_upper_triangular' 1e-10 r
 122 --   True
 123 --
 124 --   >>> let m = fromList [[2,3],[5,7]] :: Mat2 Double
 125 --   >>> let (q,r) = qr m
 126 --   >>> frobenius_norm (m - (q*r)) < 1e-12
 127 --   True
 128 --   >>> is_upper_triangular' 1e-10 r
 129 --   True
 130 --
 131 --   >>> let m = fromList [[12,-51,4],[6,167,-68],[-4,24,-41]] :: Mat3 Double
 132 --   >>> let (q,r) = qr m
 133 --   >>> frobenius_norm (m - (q*r)) < 1e-12
 134 --   True
 135 --   >>> is_upper_triangular' 1e-10 r
 136 --   True
 137 --
 138 qr :: forall m n a. (Arity m, Arity n, Eq a, Algebraic.C a, Ring.C a)
 139    => Mat m n a -> (Mat m m a, Mat m n a)
 140 qr matrix =
 141   ifoldl col_function initial_qr columns
 142   where
 143     Mat columns = transpose matrix
 144     initial_qr = (identity_matrix, matrix)
 145
 146     -- | Process the column and spit out the current QR
 147     --   factorization. In the first column, we want to get rotators
 148     --   Q12, Q13, Q14,... In the second column, we want rotators Q23,
 149     --   Q24, Q25,...
 150     col_function (q,r) col_idx col =
 151       ifoldl (f col_idx) (q,r) col
 152
 153     -- | Process the entries in a column, doing basically the same
 154     --   thing as col_dunction does. It updates the QR factorization,
 155     --   maybe, and returns the current one.
 156     f col_idx (q,r) idx _ -- ignore the current element
 157       | idx <= col_idx = (q,r) -- leave it alone
 158       | otherwise = (q*rotator, (transpose rotator)*r)
 159           where
 160             y = r !!! (idx, col_idx)
 161             rotator :: Mat m m a
 162             rotator = givens_rotator col_idx idx (r !!! (col_idx, col_idx)) y
 163
 164
 165
 166 -- | Determine the eigenvalues of the given @matrix@ using the
 167 --   iterated QR algorithm (see Golub and Van Loan, \"Matrix
 168 --   Computations\").
 169 --
 170 --   Examples:
 171 --
 172 --   >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
 173 --   >>> import Linear.Matrix ( frobenius_norm, fromList, identity_matrix )
 174 --
 175 --   >>> let m = fromList [[1,1],[-2,4]] :: Mat2 Double
 176 --   >>> let actual = eigenvalues 1000 m
 177 --   >>> let expected = fromList [[3],[2]] :: Col2 Double
 178 --   >>> frobenius_norm (actual - expected) < 1e-12
 179 --   True
 180 --
 181 --   >>> let m = identity_matrix :: Mat2 Double
 182 --   >>> let actual = eigenvalues 10 m
 183 --   >>> let expected = fromList [[1],[1]] :: Col2 Double
 184 --   >>> frobenius_norm (actual - expected) < 1e-12
 185 --   True
 186 --
 187 --   >>> let m = fromList [[0,1,0],[0,0,1],[1,-3,3]] :: Mat3 Double
 188 --   >>> let actual = eigenvalues 1000 m
 189 --   >>> let expected = fromList [[1],[1],[1]] :: Col3 Double
 190 --   >>> frobenius_norm (actual - expected) < 1e-2
 191 --   True
 192 --
 193 eigenvalues :: forall m a. (Arity m, Algebraic.C a, Eq a)
 194              => Int
 195              -> Mat (S m) (S m) a
 196              -> Mat (S m) N1 a
 197 eigenvalues iterations matrix =
 198   diagonal (ut_approximation iterations)
 199   where
 200     ut_approximation :: Int -> Mat (S m) (S m) a
 201     ut_approximation 0 = matrix
 202     ut_approximation k = rk*qk where (qk,rk) = qr (ut_approximation (k-1))
 203