From: Michael Orlitzky Date: Tue, 4 Feb 2014 22:40:59 +0000 (-0500) Subject: Add comments to Linear.QR about non-convergence. X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=bd9c9a0bc058786a45a7f679e4037a4a814d79d8;p=numerical-analysis.git Add comments to Linear.QR about non-convergence. --- diff --git a/src/Linear/QR.hs b/src/Linear/QR.hs index b628c12..13782c4 100644 --- a/src/Linear/QR.hs +++ b/src/Linear/QR.hs @@ -171,6 +171,9 @@ qr matrix = -- 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 ) @@ -198,12 +201,20 @@ 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 = - diagonal (ut_approximation iterations) - where - ut_approximation :: Int -> Mat (S m) (S m) a - ut_approximation 0 = matrix - ut_approximation k = rk*qk where (qk,rk) = qr (ut_approximation (k-1)) +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 @@ -214,6 +225,9 @@ eigenvalues iterations matrix = -- 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 ) @@ -268,8 +282,7 @@ eigenvectors_symmetric iterations matrix where (t_prev, p_prev) = tp_pair (k-1) (qk,rk) = qr t_prev - tk = rk*qk pk = p_prev*qk - + tk = rk*qk (values, vectors) = (first diagonal) (tp_pair iterations)