-- | QR factorization via Givens rotations.
--
module Linear.QR (
+ eigenvalues,
+ eigenvectors_symmetric,
givens_rotator,
qr )
where
import qualified Algebra.Ring as Ring ( C )
import qualified Algebra.Algebraic as Algebraic ( C )
-import Data.Vector.Fixed ( ifoldl )
+import Control.Arrow ( first )
+import Data.Vector.Fixed ( S, ifoldl )
import Data.Vector.Fixed.Cont ( Arity )
-import Debug.Trace
import NumericPrelude hiding ( (*) )
import Linear.Matrix (
+ Col,
Mat(..),
(*),
(!!!),
construct,
+ diagonal,
identity_matrix,
+ symmetric,
transpose )
-- 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
-- >>> 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)
+qr :: forall m n a. (Arity m, Arity n, Eq a, 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
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,
+ -- thing as col_function 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.
+ | idx <= col_idx = (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
+
+
+
+-- | Determine the eigenvalues of the given @matrix@ using the
+-- iterated QR algorithm (see Golub and Van Loan, \"Matrix
+-- Computations\").
+--
+-- Warning: this may not converge if there are repeated eigenvalues
+-- (in magnitude).
+--
+-- Examples:
+--
+-- >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
+-- >>> import Linear.Matrix ( frobenius_norm, fromList, identity_matrix )
+--
+-- >>> let m = fromList [[1,1],[-2,4]] :: Mat2 Double
+-- >>> let actual = eigenvalues 1000 m
+-- >>> let expected = fromList [[3],[2]] :: Col2 Double
+-- >>> frobenius_norm (actual - expected) < 1e-12
+-- True
+--
+-- >>> let m = identity_matrix :: Mat2 Double
+-- >>> let actual = eigenvalues 10 m
+-- >>> let expected = fromList [[1],[1]] :: Col2 Double
+-- >>> frobenius_norm (actual - expected) < 1e-12
+-- True
+--
+-- >>> let m = fromList [[0,1,0],[0,0,1],[1,-3,3]] :: Mat3 Double
+-- >>> let actual = eigenvalues 1000 m
+-- >>> let expected = fromList [[1],[1],[1]] :: Col3 Double
+-- >>> frobenius_norm (actual - expected) < 1e-2
+-- True
+--
+eigenvalues :: forall m a. (Arity m, Algebraic.C a, Eq a)
+ => Int
+ -> Mat (S m) (S m) a
+ -> Col (S m) a
+eigenvalues iterations matrix
+ | iterations < 0 = error "negative iterations requested"
+ | iterations == 0 = diagonal matrix
+ | otherwise =
+ diagonal (ut_approximation (iterations - 1))
+ where
+ ut_approximation :: Int -> Mat (S m) (S m) a
+ ut_approximation 0 = matrix
+ ut_approximation k = ut_next
+ where
+ ut_prev = ut_approximation (k-1)
+ (qk,rk) = qr ut_prev
+ ut_next = rk*qk
+
+
+
+-- | Compute the eigenvalues and eigenvectors of a symmetric matrix
+-- using an iterative QR algorithm. This is similar to what we do in
+-- 'eigenvalues' except we also return the product of all \"Q\"
+-- matrices that we have generated. This turns out to me the matrix
+-- of eigenvectors when the original matrix is symmetric. For
+-- references see Goluv and Van Loan, \"Matrix Computations\", or
+-- \"Calculation of Gauss Quadrature Rules\" by Golub and Welsch.
+--
+-- Warning: this may not converge if there are repeated eigenvalues
+-- (in magnitude).
+--
+-- Examples:
+--
+-- >>> import Linear.Matrix ( Col2, Col3, Mat2, Mat3 )
+-- >>> import Linear.Matrix ( column, frobenius_norm, fromList )
+-- >>> import Linear.Matrix ( identity_matrix, vec3d )
+-- >>> import Normed ( Normed(..) )
+--
+-- >>> let m = identity_matrix :: Mat3 Double
+-- >>> let (vals, vecs) = eigenvectors_symmetric 100 m
+-- >>> let expected_vals = fromList [[1],[1],[1]] :: Col3 Double
+-- >>> let expected_vecs = m
+-- >>> vals == expected_vals
+-- True
+-- >>> vecs == expected_vecs
+-- True
+--
+-- >>> let m = fromList [[3,2,4],[2,0,2],[4,2,3]] :: Mat3 Double
+-- >>> let (vals, vecs) = eigenvectors_symmetric 1000 m
+-- >>> let expected_vals = fromList [[8],[-1],[-1]] :: Col3 Double
+-- >>> let v0' = vec3d (2, 1, 2) :: Col3 Double
+-- >>> let v0 = (1 / (norm v0') :: Double) *> v0'
+-- >>> let v1' = vec3d (-1, 2, 0) :: Col3 Double
+-- >>> let v1 = (1 / (norm v1') :: Double) *> v1'
+-- >>> let v2' = vec3d (-4, -2, 5) :: Col3 Double
+-- >>> let v2 = (1 / (norm v2') :: Double) *> v2'
+-- >>> frobenius_norm ((column vecs 0) - v0) < 1e-12
+-- True
+-- >>> frobenius_norm ((column vecs 1) - v1) < 1e-12
+-- True
+-- >>> frobenius_norm ((column vecs 2) - v2) < 1e-12
+-- True
+--
+eigenvectors_symmetric :: forall m a. (Arity m, Algebraic.C a, Eq a)
+ => Int
+ -> Mat (S m) (S m) a
+ -> (Col (S m) a, Mat (S m) (S m) a)
+eigenvectors_symmetric iterations matrix
+ | iterations < 0 = error "negative iterations requested"
+ | iterations == 0 = (diagonal matrix, identity_matrix)
+ | not $ symmetric matrix = error "argument is not symmetric"
+ | otherwise =
+ (values, vectors)
+ where
+ -- | We think of \"T\" as an approximation to an
+ -- upper-triangular matrix from which we get our
+ -- eigenvalues. The matrix \"P\" is the product of all
+ -- previous \"Q\"s and its columns approximate the
+ -- eigenvectors.
+ tp_pair :: Int -> (Mat (S m) (S m) a, Mat (S m) (S m) a)
+ tp_pair 0 = (matrix, identity_matrix)
+ tp_pair k = (tk,pk)
+ where
+ (t_prev, p_prev) = tp_pair (k-1)
+ (qk,rk) = qr t_prev
+ pk = p_prev*qk
+ tk = rk*qk
+
+ (values, vectors) = (first diagonal) (tp_pair iterations)
projects / numerical-analysis.git / blobdiff