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   eigenvectors_symmetric,
   9   givens_rotator,
  10   qr )
  11 where
  12
  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 ( (*) )
  19
  20 import Linear.Matrix (
  21   Col,
  22   Mat(..),
  23   (*),
  24   (!!!),
  25   construct,
  26   diagonal,
  27   identity_matrix,
  28   symmetric,
  29   transpose )
  30
  31
  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.
  36 --
  37 --   Examples (Watkins, p. 193):
  38 --
  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 )
  43 --
  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
  47 --   True
  48 --
  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
  52 --   True
  53 --   >>> let (Mat rows) = rot
  54 --   >>> let (Mat cols) = transpose rot
  55 --   >>> V.map norm rows :: Vec2 Double
  56 --   fromList [1.0,1.0]
  57 --   >>> V.map norm cols :: Vec2 Double
  58 --   fromList [1.0,1.0]
  59 --
  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
  67 --   True
  68 --   >>> ((transpose rot) * m) !!! (2,0) == 0
  69 --   True
  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]
  76 --
  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 =
  80   construct f
  81   where
  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
  85
  86     f :: Int -> Int -> a
  87     f y z
  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
  94
  95
  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
 100 --   bottom.
 101 --
 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).
 105 --
 106 --   We do not require that the diagonal elements of R are positive,
 107 --   so our factorization is a little less unique than usual.
 108 --
 109 --   Examples:
 110 --
 111 --   >>> import Linear.Matrix
 112 --
 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
 119 --   True
 120 --   >>> frobenius_norm (r - ex_r) == 0
 121 --   True
 122 --   >>> let m' = q*r
 123 --   >>> frobenius_norm (m - m') < 1e-10
 124 --   True
 125 --   >>> is_upper_triangular' 1e-10 r
 126 --   True
 127 --
 128 --   >>> let m = fromList [[2,3],[5,7]] :: Mat2 Double
 129 --   >>> let (q,r) = qr m
 130 --   >>> frobenius_norm (m - (q*r)) < 1e-12
 131 --   True
 132 --   >>> is_upper_triangular' 1e-10 r
 133 --   True
 134 --
 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
 138 --   True
 139 --   >>> is_upper_triangular' 1e-10 r
 140 --   True
 141 --
 142 qr :: forall m n a. (Arity m, Arity n, Eq a, Algebraic.C a, Ring.C a)
 143    => Mat m n a -> (Mat m m a, Mat m n a)
 144 qr matrix =
 145   ifoldl col_function initial_qr columns
 146   where
 147     Mat columns = transpose matrix
 148     initial_qr = (identity_matrix, matrix)
 149
 150     -- | Process the column and spit out the current QR
 151     --   factorization. In the first column, we want to get rotators
 152     --   Q12, Q13, Q14,... In the second column, we want rotators Q23,
 153     --   Q24, Q25,...
 154     col_function (q,r) col_idx col =
 155       ifoldl (f col_idx) (q,r) col
 156
 157     -- | Process the entries in a column, doing basically the same
 158     --   thing as col_function does. It updates the QR factorization,
 159     --   maybe, and returns the current one.
 160     f col_idx (q,r) idx _ -- ignore the current element
 161       | idx <= col_idx = (q,r) -- leave it alone
 162       | otherwise = (q*rotator, (transpose rotator)*r)
 163           where
 164             y = r !!! (idx, col_idx)
 165             rotator :: Mat m m a
 166             rotator = givens_rotator col_idx idx (r !!! (col_idx, col_idx)) y
 167
 168
 169
 170 -- | Determine the eigenvalues of the given @matrix@ using the
 171 --   iterated QR algorithm (see Golub and Van Loan, \"Matrix
 172 --   Computations\").
 173 --
 174 --   Warning: this may not converge if there are repeated eigenvalues
 175 --   (in magnitude).
 176 --
 177 --   Examples:
 178 --
 179 --   >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
 180 --   >>> import Linear.Matrix ( frobenius_norm, fromList, identity_matrix )
 181 --
 182 --   >>> let m = fromList [[1,1],[-2,4]] :: Mat2 Double
 183 --   >>> let actual = eigenvalues 1000 m
 184 --   >>> let expected = fromList [[3],[2]] :: Col2 Double
 185 --   >>> frobenius_norm (actual - expected) < 1e-12
 186 --   True
 187 --
 188 --   >>> let m = identity_matrix :: Mat2 Double
 189 --   >>> let actual = eigenvalues 10 m
 190 --   >>> let expected = fromList [[1],[1]] :: Col2 Double
 191 --   >>> frobenius_norm (actual - expected) < 1e-12
 192 --   True
 193 --
 194 --   >>> let m = fromList [[0,1,0],[0,0,1],[1,-3,3]] :: Mat3 Double
 195 --   >>> let actual = eigenvalues 1000 m
 196 --   >>> let expected = fromList [[1],[1],[1]] :: Col3 Double
 197 --   >>> frobenius_norm (actual - expected) < 1e-2
 198 --   True
 199 --
 200 eigenvalues :: forall m a. (Arity m, Algebraic.C a, Eq a)
 201              => Int
 202              -> Mat (S m) (S m) a
 203              -> Col (S m) a
 204 eigenvalues iterations matrix
 205   | iterations <  0 = error "negative iterations requested"
 206   | iterations == 0 = diagonal matrix
 207   | otherwise =
 208       diagonal (ut_approximation (iterations - 1))
 209       where
 210         ut_approximation :: Int -> Mat (S m) (S m) a
 211         ut_approximation 0 = matrix
 212         ut_approximation k = ut_next
 213           where
 214             ut_prev = ut_approximation (k-1)
 215             (qk,rk) = qr ut_prev
 216             ut_next = rk*qk
 217
 218
 219
 220 -- | Compute the eigenvalues and eigenvectors of a symmetric matrix
 221 --   using an iterative QR algorithm. This is similar to what we do in
 222 --   'eigenvalues' except we also return the product of all \"Q\"
 223 --   matrices that we have generated. This turns out to me the matrix
 224 --   of eigenvectors when the original matrix is symmetric. For
 225 --   references see Goluv and Van Loan, \"Matrix Computations\", or
 226 --   \"Calculation of Gauss Quadrature Rules\" by Golub and Welsch.
 227 --
 228 --   Warning: this may not converge if there are repeated eigenvalues
 229 --   (in magnitude).
 230 --
 231 --   Examples:
 232 --
 233 --   >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
 234 --   >>> import Linear.Matrix ( column, frobenius_norm, fromList )
 235 --   >>> import Linear.Matrix ( identity_matrix, vec3d )
 236 --   >>> import Normed ( Normed(..) )
 237 --
 238 --   >>> let m = identity_matrix :: Mat3 Double
 239 --   >>> let (vals, vecs) = eigenvectors_symmetric 100 m
 240 --   >>> let expected_vals = fromList [[1],[1],[1]] :: Col3 Double
 241 --   >>> let expected_vecs = m
 242 --   >>> vals == expected_vals
 243 --   True
 244 --   >>> vecs == expected_vecs
 245 --   True
 246 --
 247 --   >>> let m = fromList [[3,2,4],[2,0,2],[4,2,3]] :: Mat3 Double
 248 --   >>> let (vals, vecs) = eigenvectors_symmetric 1000 m
 249 --   >>> let expected_vals = fromList [[8],[-1],[-1]] :: Col3 Double
 250 --   >>> let v0' = vec3d (2, 1, 2) :: Col3 Double
 251 --   >>> let v0  = (1 / (norm v0') :: Double) *> v0'
 252 --   >>> let v1' = vec3d (-1, 2, 0) :: Col3 Double
 253 --   >>> let v1  = (1 / (norm v1') :: Double) *> v1'
 254 --   >>> let v2' = vec3d (-4, -2, 5) :: Col3 Double
 255 --   >>> let v2  = (1 / (norm v2') :: Double) *> v2'
 256 --   >>> frobenius_norm ((column vecs 0) - v0) < 1e-12
 257 --   True
 258 --   >>> frobenius_norm ((column vecs 1) - v1) < 1e-12
 259 --   True
 260 --   >>> frobenius_norm ((column vecs 2) - v2) < 1e-12
 261 --   True
 262 --
 263 eigenvectors_symmetric :: forall m a. (Arity m, Algebraic.C a, Eq a)
 264                        => Int
 265                        -> Mat (S m) (S m) a
 266                        -> (Col (S m) a, Mat (S m) (S m) a)
 267 eigenvectors_symmetric iterations matrix
 268   | iterations <  0 = error "negative iterations requested"
 269   | iterations == 0 = (diagonal matrix, identity_matrix)
 270   | not $ symmetric matrix = error "argument is not symmetric"
 271   | otherwise =
 272       (values, vectors)
 273       where
 274         -- | We think of \"T\" as an approximation to an
 275         -- upper-triangular matrix from which we get our
 276         -- eigenvalues. The matrix \"P\" is the product of all
 277         -- previous \"Q\"s and its columns approximate the
 278         -- eigenvectors.
 279         tp_pair :: Int -> (Mat (S m) (S m) a, Mat (S m) (S m) a)
 280         tp_pair 0 = (matrix, identity_matrix)
 281         tp_pair k = (tk,pk)
 282                     where
 283                       (t_prev, p_prev) = tp_pair (k-1)
 284                       (qk,rk) = qr t_prev
 285                       pk = p_prev*qk
 286                       tk = rk*qk
 287
 288         (values, vectors) = (first diagonal) (tp_pair iterations)