<html>
<head>
<BASE HREF="http://www.numeric-quest.com/haskell/Orthogonals.html">
<title>
Indexless linear algebra algorithms
</title>
</head>
<body>
<ul>
<center>
<h1>
***
</h1>
<h1>
Indexless linear algebra algorithms
</h1>
<b>
<br>
Orthogonalization, linear equations, eigenvalues and eigenvectors
<br>
Literate Haskell module <i>Orthogonals.lhs</i>
</b>
<p>
Jan Skibinski, <a href="http://www.numeric-quest.com/news/">
Numeric Quest Inc.</a>, Huntsville, Ontario, Canada
<p>
1998.09.19, last modified 1998.12.28
</center>
<hr>
<p>
It has been argued that the functional paradigm offers more
support for scientific computing than the traditional imperative
programming, such as greater similarity of functional implementation
to mathematical specification of a problem. However, efficiency
of scientific algorithms implemented in Haskell is very low compared
to efficiencies of C or Fortran implementations - notwithstanding
the exceptional descriptive power of Haskell.
<dd>
It has been also argued that tradition and inertia are partially
responsible for this sore state and that many functional algorithms
are direct translations of their imperative counterparts.
<dd>
Arrays - with their indexing schemes and destructive updating
are basic tools of imperative programming. But pure functional
languages, which prohibit variable reassignments, cannot compete
with imperative languages by using the same tools and following
similar reasoning and patterns - unless the functional arrays
themselves are designed with performance in mind. This is
a case with Clean, where efficiency of one kind of their arrays
-- strict unboxed array, approaches efficiency of C.
<dd>
But this has not been done for Haskell arrays yet. They are
lazy, boxed and use auxilliary association lists (index, value)
for initialization -- the latter being mostly responsible for
low efficiency of those algorithms that create many interim
arrays.
<dd>
It appears, that -- as long as indexing scheme is not used
for lookups and updates -- Haskell lists are more efficient
than arrays -- at least at the currents state of Haskell.
<p>
With this in mind, we are attempting to demonstrate here
that the indexing traps can be successfully avoided.
This module implements afresh several typical problems from linear
algebra. Standard Haskell lists are employed instead of arrays
and not a single algorithm ever uses indices for lookups
or updates.
<dd>
We do not claim high efficiency of these algorithms; consider
them exploratory. However, we do claim that the clarity of
these algorithms is significantly better than of those functionally
similar algorithms that employ indexing schemes.
<p>
Two major algorithms have been invented and implemented in Haskell:
one for solving systems of linear equations and one for finding
eigenvalues and eigenvectors of almost any type of a square matrix.
This includes symmetric, hermitian, general complex or nonsymmetric
matrices with real eigenvalues.
<dd>
Amazingly, both methods are based on the same factorization, akin
to QR method, but not exactly the same as the standard QR one.
A simple trick allows to extend this method to nonsymmetric real
matrices with complex eigenvalues and thus one method applies to
all types of matrices.
It follows that the eigenvalue/eigenvector problem can be consistently
treated all across the board. In addition, no administrative
(housekeeping) boring trivia is required here and that helps to
clearly explain the mechanisms employed.
</i>
<p>
<hr>
<p>
<b>
Contents
</b>
<p>
<ul>
<p><li>
Notation
<p><li>
Scalar products and vector normalization
<ul>
<li><b>
bra_ket</b>, scalar product
<li><b>
sum_product</b>, a cousin of bra_ket
<li><b>
norm</b>, vector norm
<li><b>
normalized</b>, vector normalized to one
</ul>
<p><li>
Transposition and adjoining of matrices
<ul>
<li><b>
transposed</b>, transposed matrix
<li><b>
adjoint</b>, transposed and conjugated matrix
</ul>
<p><li>
Products involving matrices
<ul>
<li><b>
matrix_matrix</b>, product of two matrices as list of rows
<li><b>
matrix_matrix'</b>, product of two matrices as list of columns
<li><b>
triangle_matrix'</b>, upper triangular matrix times square matrix
<li><b>
matrix_ket</b>, matrix times ket vector
<li><b>
bra_matrix</b>, bra vector times matrix
<li><b>
bra_matrix_ket</b>, matrix multiplied on both sides by vectors
<li><b>
scalar_matrix</b>, scalar times matrix
</ul>
<p><li>
Orthogonalization process
<ul>
<li><b>
orthogonals</b>, set of orthogonal vectors
<li><b>
gram_schmidt</b>, vector perpendicular to a hyperplane
</ul>
<p><li>
Solutions of linear equations by orthogonalization
<ul>
<li><b>
one_ket_triangle</b>, triangularization of one vector equation
<li><b>
one_ket_solution</b>, solution for one unknown vector
<li><b>
many_kets_triangle</b>, triangularization of several vector equations
<li><b>
many_kets_solution</b>, solution for several unknown vectors
</ul>
<p><li>
Matrix inversion
<ul>
<li><b>
inverse</b>, inverse of a matrix
</ul>
<p><li>
QR factorization of matrices provided by "many_kets_triangle"
<ul>
<li><b>
factors_QR</b>, QR alike factorization of matrices
<li><b>
determinant</b>, computation of the determinant based on the QR factorization
</ul>
<p><li>
Similarity transformations and eigenvalues
<ul>
<li><b>
similar_to</b>, matrix obtained by similarity transformation
<li><b>
iterated_eigenvalues</b>, list of approximations of eigenvalues
<li><b>
eigenvalues</b>, final approximation of eigenvalues
</ul>
<p><li>
Preconditioning of real nonsymmetric matrices
<ul>
<li><b>
add_to_diagonal</b>, simple preconditioning method
</ul>
<p><li>
Examples of iterated eigenvalues
<ul>
<li>
Symmetric real matrix
<li>
Hermitian complex matrix
<li>
General complex matrix
<li>
Nonsymmetric real matrix with real eigenvalues
<li>
Nonsymmetric real matrix with complex eigenvalues
</ul>
<p><li>
Eigenvectors for distinct eigenvalues
<ul>
<li><b>
eigenkets</b>, eigenvectors for distinct eigenvalues
</ul>
<p><li>
Eigenvectors for degenerated eigenvalues
<ul>
<li><b>
eigenket'</b>, eigenvector based on a trial vector
</ul>
<p><li>
Auxiliary functions
<ul>
<li><b>
unit_matrix</b>, a unit matrix with 1's on a diagonal
<li><b>
unit_vector</b>, a vector with one non-zero componenet
<li><b>
diagonals</b>, vector made of a matrix diagonal
</ul>
</ul>
<p>
<hr>
<p>
<b>
Notation
</b>
<p>
What follows is written in Dirac's notation, as used
in Quantum Mechanics. Matrices are represented by capital
letters, while vectors come in two varieties:
<ul>
<p><li>
Bra vector x, written < x |, is represented by one-row matrix
<p><li> Ket vector y, written | y >, is represented by one-column matrix
</ul>
<p>
Bra vectors can be obtained from ket vectors by transposition
and conjugation of their components. Conjugation is only
important for complex vectors.
<p>
Scalar product of two vectors | x > and | y > is written
as
<pre>
< x | y >
</pre>
which looks like a bracket and is sometimes called a "bra_ket".
This justifies "bra" and "ket" names introduced by Dirac. There
is a good reason for conjugating the components of "bra-vector":
the scalar product of
<pre>
< x | x >
</pre>
should be a square of the norm of the vector "x", and that
means that it should be represented by a real number, or complex
number but with its imaginary part equal to zero.
<p>
<hr>
<p>
<pre>
> module Orthogonals where
> import Complex
> import Ratio
> import qualified List
</pre>
<b>
Scalar product and vector normalization
</b>
<p>
The scalar product "bra_ket" is a basis of many algorithms
presented here.
<pre>
> bra_ket :: (Scalar a, Num a) => [a] -> [a] -> a
> bra_ket u v =
> --
> -- Scalar product of two vectors u and v,
> -- or < u | v > in Dirac's notation.
> -- This is equally valid for both: real and complex vectors.
> --
> sum_product u (map coupled v)
</pre>
Notice the call to function "coupled" in the above implementation
of scalar product. This function conjugates its argument
if it is complex, otherwise does not change it. It is defined
in the class Scalar - specifically designed for this purpose
mainly.
<dd>
This class also defines a norm of a vector that might be used
by some algorithms. So far we have been able to avoid this.
<pre>
> class Scalar a where
> coupled :: a->a
> norm :: [a] -> a
> almostZero :: a -> Bool
> scaled :: [a] -> [a]
> instance Scalar Double where
> coupled x = x
> norm u = sqrt (bra_ket u u)
> almostZero x = (abs x) < 1.0e-8
> scaled = scaled'
> instance Scalar Float where
> coupled x = x
> norm u = sqrt (bra_ket u u)
> almostZero x = (abs x) < 1.0e-8
> scaled = scaled'
> instance (Integral a) => Scalar (Ratio a) where
> coupled x = x
> -- norm u = fromDouble ((sqrt (bra_ket u u))::Double)
> -- Intended hack to silently convert to and from Double.
> -- But I do not know how to declare it properly.
> --
> -- Our type Fraction, when used instead of Ratio a, has its own
> -- definition of sqrt. No hack would be needed here.
> almostZero x = abs x < 1e-8
> scaled = scaled'
> instance (RealFloat a) => Scalar (Complex a) where
> coupled (x:+y) = x:+(-y)
> norm u = sqrt (realPart (bra_ket u u)) :+ 0
> almostZero z = (realPart (abs z)) < 1.0e-8
> scaled u = [(x/m):+(y/m) | x:+y <- u]
> where m = maximum [max (abs x) (abs y) | x:+y <- u]
> norm1 :: (Num a) => [a] -> a
> norm1 = sum . map abs
> norminf :: (Num a, Ord a) => [a] -> a
> norminf = maximum . map abs
> matnorm1 :: (Num a, Ord a) => [[a]] -> a
> matnorm1 = matnorminf . transposed
> matnorminf :: (Num a, Ord a) => [[a]] -> a
> matnorminf = maximum . map norm1
</pre>
But we also need a slightly different definition of
scalar product that will appear in multiplication of matrices
by vectors (or vice versa): a straightforward accumulated product
of two lists, where no complex conjugation takes place.
We will call it a 'sum_product".
<pre>
> sum_product :: Num a => [a] -> [a] -> a
> sum_product u v =
> --
> -- Similar to scalar product but without
> -- conjugations of | u > components
> -- Used in matrix-vector or vector-matrix products
> --
> sum (zipWith (*) u v)
</pre>
Some algorithms might need vectors normalized to one, although
we'll try to avoid the normalizations due to its high cost
or its inapplicability to rational numbers. Instead, we wiil
scale vectors by their maximal components.
<pre>
> normalized :: (Scalar a, Fractional a) => [a] -> [a]
> normalized u =
> [uk/n | uk <- u]
> where
> n = norm u
> scaled' u =
> [uk/um | uk <- u]
> where
> um = maximum [abs uk| uk <- u]
</pre>
<hr>
<p>
<b>
Transposition and adjoining of matrices
</b>
<p>
Matrices are represented here by lists of lists.
Function "transposed" converts from row-wise to column-wise
representation, or vice versa.
<dd>
When transposition is combined with complex conjugation
the resulting matrix is called "adjoint".
<p>
A square matrix is called symmetric if it is equal to its transpose
<pre>
A = A<sup>T</sup>
</pre>
It is called Hermitian, or self-adjoint, if it equals to
its adjoint
<pre>
A = A<sup>+</sup>
> transposed :: [[a]] -> [[a]]
> transposed a
> | null (head a) = []
> | otherwise = ([head mi| mi <- a])
> :transposed ([tail mi| mi <- a])
> adjoint :: Scalar a => [[a]] -> [[a]]
> adjoint a
> | null (head a) = []
> | otherwise = ([coupled (head mi)| mi <- a])
> :adjoint ([tail mi| mi <- a])
</pre>
<p>
<hr>
<p>
<b>
Linear combination and sum of two matrices
</b>
<p>
One can form a linear combination of two matrices, such
as
<pre>
C = alpha A + beta B
where
alpha and beta are scalars
</pre>
The most generic form of any combination, not necessary
linear, of components of two matrices is given by "matrix_zipWith"
function below, which accepts a function "f" describing such
combination. For the linear combination with two scalars
the function "f" could be defined as:
<pre>
f alpha beta a b = alpha*a + beta*b
</pre>
For a straightforward addition of two matrices this auxiliary
function is simply "(+)".
<pre>
> matrix_zipWith f a b =
> --
> -- Matrix made of a combination
> -- of matrices a and b - as specified by f
> --
> [zipWith f ak bk | (ak,bk) <- zip a b]
> add_matrices a b = matrix_zipWith (+)
</pre>
<p>
<hr>
<p>
<b>
Products involving matrices
</b>
<p>
Variety of products involving matrices can be defined.
Our Haskell implementation is based on lists of lists
and therefore is open to interpretation: sublists
can either represent the rows or the columns of a matrix.
<dd>
The following definitions are somehow arbitrary, since
one can choose alternative interpretations of lists
representing matrices.
<p>
<b>
C = A B
</b>
<p>
Inner product of two matrices A B can be expressed quite simply,
providing that matrix A is represented by a list of rows
and B - by a list of columns. Function "matrix_matrix"
answers list of rows, while "matrix_matrix'" - list
of columns.
<dd>
Major algorithms of this module make use of "triangle_matrix'",
which calculates a product of upper triangular matrix
with square matrix and returns a rectangular list of columns.
<pre>
> matrix_matrix :: Num a => [[a]] -> [[a]] -> [[a]]
> matrix_matrix a b
> --
> -- A matrix being an inner product
> -- of matrices A and B, where
> -- A is represented by a list of rows a
> -- B is represented by a list of columns b
> -- result is represented by list of rows
> -- Require: length of a is equal of length of b
> -- Require: all sublists are of equal length
>
> | null a = []
> | otherwise = ([sum_product (head a) bi | bi <- b])
> : matrix_matrix (tail a) b
> matrix_matrix' a b
> --
> -- Similar to "matrix_matrix"
> -- but the result is represented by
> -- a list of columns
> --
> | null b = []
> | otherwise = ([sum_product ai (head b) | ai <- a])
> : matrix_matrix' a (tail b)
> triangle_matrix' :: Num a => [[a]] -> [[a]] -> [[a]]
> triangle_matrix' r q =
> --
> -- List of columns of of a product of
> -- upper triangular matrix R and square
> -- matrix Q
> -- where
> -- r is a list of rows of R
> -- q is a list of columns of A
> --
> [f r qk | qk <- q]
> where
> f t u
> | null t = []
> | otherwise = (sum_product (head t) u)
> : (f (tail t) (tail u))
</pre>
<b>
| u > = A | v >
</b>
<p>
Product of a matrix and a ket-vector is another
ket vector. The following implementation assumes
that list "a" represents rows of matrix A.
<pre>
> matrix_ket :: Num a => [[a]] -> [a] -> [a]
> matrix_ket a v = [sum_product ai v| ai <- a]
</pre>
<b>
< u | = < v | A
</b>
<p>
Bra-vector multiplied by a matrix produces
another bra-vector. The implementation below
assumes that list "a" represents columns
of matrix A. It is also assumed that vector
"v" is given in its standard "ket" representation,
therefore the definition below uses "bra_ket"
instead of "sum_product".
<pre>
> bra_matrix :: (Scalar a, Num a) => [a] -> [[a]] -> [a]
> bra_matrix v a = [bra_ket v ai | ai <- a]
</pre>
<b>
alpha = < u | A | v >
</b>
<p>
This kind of product results in a scalar and is often
used to define elements of a new matrix, such as
<pre>
B[i,j] = < ei | A | ej >
</pre>
The implementation below assumes that list "a" represents
rows of matrix A.
<pre>
> bra_matrix_ket :: (Scalar a, Num a) => [a] -> [[a]] -> [a] -> a
> bra_matrix_ket u a v =
> bra_ket u (matrix_ket a v)
</pre>
<b>
B = alpha A
</b>
<p>
Below is a function which multiplies matrix by a scalar:
<pre>
> scalar_matrix :: Num a => a -> [[a]] -> [[a]]
> scalar_matrix alpha a =
> [[alpha*aij| aij <- ai] | ai<-a]
</pre>
<p>
<hr>
<p>
<b>
Orthogonalization process
</b>
<p>
Gram-Schmidt orthogonalization procedure is used here
for calculation of sets of mutually orthogonal vectors.
<dd>
Function "orthogonals" computes a set of mutually orthogonal
vectors - all orthogonal to a given vector. Such set plus
the input vector form a basis of the vector space. Another
words, they are the base vectors, although we cannot call them
unit vectors since we do not normalize them for two reasons:
<ul>
<li>
None of the algorithms presented here needs this -- quite
costly -- normalization.
<li>
Some algorithms can be used either with doubles or with
rationals. The neat output of the latter is sometimes desirable
for pedagogical or accuracy reasons. But normalization requires "sqrt"
function, which is not defined for rational numbers. We could
use our module Fraction instead, where "sqrt" is defined,
but we'll leave it for a future revision of this module.
</ul>
<p>
Function "gram_schmidt" computes one vector - orthogonal
to an incomplete set of orthogonal vectors, which form a hyperplane
in the vector space. Another words, "gram_schmidt" vector is
perpendicular to such a hyperlane.
<pre>
> orthogonals :: (Scalar a, Fractional a) => [a] -> [[a]]
> orthogonals x =
> --
> -- List of (n-1) linearly independent vectors,
> -- (mutually orthogonal) and orthogonal to the
> -- vector x, but not normalized,
> -- where
> -- n is a length of x.
> --
> orth [x] size (next (-1))
> where
> orth a n m
> | n == 1 = drop 1 (reverse a)
> | otherwise = orth ((gram_schmidt a u ):a) (n-1) (next m)
> where
> u = unit_vector m size
> size = length x
> next i = if (i+1) == k then (i+2) else (i+1)
> k = length (takeWhile (== 0) x) -- first non-zero component of x
> gram_schmidt :: (Scalar a, Fractional a) => [[a]] -> [a] -> [a]
> gram_schmidt a u =
> --
> -- Projection of vector | u > on some direction
> -- orthogonal to the hyperplane spanned by the list 'a'
> -- of mutually orthogonal (linearly independent)
> -- vectors.
> --
> gram_schmidt' a u u
> where
> gram_schmidt' a u v
> | null a = v
> | all (== 0) e = gram_schmidt' (tail a) u v
> | otherwise = gram_schmidt' (tail a) u v'
> where
> v' = vectorCombination v (-(bra_ket e u)/(bra_ket e e)) e
> e = head a
> vectorCombination x c y
> | null x = []
> | null y = []
> | otherwise = (head x + c * (head y))
> : (vectorCombination (tail x) c (tail y))
</pre>
<p>
<hr>
<p>
<b>
Solutions of linear equations by orthogonalization
</b>
<p>
A matrix equation for unknown vector | x >
<pre>
A | x > = | b >
</pre>
can be rewritten as
<pre>
x1 | 1 > + x2 | 2 > + x3 | 3 > + ... + xn | n > = | b > (7.1)
where
| 1 >, | 2 >... represent columns of the matrix A
</pre>
For any n-dimensional vector, such as "1", there exist
n-1 linearly independent vectors "ck" that are orthogonal to "1";
that is, each satisfies the relation:
<pre>
< ck | 1 > = 0, for k = 1...m, where m = n - 1
</pre>
If we could find all such vectors, then we could multiply
the equation (7.1) by each of them, and end up with m = n-1
following equations
<pre>
< c1 | 2 > x2 + < c1 | 3 > x3 + ... < c1 | n > xn = < c1 | b >
< c2 | 2 > x2 + < c2 | 3 > x3 + ... < c2 | n > xn = < c2 | b >
.......
< cm | 2 > x2 + < cm | 3 > x3 + ... < cm | n > xn = < cm | b >
</pre>
But the above is nothing more than a new matrix equation
<pre>
A' | x' > = | b' >
or
x2 | 2'> + x3 | 3'> .... + xn | n'> = | b'>
where
primed vectors | 2' >, etc. are the columns of the new
matrix A'.
</pre>
with the problem dimension reduced by one.
<dd>
Taking as an example a four-dimensional problem and writing
down the successive transformations of the original equation
we will end up with the following triangular pattern made of
four vector equations:
<pre>
x1 | 1 > + x2 | 2 > + x3 | 3 > + x4 | 4 > = | b >
x2 | 2'> + x3 | 3'> + x4 | 4'> = | b'>
x3 | 3''> + x4 | 4''> = | b''>
x4 | 4'''> = | b'''>
</pre>
But if we premultiply each vector equation by a non-zero vector
of our choice, say < 1 | , < 2' |, < 3'' |, and < 4''' | - chosen
correspondingly for equations 1, 2, 3 and 4, then the above
system of vector equations will be converted to much simpler
system of scalar equations. The result is
shown below in matrix representation:
<pre>
| p11 p12 p13 p14 | | x1 | = | q1 |
| 0 p22 p23 p24 | | x2 | = | q2 |
| 0 0 p33 p34 | | x3 | = | q3 |
| 0 0 0 p44 | | x4 | = | q4 |
</pre>
In effect, we have triangularized our original matrix A.
Below is a function that does that for any problem size:
<pre>
> one_ket_triangle :: (Scalar a, Fractional a) => [[a]] -> [a] -> [([a],a)]
> one_ket_triangle a b
> --
> -- List of pairs: (p, q) representing
> -- rows of triangular matrix P and of vector | q >
> -- in the equation P | x > = | q >, which
> -- has been obtained by linear transformation
> -- of the original equation A | x > = | b >
> --
> | null a = []
> | otherwise = (p,q):(one_ket_triangle a' b')
> where
> p = [bra_ket u ak | ak <- a]
> q = bra_ket u b
> a' = [[bra_ket ck ai | ck <- orth] | ai <- v]
> b' = [ bra_ket ck b | ck <- orth]
> orth = orthogonals u
> u = head a
> v = tail a
</pre>
The triangular system of equations can be easily solved by
successive substitutions - starting with the last equation.
<pre>
> one_ket_solution :: (Fractional a, Scalar a) => [[a]] -> [a] -> [a]
> one_ket_solution a b =
> --
> -- List representing vector |x>, which is
> -- a solution of the matrix equation
> -- A |x> = |b>
> -- where
> -- a is a list of columns of matrix A
> -- b is a list representing vector |b>
> --
> solve' (unzip (reverse (one_ket_triangle a b))) []
> where
> solve' (a, b) xs
> | null a = xs
> | otherwise = solve' ((tail a), (tail b)) (x:xs)
> where
> x = (head b - (sum_product (tail u) xs))/(head u)
> u = head a
</pre>
The triangularization procedure can be easily extended
to a list of several ket-vectors | b > on the right hand
side of the original equation A | x > = | b > -- instead
of just one:
<pre>
> many_kets_triangle :: (Scalar a, Fractional a) => [[a]] -> [[a]] -> [([a],[a])]
> many_kets_triangle a b
> --
> -- List of pairs: (p, q) representing
> -- rows of triangular matrix P and of rectangular matrix Q
> -- in the equation P X = Q, which
> -- has been obtained by linear transformation
> -- of the original equation A X = B
> -- where
> -- a is a list of columns of matrix A
> -- b is a list of columns of matrix B
> --
> | null a = []
> | otherwise = (p,q):(many_kets_triangle a' b')
> where
> p = [bra_ket u ak | ak <- a]
> q = [bra_ket u bk | bk <- b]
> a' = [[bra_ket ck ai | ck <- orth] | ai <- v]
> b' = [[bra_ket ck bi | ck <- orth] | bi <- b]
> orth = orthogonals u
> u = head a
> v = tail a
</pre>
Similarly, function 'one_ket_solution' can be generalized
to function 'many_kets_solution' that handles cases with
several ket-vectors on the right hand side.
<pre>
> many_kets_solution a b =
> --
> -- List of columns of matrix X, which is
> -- a solution of the matrix equation
> -- A X = B
> -- where
> -- a is a list of columns of matrix A
> -- b is a list of columns of matrix B
> --
> solve' p q emptyLists
> where
> (p, q) = unzip (reverse (many_kets_triangle a b))
> emptyLists = [[] | k <- [1..(length (head q))]]
> solve' a' b' x
> | null a' = x
> | otherwise = solve' (tail a') (tail b')
> [(f vk xk):xk | (xk, vk) <- (zip x v)]
> where
> f vk xk = (vk - (sum_product (tail u) xk))/(head u)
> u = head a'
> v = head b'
</pre>
<p>
<hr>
<p>
<b>
Matrix inversion
</b>
<p>
Function 'many_kets_solution' can be used to compute
inverse of matrix A by specializing matrix B to a unit
matrix I:
<pre>
A X = I
</pre>
It follows that matrix X is an inverse of A; that is X = A<sup>-1</sup>.
<pre>
> inverse :: (Scalar a, Fractional a) => [[a]] -> [[a]]
> inverse a = many_kets_solution a (unit_matrix (length a))
> --
> -- List of columns of inverse of matrix A
> -- where
> -- a is list of columns of A
</pre>
<p>
<hr>
<p>
<b>
QR factorization of matrices
</b>
<p>
The process described above and implemented by
'many_kets_triangle' function transforms the equation
<pre>
A X = B
</pre>
into another equation for the same matrix X
<pre>
R X = S
</pre>
where R is an upper triangular matrix. All operations
performed on matrices A and B during this process are linear,
and therefore we should be able to find a square matrix Q
that describes the entire process in one step. Indeed, assuming
that matrix A can be decomposed as a product of unknown matrix Q
and triangular matrix R and that Q<sup>-1</sup> is an inverse of matrix Q
we can reach the last equation by following these steps:
<pre>
A X = B
(Q R) X = B
Q<sup>-1</sup> Q R X = Q<sup>-1</sup> B
R X = S
</pre>
It follows that during this process a given matrix B
transforms to matrix S, as delivered by 'many_kets_triangle':
<pre>
S = Q<sup>-1</sup> B
</pre>
from which the inverse of Q can be found:
<pre>
Q<sup>-1</sup> = S B<sup>-1</sup>
</pre>
Having a freedom of choice of the right hand side matrix B
we can choose the unit matrix I in place of B, and therefore
simplify the definition of Q<sup>-1</sup>:
<pre>
Q<sup>-1</sup> = S, if B is unit matrix
</pre>
It follows that any non-singular matrix A can be decomposed
as a product of a matrix Q and a triangular matrix R
<pre>
A = Q R
</pre>
where matrices Q<sup>-1</sup> and R are delivered by "many_kets_triangle"
as a result of triangularization process of equation:
<pre>
A X = I
</pre>
The function below extracts a pair of matrices Q and R
from the answer provided by "many_kets_triangle".
During this process it inverts matrix Q<sup>-1</sup> to Q.
This factorization will be used by a sequence of similarity
transformations to be defined in the next section.
<pre>
> factors_QR :: (Fractional a, Scalar a) => [[a]] -> ([[a]],[[a]])
> factors_QR a =
> --
> -- A pair of matrices (Q, R), such that
> -- A = Q R
> -- where
> -- R is upper triangular matrix in row representation
> -- (without redundant zeros)
> -- Q is a transformation matrix in column representation
> -- A is square matrix given as columns
> --
> (inverse (transposed q1),r)
> where
> (r, q1) = unzip (many_kets_triangle a (unit_matrix (length a)))
</pre>
<p>
<hr>
<p>
<b>
Computation of the determinant
</b>
<!-- added by Henning Thielemann -->
<pre>
> determinant :: (Fractional a, Scalar a) => [[a]] -> a
> determinant a =
> let (q,r) = factors_QR a
> -- matrix Q is not normed so we have to respect the norms of its rows
> in product (map norm q) * product (map head r)
</pre>
Naive division-free computation of the determinant by expanding the first column.
It consumes n! multiplications.
<pre>
> determinantNaive :: (Num a) => [[a]] -> a
> determinantNaive [] = 1
> determinantNaive m =
> sum (alternate
> (zipWith (*) (map head m)
> (map determinantNaive (removeEach (map tail m)))))
</pre>
Compute the determinant with about n^4 multiplications
without division according to the clow decomposition algorithm
of Mahajan and Vinay, and Berkowitz
as presented by Günter Rote:
<a href="http://page.inf.fu-berlin.de/~rote/Papers/pdf/Division-free+algorithms.pdf">
Division-Free Algorithms for the Determinant and the Pfaffian:
Algebraic and Combinatorial Approaches</a>.
<pre>
> determinantClow :: (Num a) => [[a]] -> a
> determinantClow [] = 1
> determinantClow m =
> let lm = length m
> in parityFlip lm (last (newClow m
> (nest (lm-1) (longerClow m)
> (take lm (iterate (0:) [1])))))
</pre>
Compute the weights of all clow sequences
where the last clow is closed and a new one is started.
<pre>
> newClow :: (Num a) => [[a]] -> [[a]] -> [a]
> newClow a c =
> scanl (-) 0
> (sumVec (zipWith (zipWith (*)) (List.transpose a) c))
</pre>
Compute the weights of all clow sequences
where the last (open) clow is extended by a new arc.
<pre>
> extendClow :: (Num a) => [[a]] -> [[a]] -> [[a]]
> extendClow a c =
> map (\ai -> sumVec (zipWith scaleVec ai c)) a
</pre>
Given the matrix of all weights of clows of length l
compute the weight matrix for all clows of length (l+1).
Take the result of 'newClow' as diagonal
and the result of 'extendClow' as lower triangle
of the weight matrix.
<pre>
> longerClow :: (Num a) => [[a]] -> [[a]] -> [[a]]
> longerClow a c =
> let diagonal = newClow a c
> triangle = extendClow a c
> in zipWith3 (\i t d -> take i t ++ [d]) [0 ..] triangle diagonal
</pre>
Auxiliary functions for the clow determinant.
<pre>
> {- | Compositional power of a function,
> i.e. apply the function n times to a value. -}
> nest :: Int -> (a -> a) -> a -> a
> nest 0 _ x = x
> nest n f x = f (nest (n-1) f x)
>
> {- successively select elements from xs and remove one in each result list -}
> removeEach :: [a] -> [[a]]
> removeEach xs =
> zipWith (++) (List.inits xs) (tail (List.tails xs))
>
> alternate :: (Num a) => [a] -> [a]
> alternate = zipWith id (cycle [id, negate])
>
> parityFlip :: Num a => Int -> a -> a
> parityFlip n x = if even n then x else -x
>
> {-| Weight a list of numbers by a scalar. -}
> scaleVec :: (Num a) => a -> [a] -> [a]
> scaleVec k = map (k*)
>
> {-| Add corresponding numbers of two lists. -}
> {- don't use zipWith because it clips to the shorter list -}
> addVec :: (Num a) => [a] -> [a] -> [a]
> addVec x [] = x
> addVec [] y = y
> addVec (x:xs) (y:ys) = x+y : addVec xs ys
>
> {-| Add some lists. -}
> sumVec :: (Num a) => [[a]] -> [a]
> sumVec = foldl addVec []
</pre>
<p>
<hr>
<p>
<b>
Similarity transformations and eigenvalues
</b>
<p>
Two n-square matrices A and B are called similar if there
exists a non-singular matrix S such that:
<pre>
B = S<sup>-1</sup> A S
</pre>
It can be proven that:
<ul>
<li>
Any two similar matrices have the same eigenvalues
<li>
Every n-square matrix A is similar to a triangular matrix
whose diagonal elements are the eigenvalues of A.
</ul>
<p>
If matrix A can be transformed to a triangular or a diagonal
matrix Ak by a sequence of similarity transformations then
the eigenvalues of matrix A are the diagonal elements of Ak.
<p>
Let's construct the sequence of matrices similar to A
<pre>
A, A1, A2, A3...
</pre>
by the following iterations - each of which factorizes a matrix
by applying the function 'factors_QR' and then forms a product
of the factors taken in the reverse order:
<pre>
A = Q R = Q (R Q) Q<sup>-1</sup> = Q A1 Q<sup>-1</sup> =
= Q (Q1 R1) Q<sup>-1</sup> = Q Q1 (R1 Q1) Q1<sup>-1</sup> Q<sup>-1</sup> = Q Q1 A2 Q1<sup>-1</sup> Q<sup>-1</sup> =
= Q Q1 (Q2 R2) Q1<sup>-1</sup> Q<sup>-1</sup> = ...
</pre>
We are hoping that after some number of iterations some matrix
Ak would become triangular and therefore its diagonal
elements could serve as eigenvalues of matrix A. As long as
a matrix has real eigenvalues only, this method should work well.
This applies to symmetric and hermitian matrices. It appears
that general complex matrices -- hermitian or not -- can also
be handled this way. Even more, this method also works for some
nonsymmetric real matrices, which have real eigenvalues only.
<dd>
The only type of matrices that cannot be treated by this algorithm
are real nonsymmetric matrices, whose some eigenvalues are complex.
There is no operation in the process that converts real elements
to complex ones, which could find their way into diagonal
positions of a triangular matrix. But a simple preconditioning
of a matrix -- described in the next section -- replaces
a real matrix by a complex one, whose eigenvalues are related
to the eigenvalues of the matrix being replaced. And this allows
us to apply the same method all across the board.
<dd>
It is worth noting that a process known in literature as QR
factorization is not uniquely defined and different algorithms
are employed for this. The algorithms using QR factorization
apply only to symmetric or hermitian matrices, and Q matrix
must be either orthogonal or unitary.
<dd>
But our transformation matrix Q is not orthogonal nor unitary,
although its first row is orthogonal to all other rows. In fact,
this factorization is only similar to QR factorization. We just
keep the same name to help identify a category of the methods
to which it belongs.
<dd>
The same factorization is used for tackling two major problems:
solving the systems of linear equations and finding the eigenvalues
of matrices.
<dd>
Below is the function 'similar_to', which makes a new matrix that is
similar to a given matrix by applying our similarity transformation.
<dd>
Function 'iterated_eigenvalues' applies this transformation n
times - storing diagonals of each new matrix as approximations of
eigenvalues.
<dd>
Function 'eigenvalues' follows the same process but reports the last
approximation only.
<pre>
> similar_to :: (Fractional a, Scalar a) => [[a]] -> [[a]]
> similar_to a =
> --
> -- List of columns of matrix A1 similar to A
> -- obtained by factoring A as Q R and then
> -- forming the product A1 = R Q = (inverse Q) A Q
> -- where
> -- a is list of columns of A
> --
> triangle_matrix' r q
> where
> (q,r) = factors_QR a
> iterated_eigenvalues a n
> --
> -- List of vectors representing
> -- successive approximations of
> -- eigenvalues of matrix A
> -- where
> -- a is a list of columns of A
> -- n is a number of requested iterations
> --
> | n == 0 = []
> | otherwise = (diagonals a)
> : iterated_eigenvalues (similar_to a) (n-1)
> eigenvalues a n
> --
> -- Eigenvalues of matrix A
> -- obtained by n similarity iterations
> -- where
> -- a are the columns of A
> --
> | n == 0 = diagonals a
> | otherwise = eigenvalues (similar_to a) (n-1)
</pre>
<p>
<hr>
<p>
<b>
Preconditioning of real nonsymmetric matrices
</b>
<p>
As mentioned above, our QR-like factorization method works
well with almost all kind of matrices, but with the exception
of a class of real nonsymmetric matrices that have
complex eigenvalues.
<dd>
There is no mechanism in that method that would be able to
produce complex eigenvalues out of the real components of
this type of nonsymmetric matrices. Simple trivial replacement
of real components of a matrix by its complex counterparts
does not work because zero-valued imaginary components do
not contribute in any way to production of nontrivial
imaginary components during the factorization process.
<dd>
What we need is a trick that replaces real nonsymmetric matrix
by a nontrivial complex matrix in such a way that the results
of such replacements could be undone when the series of
similarity transformations finally produced the expected
effect in a form of a triangular matrix.
<dd>
The practical solution is surprisingly simple:
it's suffice to add any complex number, such as "i", to the
main diagonal of a matrix, and when triangularization is done
-- subtract it back from computed eigenvalues.
The explanation follows.
<p>
Consider the eigenproblem for real and nonsymmetric matrix A.
<pre>
A | x > = a | x >
</pre>
Let us now define a new complex matrix B, such that:
<pre>
B = A + alpha I
where
I is a unit matrix and alpha is a complex scalar
</pre>
It is obvious that matrices A and B commute; that is:
<pre>
A B = B A
</pre>
It can be proven that if two matrices commute then they
have the same eigenvectors. Therefore we can use vector
| x > of matrix A as an eigenvector of B:
<pre>
B | x > = b | x >
B | x > = A | x > + alpha I | x >
= a | x > + alpha | x >
= (a + alpha) | x >
</pre>
It follows that eigenvalues of B are related to the eigenvalues
of A by:
<pre>
b = a + alpha
</pre>
After eigenvalues of complex matrix B have been succesfully
computed, all what remains is to subtract "alpha" from them
all to obtain eigenvalues of A. And nothing has to be done
to eigenvectors of B - they are the same for A as well.
Simple and elegant!
<p>
Below is an auxiliary function that adds a scalar to the
diagonal of a matrix:
<pre>
> add_to_diagonal :: Num a => a -> [[a]] -> [[a]]
> add_to_diagonal alpha a =
> --
> -- Add constant alpha to diagonal of matrix A
> --
> [f ai ni | (ai,ni) <- zip a [0..(length a -1)]]
> where
> f b k = p++[head q + alpha]++(tail q)
> where
> (p,q) = splitAt k b
>
</pre>
<p>
<hr>
<p>
<b>
Examples of iterated eigenvalues
</b>
<p>
Here is an example of a symmetric real matrix with results
of application of function 'iterated_eigenvalues'.
<pre>
| 7 -2 1 |
|-2 10 -2 |
| 1 -2 7 |
[[7.0, 10.0, 7.0],
[8.66667, 9.05752, 6.27582],
[10.7928, 7.11006, 6.09718],
[11.5513, 6.40499, 6.04367],
[11.7889, 6.18968, 6.02142],
[11.8943, 6.09506, 6.01068],
[11.9468, 6.04788, 6.00534],
[11.9733, 6.02405, 6.00267],
[11.9866, 6.01206, 6.00134],
[11.9933, 6.00604, 6.00067],
[11.9966, 6.00302, 6.00034],
[11.9983, 6.00151, 6.00017],
[11.9992, 6.00076, 6.00008],
[11.9996, 6.00038, 6.00004],
[11.9998, 6.00019, 6.00002],
[11.9999, 6.00010, 6.00001],
[11.9999, 6.00005, 6.00001]]
The true eigenvalues are:
12, 6, 6
</pre>
Here is an example of a hermitian matrix. (Eigenvalues of hermitian
matrices are real.) The algorithm works well and converges fast.
<pre>
| 2 0 i|
[ 0 1 0 |
[ -i 0 2 |
[[2.8 :+ 0.0, 1.0 :+ 0.0, 1.2 :+ 0.0],
[2.93979 :+ 0.0, 1.0 :+ 0.0, 1.06021 :+ 0.0],
[2.97972 :+ 0.0, 1.0 :+ 0.0, 1.02028 :+ 0.0],
[2.9932 :+ 0.0, 1.0 :+ 0.0, 1.0068 :+ 0.0],
[2.99773 :+ 0.0, 1.0 :+ 0.0, 1.00227 :+ 0.0],
[2.99924 :+ 0.0, 1.0 :+ 0.0, 1.00076 :+ 0.0],
[2.99975 :+ 0.0, 1.0 :+ 0.0, 1.00025 :+ 0.0],
[2.99992 :+ 0.0, 1.0 :+ 0.0, 1.00008 :+ 0.0],
[2.99997 :+ 0.0, 1.0 :+ 0.0, 1.00003 :+ 0.0],
[2.99999 :+ 0.0, 1.0 :+ 0.0, 1.00001 :+ 0.0],
[3.0 :+ 0.0, 1.0 :+ 0.0, 1.0 :+ 0.0],
[3.0 :+ 0.0, 1.0 :+ 0.0, 1.0 :+ 0.0],
[3.0 :+ 0.0, 1.0 :+ 0.0, 1.0 :+ 0.0]]
</pre>
Here is another example: this is a complex matrix and it is not
even hermitian. Yet, the algorithm still works, although its
fluctuates around true values.
<pre>
| 2-i 0 i |
| 0 1+i 0 |
| i 0 2-i |
[[2.0 :+ (-1.33333), 1.0 :+ 1.0, 2.0 :+ (-0.666667)],
[1.89245 :+ (-1.57849), 1.0 :+ 1.0, 2.10755 :+ (-0.421509)],
[1.81892 :+ (-1.80271), 1.0 :+ 1.0, 2.18108 :+ (-0.197289)],
[1.84565 :+ (-1.99036), 1.0 :+ 1.0, 2.15435 :+ (-0.00964242)],
[1.93958 :+ (-2.07773), 1.0 :+ 1.0, 2.06042 :+ 0.0777281],
[2.0173 :+ (-2.06818), 1.0 :+ 1.0, 1.9827 :+ 0.0681793],
[2.04357 :+ (-2.02437), 1.0 :+ 1.0, 1.95643 :+ 0.0243654],
[2.03375 :+ (-1.99072), 1.0 :+ 1.0, 1.96625 :+ (-0.00928429)],
[2.01245 :+ (-1.97875), 1.0 :+ 1.0, 1.98755 :+ (-0.0212528)],
[1.99575 :+ (-1.98307), 1.0 :+ 1.0, 2.00425 :+ (-0.0169263)],
[1.98938 :+ (-1.99359), 1.0 :+ 1.0, 2.01062 :+ (-0.00640583)],
[1.99145 :+ (-2.00213), 1.0 :+ 1.0, 2.00855 :+ 0.00212504],
[1.9968 :+ (-2.00535), 1.0 :+ 1.0, 2.0032 :+ 0.00535265],
[2.00108 :+ (-2.00427), 1.0 :+ 1.0, 1.99892 :+ 0.0042723],
[2.00268 :+ (-2.00159), 1.0 :+ 1.0, 1.99732 :+ 0.00158978],
[2.00213 :+ (-1.99946), 1.0 :+ 1.0, 1.99787 :+ (-0.000541867)],
[2.00079 :+ (-1.99866), 1.0 :+ 1.0, 1.9992 :+ (-0.00133514)],
[1.99973 :+ (-1.99893), 1.0 :+ 1.0, 2.00027 :+ (-0.00106525)],
[1.99933 :+ (-1.9996) , 1.0 :+ 1.0, 2.00067 :+ (-0.000397997)],
[1.99947 :+ (-2.00013), 1.0 :+ 1.0, 2.00053 :+ 0.000134972]]
The true eigenvalues are
2 - 2i, 1 + i, 2
</pre>
Some nonsymmetric real matrices have all real eigenvalues and
our algorithm still works for such cases. Here is one
such an example, which traditionally would have to be treated
by one of the Lanczos-like algorithms, specifically designed
for nonsymmetric real matrices. Evaluation of
<br>
<i>
iterated_eigenvalues [[2,1,1],[-2,1,3],[3,1,-1::Double]] 20
</i>
<br>
gives the following results
<pre>
[[3.0, -0.70818,-0.291815],
[3.06743, -3.41538, 2.34795],
[3.02238, -1.60013, 0.577753],
[3.00746, -2.25793, 1.25047],
[3.00248, -1.88764, 0.885154],
[3.00083, -2.06025, 1.05943],
[3.00028, -1.97098, 0.970702],
[3.00009, -2.0148, 1.01471],
[3.00003, -1.99268, 0.992648],
[3.00001, -2.00368, 1.00367],
[3.0, -1.99817, 0.998161],
[3.0, -2.00092, 1.00092],
[3.0, -1.99954, 0.99954],
[3.0, -2.00023, 1.00023],
[3.0, -1.99989, 0.999885],
[3.0, -2.00006, 1.00006],
[3.0, -1.99997, 0.999971],
[3.0, -2.00001, 1.00001],
[3.0, -1.99999, 0.999993],
[3.0, -2.0, 1.0]]
The true eigenvalues are:
3, -2, 1
</pre>
Finally, here is a case of a nonsymmetric real matrix with
complex eigenvalues:
<pre>
| 2 -3 |
| 1 0 |
</pre>
The direct application of "iterated_eigenvalues" would
fail to produce expected eigenvalues:
<pre>
1 + i sqrt(2) and 1 - i sqrt (2)
</pre>
But if we first precondition the matrix by adding "i" to its diagonal:
<pre>
| 2+i -3|
| 1 i|
</pre>
and then compute its iterated eigenvalues:
<br>
<i>
iterated_eigenvalues [[2:+1,1],[-3,0:+1]] 20
</i>
<br>
then the method will succeed. Here are the results:
<pre>
[[1.0 :+ 1.66667, 1.0 :+ 0.333333 ],
[0.600936 :+ 2.34977, 1.39906 :+ (-0.349766)],
[0.998528 :+ 2.59355, 1.00147 :+ (-0.593555)],
[1.06991 :+ 2.413, 0.93009 :+ (-0.412998)],
[1.00021 :+ 2.38554, 0.99979 :+ (-0.385543)],
[0.988004 :+ 2.41407, 1.012 :+ (-0.414074)],
[0.999963 :+ 2.41919, 1.00004 :+ (-0.419191)],
[1.00206 :+ 2.41423, 0.99794 :+ (-0.414227)],
[1.00001 :+ 2.41336, 0.99999 :+ (-0.413361)],
[0.999647 :+ 2.41421, 1.00035 :+ (-0.414211)],
[0.999999 :+ 2.41436, 1.0 :+ (-0.41436) ],
[1.00006 :+ 2.41421, 0.99993 :+ (-0.414214)],
[1.0 :+ 2.41419, 1.0 :+ (-0.414188)],
[0.99999 :+ 2.41421, 1.00001 :+ (-0.414213)],
[1.0 :+ 2.41422, 1.0 :+ (-0.414218)],
[1.0 :+ 2.41421, 0.99999 :+ (-0.414213)],
[1.0 :+ 2.41421, 1.0 :+ (-0.414212)],
[1.0 :+ 2.41421, 1.0 :+ (-0.414213)],
[1.0 :+ 2.41421, 1.0 :+ (-0.414213)],
[1.0 :+ 2.41421, 1.0 :+ (-0.414213)]]
</pre>
After subtracting "i" from the last result, we will get
what is expected.
<p>
<hr>
<p>
<b>
Eigenvectors for distinct eigenvalues
</b>
<p>
Assuming that eigenvalues of matrix A are already found
we may now attempt to find the corresponding aigenvectors
by solving the following homogeneous equation
<pre>
(A - a I) | x > = 0
</pre>
for each eigenvalue "a". The matrix
<pre>
B = A - a I
</pre>
is by definition singular, but in most cases it can be
triangularized by the familiar "factors_QR" procedure.
<pre>
B | x > = Q R | x > = 0
</pre>
It follows that the unknown eigenvector | x > is one of
the solutions of the homogeneous equation:
<pre>
R | x > = 0
</pre>
where R is a singular, upper triangular matrix with at least one
zero on its diagonal.
<dd>
If | x > is a solution we seek, so is its scaled version
alpha | x >. Therefore we have some freedom of scaling choice.
Since this is a homogeneous equation, one of the components
of | x > can be freely chosen, while the remaining components
will depend on that choice.
</pre>
To solve the above, we will be working from the bottom up of
the matrix equation, as illustrated in the example below:
<pre>
| 0 1 1 3 | | x1 |
| 0 1 1 2 | | x2 | /\
| 0 0 2 4 | | x3 | = 0 ||
| 0 0 0 0 | | x4 | ||
</pre>
Recall that the diagonal elements of any triangular matrix
are its eigenvalues.
Our example matrix has three distinct eigenvalues:
0, 1, 2. The eigenvalue 0 has degree of degeneration two.
Presence of degenerated eigenvalues complicates
the solution process. The complication arises when we have to
make our decision about how to solve the trivial scalar equations
with zero coefficients, such as
<pre>
0 * x4 = 0
</pre>
resulting from multiplication of the bottom row by vector | x >.
Here we have two choices: "x4" could be set to 0, or to any
nonzero number 1, say. By always choosing the "0" option
we might end up with the all-zero trivial vector -- which is
obviously not what we want. Persistent choice of the "1" option,
might lead to a conflict between some of the equations, such as
the equations one and four in our example.
<p>
So the strategy is as follows.
<p>
If there is at least one zero on the diagonal, find the topmost
row with zero on the diagonal and choose for it the solution "1".
Diagonal zeros in other rows would force the solution "0".
If the diagonal element is not zero than simply solve
an arithmetic equation that arises from the substitutions of
previously computed components of the eigenvector. Since certain
inaccuracies acumulate during QR factorization, set to zero all
very small elements of matrix R.
<p>
By applying this strategy to our example we'll end up with the
eigenvector
<pre>
< x | = [1, 0, 0, 0]
</pre>
<p>
If the degree of degeneration of an eigenvalue of A is 1 then the
corresponding eigenvector is unique -- subject to scaling.
Otherwise an eigenvector found by this method is one of many
possible solutions, and any linear combination of such solutions
is also an eigenvector. This method is not able to find more than one
solution for degenerated eigenvalues. An alternative method, which
handles degenerated cases, will be described in the next section.
<p>
The function below calculates eigenvectors corresponding to
distinct selected eigenvalues of any square matrix A, provided
that the singular matrix B = A - a I can still be factorized as Q R,
where R is an upper triangular matrix.
<pre>
> eigenkets a u
> --
> -- List of eigenkets of a square matrix A
> -- where
> -- a is a list of columns of A
> -- u is a list of eigenvalues of A
> -- (This list does not need to be complete)
> --
> | null u = []
> | not (null x') = x':(eigenkets a (tail u))
> | otherwise = (eigenket_UT (reverse b) d []):(eigenkets a (tail u))
> where
> a' = add_to_diagonal (-(head u)) a
> x' = unit_ket a' 0 (length a')
> b = snd (factors_QR a')
> d = discriminant [head bk | bk <- b] 1
> discriminant u n
> | null u = []
> | otherwise = x : (discriminant (tail u) m)
> where
> (x, m)
> | (head u) == 0 = (n, 0)
> | otherwise = (n, n)
> eigenket_UT b d xs
> | null b = xs
> | otherwise = eigenket_UT (tail b) (tail d) (x:xs)
> where
> x = solve_row (head b) (head d) xs
>
> solve_row u n x
> | almostZero p = n
> | otherwise = q/p
> where
> p = head u
> q
> | null x = 0
> | otherwise = -(sum_product (tail u) x)
>
> unit_ket a' m n
> | null a' = []
> | all (== 0) (head a') = unit_vector m n
> | otherwise = unit_ket (tail a') (m+1) n
</pre>
<p>
<hr>
<p>
<b>
Eigenvectors for degenerated eigenvalues
</b>
<p>
Few facts:
<ul>
<li>
Eigenvectors of a general matrix A, which does not have any
special symmetry, are not generally orthogonal. However, they
are orthogonal, or can be made orthogonal, to another set of
vectors that are eigenvectors of adjoint matrix A<sup>+</sup>;
that is the matrix obtained by complex conjugation and transposition
of matrix A.
<li>
Eigenvectors corresponding to nondegenerated eigenvalues of
hermitian or symmetric matrix are orthogonal.
<li>
Eigenvectors corresponding to degenerated eigenvalues are - in
general - neither orthogonal among themselves, nor orthogonal
to the remaining eigenvectors corresponding to other
eigenvalues. But since any linear combination of such degenerated
eigenvectors is also an eigenvector, we can orthogonalize
them by Gram-Schmidt orthogonalization procedure.
</ul>
Many practical applications deal solely with hermitian
or symmetric matrices, and for such cases the orthogonalization
is not only possible, but also desired for variety of reasons.
<dd>
But the method presented in the previous section is not able
to find more than one eigenvector corresponding to a degenerated
eigenvalue. For example, the symmetric matrix
<pre>
| 7 -2 1 |
A = | -2 10 -2 |
| 1 -2 7 |
</pre>
has two distinct eigenvalues: 12 and 6 -- the latter
being degenerated with degree of two. Two corresponding
eigenvectors are:
<pre>
< x1 | = [1, -2, 1] -- for 12
< x2 | = [1, 1, 1] -- for 6
</pre>
It happens that those vectors are orthogonal, but this is
just an accidental result. However, we are missing a third
distinct eigenvector. To find it we need another method.
One possibility is presented below and the explanation
follows.
<pre>
> eigenket' a alpha eps x' =
> --
> -- Eigenket of matrix A corresponding to eigenvalue alpha
> -- where
> -- a is a list of columns of matrix A
> -- eps is a trial inaccuracy factor
> -- artificially introduced to cope
> -- with singularities of A - alpha I.
> -- One might try eps = 0, 0.00001, 0.001, etc.
> -- x' is a trial eigenvector
> --
> scaled [xk' - dk | (xk', dk) <- zip x' d]
> where
> b = add_to_diagonal (-alpha*(1+eps)) a
> d = one_ket_solution b y
> y = matrix_ket (transposed b) x'
</pre>
Let us assume a trial vector | x' >, such that
<pre>
| x' > = | x > + | d >
where
| x > is an eigenvector we seek
| d > is an error of our estimation of | x >
</pre>
We first form a matrix B, such that:
<pre>
B = A - alpha I
</pre>
and multiply it by the trial vector | x' >, which
results in a vector | y >
<pre>
B | x' > = |y >
</pre>
On another hand:
<pre>
B | x' > = B | x > + B | d > = B | d >
because
B | x > = A | x > - alpha | x > = 0
</pre>
Comparing both equations we end up with:
<pre>
B | d > = | y >
</pre>
that is: with the system of linear equations for unknown error | d >.
Finally, we subtract error | d > from our trial vector | x' >
to obtain the true eigenvector | x >.
<p>
But there is some problem with this approach: matrix B is
by definition singular, and as such, it might be difficult
to handle. One of the two processes might fail, and their failures
relate to division by zero that might happen during either the
QR factorization, or the solution of the triangular system of equations.
<p>
But if we do not insist that matrix B should be exactly singular,
but almost singular:
<pre>
B = A - alpha (1 + eps) I
</pre>
then this method might succeed. However, the resulting eigenvector
will be the approximation only, and we would have to experiment
a bit with different values of "eps" to extrapolate the true
eigenvector.
<p>
The trial vector | x' > can be chosen randomly, although some
choices would still lead to singularity problems. Aside from
this, this method is quite versatile, because:
<ul>
<li>
Any random vector | x' > leads to the same eigenvector
for nondegenerated eigenvalues,
<li>
Different random vectors | x' >, chosen for degenerated
eigenvalues, produce -- in most cases -- distinct eigenvectors.
And this is what we want. If we need it, we can the always
orthogonalize those eigenvectors either internally (always
possible) or externally as well (possible only for hermitian
or symmetric matrices).
</ul>
It might be instructive to compute the eigenvectors for
the examples used in demonstration of computation of eigenvalues.
We'll leave to the reader, since this module is already too obese.
<p>
<hr>
<p>
<b>
Auxiliary functions
</b>
<p>
The functions below are used in the main algorithms of
this module. But they can be also used for testing. For example,
the easiest way to test the usage of resources is to use easily
definable unit matrices and unit vectors, as in:
<pre>
one_ket_solution (unit_matrix n::[[Double]])
(unit_vector 0 n::[Double])
where n = 20, etc.
> unit_matrix :: Num a => Int -> [[a]]
> unit_matrix m =
> --
> -- Unit square matrix of with dimensions m x m
> --
> [g 0 k | k <- [0..(m-1)]]
> where
> g i k
> | i == m = []
> | i == k = 1:(g (i+1) k)
> | otherwise = 0:(g (i+1) k)
>
> unit_vector :: Num a => Int -> Int -> [a]
> unit_vector i m =
> --
> -- Unit vector of length m
> -- with 1 at position i, zero otherwise
> [g i k| k <- [0..(m-1)]]
> where
> g i k
> | i == k = 1
> | otherwise = 0
> diagonals :: [[a]] -> [a]
> diagonals a =
> --
> -- Vector made of diagonal components
> -- of square matrix a
> --
> diagonals' a 0
> where
> diagonals' a n
> | null a = []
> | otherwise = (head (drop n (head a)))
> :(diagonals' (tail a) (n+1))
</pre>
<pre>
-----------------------------------------------------------------------------
--
-- Copyright:
--
-- (C) 1998 Numeric Quest Inc., All rights reserved
--
-- Email:
--
-- jans@numeric-quest.com
--
-- License:
--
-- GNU General Public License, GPL
--
-----------------------------------------------------------------------------
</pre>
</ul>
</body>
<SCRIPT language="Javascript">
<!--
// FILE ARCHIVED ON 20010628005806 AND RETRIEVED FROM THE
// INTERNET ARCHIVE ON 20030626101500.
// JAVASCRIPT APPENDED BY WAYBACK MACHINE, COPYRIGHT INTERNET ARCHIVE.
// ALL OTHER CONTENT MAY ALSO BE PROTECTED BY COPYRIGHT (17 U.S.C.
// SECTION 108(a)(3)).
var sWayBackCGI = "http://web.archive.org/web/20010628005806/";
function xLateUrl(aCollection, sProp) {
var i = 0;
for(i = 0; i < aCollection.length; i++)
if (aCollection[i][sProp].indexOf("mailto:") == -1 &&
aCollection[i][sProp].indexOf("javascript:") == -1)
aCollection[i][sProp] = sWayBackCGI + aCollection[i][sProp];
}
if (document.links) xLateUrl(document.links, "href");
if (document.images) xLateUrl(document.images, "src");
if (document.embeds) xLateUrl(document.embeds, "src");
if (document.body && document.body.background)
document.body.background = sWayBackCGI + document.body.background;
//-->
</SCRIPT>
</html>