From: Michael Orlitzky Date: Sun, 2 Feb 2014 17:12:09 +0000 (-0500) Subject: Add the Linear.QR module and update the ghci/cabal files. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=e01a21693319ed1bf888d543caf817e1c803a88b;p=numerical-analysis.git Add the Linear.QR module and update the ghci/cabal files. --- diff --git a/.ghci b/.ghci index b8019ba..942d73d 100644 --- a/.ghci +++ b/.ghci @@ -8,6 +8,7 @@ src/Integration/Trapezoid.hs src/Linear/Iteration.hs src/Linear/Matrix.hs + src/Linear/QR.hs src/Linear/System.hs src/Linear/Vector.hs src/Misc.hs @@ -22,6 +23,7 @@ import Integration.Simpson import Integration.Trapezoid import Linear.Iteration import Linear.Matrix +import Linear.QR import Linear.System import Linear.Vector import Misc diff --git a/numerical-analysis.cabal b/numerical-analysis.cabal index 2bc7468..fe47ead 100644 --- a/numerical-analysis.cabal +++ b/numerical-analysis.cabal @@ -17,9 +17,19 @@ description: data-files: makefile library - exposed-modules: Integration.Simpson, Integration.Trapezoid, - Linear.Iteration, Linear.Matrix, Linear.System, Linear.Vector, - Misc, Normed, ODE.IVP, Polynomials.Orthogonal, Roots.Simple, + exposed-modules: + Integration.Simpson, + Integration.Trapezoid, + Linear.Iteration, + Linear.Matrix, + Linear.QR, + Linear.System, + Linear.Vector, + Misc, + Normed, + ODE.IVP, + Polynomials.Orthogonal, + Roots.Simple, Roots.Fast build-depends: diff --git a/src/Linear/QR.hs b/src/Linear/QR.hs new file mode 100644 index 0000000..58027bb --- /dev/null +++ b/src/Linear/QR.hs @@ -0,0 +1,92 @@ +{-# LANGUAGE NoImplicitPrelude #-} +{-# LANGUAGE ScopedTypeVariables #-} + +-- | QR factorization via Givens rotations. +-- +module Linear.QR ( + givens_rotator, + qr ) +where + +import qualified Algebra.Ring as Ring ( C ) +import qualified Algebra.Algebraic as Algebraic ( C ) +import Data.Vector.Fixed ( ifoldl ) +import Data.Vector.Fixed.Cont ( Arity ) +import NumericPrelude hiding ( (*) ) + +import Linear.Matrix ( + Mat(..), + (*), + (!!!), + construct, + identity_matrix, + transpose ) + + +-- | Construct a givens rotation matrix that will operate on row @i@ +-- and column @j@. This is done to create zeros in some column of +-- the target matrix. You must also supply that column's @i@th and +-- @j@th entries as arguments. +-- +-- Examples (Watkins, p. 193): +-- +-- >>> import Linear.Matrix ( Mat2, fromList ) +-- >>> let m = givens_rotator 0 1 1 1 :: Mat2 Double +-- >>> let m2 = fromList [[1, -1],[1, 1]] :: Mat2 Double +-- >>> m == (1 / (sqrt 2) :: Double) *> m2 +-- True +-- +givens_rotator :: forall m a. (Ring.C a, Algebraic.C a, Arity m) + => Int -> Int -> a -> a -> Mat m m a +givens_rotator i j xi xj = + construct f + where + xnorm = sqrt $ xi^2 + xj^2 + c = xi / xnorm + s = xj / xnorm + + f :: Int -> Int -> a + f y z + | y == i && z == i = c + | y == j && z == j = c + | y == i && z == j = negate s + | y == j && z == i = s + | y == z = fromInteger 1 + | otherwise = fromInteger 0 + + +-- | Compute the QR factorization of a matrix (using Givens +-- rotations). This is accomplished with two folds: the first one +-- traverses the columns of the matrix from left to right, and the +-- second traverses the entries of the column from top to +-- bottom. +-- +-- The state that is passed through the fold is the current (q,r) +-- factorization. We keep the pair updated by multiplying @q@ and +-- @r@ by the new rotator (or its transpose). +-- +qr :: forall m n a. (Arity m, Arity n, Algebraic.C a, Ring.C a) + => Mat m n a -> (Mat m m a, Mat m n a) +qr matrix = + ifoldl col_function initial_qr columns + where + Mat columns = transpose matrix + initial_qr = (identity_matrix, matrix) + + -- | Process the column and spit out the current QR + -- factorization. In the first column, we want to get rotators + -- Q12, Q13, Q14,... In the second column, we want rotators Q23, + -- Q24, Q25,... + col_function (q,r) col_idx col = + ifoldl (f col_idx) (q,r) col + + -- | Process the entries in a column, doing basically the same + -- thing as col_dunction does. It updates the QR factorization, + -- maybe, and returns the current one. + f col_idx (q,r) idx x + | idx <= col_idx = (q,r) -- leave it alone. + | otherwise = + (q*rotator, (transpose rotator)*r) + where + rotator :: Mat m m a + rotator = givens_rotator col_idx idx (r !!! (idx, col_idx)) x