{-# 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 Normed ( Normed(..) ) -- >>> import Linear.Vector ( Vec2, Vec3 ) -- >>> import Linear.Matrix ( Mat2, Mat3, fromList, frobenius_norm ) -- >>> import qualified Data.Vector.Fixed as V ( map ) -- -- >>> 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 -- -- >>> let m = fromList [[2,3],[5,7]] :: Mat2 Double -- >>> let rot = givens_rotator 0 1 2.0 5.0 :: Mat2 Double -- >>> ((transpose rot) * m) !!! (1,0) < 1e-12 -- True -- >>> let (Mat rows) = rot -- >>> let (Mat cols) = transpose rot -- >>> V.map norm rows :: Vec2 Double -- fromList [1.0,1.0] -- >>> V.map norm cols :: Vec2 Double -- fromList [1.0,1.0] -- -- >>> let m = fromList [[12,-51,4],[6,167,-68],[-4,24,-41]] :: Mat3 Double -- >>> let rot = givens_rotator 1 2 6 (-4) :: Mat3 Double -- >>> let ex_rot_r1 = [1,0,0] :: [Double] -- >>> let ex_rot_r2 = [0,0.83205,-0.55470] :: [Double] -- >>> let ex_rot_r3 = [0, 0.55470, 0.83205] :: [Double] -- >>> let ex_rot = fromList [ex_rot_r1, ex_rot_r2, ex_rot_r3] :: Mat3 Double -- >>> frobenius_norm ((transpose rot) - ex_rot) < 1e-4 -- True -- >>> ((transpose rot) * m) !!! (2,0) == 0 -- True -- >>> let (Mat rows) = rot -- >>> let (Mat cols) = transpose rot -- >>> V.map norm rows :: Vec3 Double -- fromList [1.0,1.0,1.0] -- >>> V.map norm cols :: Vec3 Double -- fromList [1.0,1.0,1.0] -- givens_rotator :: forall m a. (Eq 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 = if xnorm == (fromInteger 0) then (fromInteger 1) else xi / xnorm s = if xnorm == (fromInteger 0) then (fromInteger 0) else 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). -- -- We do not require that the diagonal elements of R are positive, -- so our factorization is a little less unique than usual. -- -- Examples: -- -- >>> import Linear.Matrix -- -- >>> let m = fromList [[1,2],[1,3]] :: Mat2 Double -- >>> let (q,r) = qr m -- >>> let c = (1 / (sqrt 2 :: Double)) -- >>> let ex_q = c *> (fromList [[1,-1],[1,1]] :: Mat2 Double) -- >>> let ex_r = c *> (fromList [[2,5],[0,1]] :: Mat2 Double) -- >>> frobenius_norm (q - ex_q) == 0 -- True -- >>> frobenius_norm (r - ex_r) == 0 -- True -- >>> let m' = q*r -- >>> frobenius_norm (m - m') < 1e-10 -- True -- >>> is_upper_triangular' 1e-10 r -- True -- -- >>> let m = fromList [[2,3],[5,7]] :: Mat2 Double -- >>> let (q,r) = qr m -- >>> frobenius_norm (m - (q*r)) < 1e-12 -- True -- >>> is_upper_triangular' 1e-10 r -- True -- -- >>> let m = fromList [[12,-51,4],[6,167,-68],[-4,24,-41]] :: Mat3 Double -- >>> let (q,r) = qr m -- >>> frobenius_norm (m - (q*r)) < 1e-12 -- True -- >>> is_upper_triangular' 1e-10 r -- True -- qr :: forall m n a. (Arity m, Arity n, Eq a, Algebraic.C a, Ring.C a, Show 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 _ -- ignore the current element | idx <= col_idx = (q,r) -- trace ("---------------\nidx: " ++ (show idx) ++ ";\ncol_idx: " ++ (show col_idx) ++ "; leaving it alone") (q,r) -- leave it alone. | otherwise = (q*rotator, (transpose rotator)*r) -- trace ("---------------\nidx: " ++ (show idx) ++ ";\ncol_idx: " ++ (show col_idx) ++ ";\nupdating Q and R;\nq: " ++ (show q) ++ ";\nr " ++ (show r) ++ ";\nnew q: " ++ (show $ q*rotator) ++ ";\nnew r: " ++ (show $ (transpose rotator)*r) ++ ";\ny: " ++ (show y) ++ ";\nr[i,j]: " ++ (show (r !!! (col_idx, col_idx)))) -- (q*rotator, (transpose rotator)*r) where y = r !!! (idx, col_idx) rotator :: Mat m m a rotator = givens_rotator col_idx idx (r !!! (col_idx, col_idx)) y