src/Linear/QR.hs

   1 {-# LANGUAGE NoImplicitPrelude #-}
   2 {-# LANGUAGE ScopedTypeVariables #-}
   3
   4 -- | QR factorization via Givens rotations.
   5 --
   6 module Linear.QR (
   7   givens_rotator,
   8   qr )
   9 where
  10
  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 ( (*) )
  16
  17 import Linear.Matrix (
  18   Mat(..),
  19   (*),
  20   (!!!),
  21   construct,
  22   identity_matrix,
  23   transpose )
  24
  25
  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.
  30 --
  31 --   Examples (Watkins, p. 193):
  32 --
  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
  37 --   True
  38 --
  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 =
  42   construct f
  43   where
  44     xnorm = sqrt $ xi^2 + xj^2
  45     c = xi / xnorm
  46     s = xj / xnorm
  47
  48     f :: Int -> Int -> a
  49     f y z
  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
  56
  57
  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
  62 --   bottom.
  63 --
  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).
  67 --
  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)
  70 qr matrix =
  71   ifoldl col_function initial_qr columns
  72   where
  73     Mat columns = transpose matrix
  74     initial_qr = (identity_matrix, matrix)
  75
  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,
  79     --   Q24, Q25,...
  80     col_function (q,r) col_idx col =
  81       ifoldl (f col_idx) (q,r) col
  82
  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.
  86     f col_idx (q,r) idx x
  87       | idx <= col_idx = (q,r) -- leave it alone.
  88       | otherwise =
  89           (q*rotator, (transpose rotator)*r)
  90           where
  91             rotator :: Mat m m a
  92             rotator = givens_rotator col_idx idx (r !!! (idx, col_idx)) x