------------------------------------------------------------------------------
-- Haskell module:      LinearAlgorithms
-- Date:                initialized 2001-03-25, last modified 2001-04-01
-- Author:              Jan Skibinski, Numeric Quest Inc.
-- Location:            http://www.numeric-quest.com/haskell/LinearAlgorithms.hs
-- See also:            http://www.numeric-quest.com/haskell/Orthogonals.html
--
-- Description:
-- This module provides several _selected_ linear algebra algorithms,
-- supporting computation of eigenvalues and eigenvectors of dense
-- matrices of small size. This module is to be utilized by module
-- Eigensystem, which redefines the eigenproblems in terms of
-- linear operators (maps) and abstract Dirac vectors.

-- Here is a list of implemented algorithms:
--
-- + triangular         A => R          where R is upper triangular
-- + triangular2        A => (R, Q)     such that R = Q' A Q
--
-- + tridiagonal        H => T          where H is Hermitian and T is
-- + tridiagonal2       H => (T, Q)     tridiagonal, such that T = Q' H Q
--
-- + subsAnnihilator    A => Q  such that Q A has zeroed subdiagonals
-- + reflection         x => y  where y is a complex reflection of x
--
-- Other algoritms, such as solution of linear equations are, at this time,
-- imported from module Orthogonals. The latter also deals with triangulization,
-- so you can compare the results from two different approaches:
-- orthogonalization vs. Householder reduction used in this module.
-- In essence the former method is a bit faster but overflows for large
-- number of iterations since, for typing reasons - its algorithms
-- avoid the normalization of vectors.
-- For full documentation of this module, and for references and the license,
-- go to the bottom of the page.
----------------------------------------------------------------------------

module LinearAlgorithms (
        triangular,
        triangular2,
        tridiagonal,
        tridiagonal2,
        Scalar,) where

import Complex
import Orthogonals hiding (Scalar)

type Scalar = Complex Double

----------------------------------------------------------------------------
-- Category: Iterative triangularization
--
--   triangular         A => R          where R is upper triangular
--   triangular2        A => (R, Q)     such that R = Q' A Q
----------------------------------------------------------------------------

mult :: [[Scalar]] -> [[Scalar]] -> [[Scalar]]
a `mult` b
    --  A matrix-product of matrices 'a' and 'b'
    --          C = A B
    --  where all matrices are represented as lists
    --  of scalar columns
        = matrix_matrix' (transposed a) b

triangular :: Int -> [[Scalar]] -> [[Scalar]]
triangular n a
    --  A (hopefully) triangular matrix R = Q' A Q obtained by
    --  'n' similarity transformations S(k) of matrix A:
    --          Q = S1 S2 S3 ....
    --
    -- If matrix A is Hermitian then the result is close
    -- to a diagonal matrix for sufficiently large n.
    | n == 0    = a
    | otherwise = triangular (n - 1) a1
    where
        a1  = (q' `mult` a ) `mult` q
        q'  = subsAnnihilator 0 a
        q   = adjoint q'


triangular2 :: Int -> [[Scalar]] -> ([[Scalar]], [[Scalar]])
triangular2 n a
    --  A pair of matrices (R, Q) obtained by 'n'
    --  similarity transformations, where R = Q' A Q
    --  is a (hopefully) triangular matrix, or diagonal
    --  if A is Hermitian. The transformation matrix Q
    --  is required for computation of eigenvectors
    --  of A.
    = triangular2' n a (unit_matrix n)
    where
        triangular2' o b p
            | o == 0    = (b, p)
            | otherwise = triangular2' (o - 1) b1 p1
            where
                b1 = (q' `mult` b ) `mult` q
                p1 = p `mult` q
                q' = subsAnnihilator 0 b
                q  = adjoint q'


----------------------------------------------------------------------------
-- Category: Tridiagonalization of a Hermitian matrix
--
-- + tridiagonal        H -> T  where H is Hermitian and T is tridiagonal
-- + tridiagonal2       H -> (T, Q)     such that T = Q' H Q
----------------------------------------------------------------------------


tridiagonal :: [[Scalar]] -> [[Scalar]]
tridiagonal h
    --  A tridiagonal matrix T = Q' H Q, obtained from Hermitian
    --  matrix H by a finite number of elementary similarity
    --  transformations (Householder reductions).
    | n < 3             = h
    | otherwise         = f (tail es) h 1
    where
        n       = length h
        es      = unit_matrix n

        f bs a k
            | length bs == 1    = a
            | otherwise         = f (tail bs)  a1 (k+1)
            where
                a1      = (q' `mult` a) `mult` q
                q'      = [r e | e <- es]
                q       = adjoint q'
                r       = reflection u (head bs)
                u       = replicate k 0 ++ drop k (a!!(k-1))


tridiagonal2 :: [[Scalar]] -> ([[Scalar]], [[Scalar]])
tridiagonal2 h
    --  A pair (T, Q) of matrices, obtained from
    --  similarity transformation of Hermitian matrix H
    --  where T = Q' H Q is a tridiagonal matrix and Q is unitary
    --  transformation made of a finite product of
    --  elementary Householder reductions.
    | n < 3             = (h, es)
    | otherwise         = f (tail es) h es 1
    where
        n       = length h
        es      = unit_matrix n

        f bs a p k
            | length bs == 1    = (a, p)
            | otherwise         = f (tail bs) a1 p1 (k+1)
            where
                a1      = (q' `mult` a) `mult` q
                q'      = [r e | e <- es]
                q       = adjoint q'
                p1      = p `mult` q
                r       = reflection u (head bs)
                u       = replicate k 0 ++ drop k (a!!(k-1))


----------------------------------------------------------------------------
-- Category: Elementary unitary transformations
--
-- + subsAnnihilator    A => Q  such that Q A has zeroed subdiagonals
-- + reflection         x => y  where y is a complex reflection of x
----------------------------------------------------------------------------

subsAnnihilator :: Int -> [[Scalar]] -> [[Scalar]]
subsAnnihilator k a
    --  A unitary matrix Q' transforming any n x n
    --  matrix A to an upper matrix B, which has
    --  zero values below its 'k'-th subdiagonal
    --  (annihilates all subdiagonals below k-th)
    --          B = Q' A
    --  where
    --      'a' is a list of columns of matrix A
    --
    --  If k=0 then B is an upper triangular matrix,
    --  if k=1 then B is an upper Hessenberg matrix.
    --  The transformation Q is built from n - k - 1
    --  elementary Householder transformations of
    --  the first n-k-1 columns of iteratively transformed
    --  matrix A.
    | n < 2 + k         = es
    | otherwise         = f (drop k es) a1 es k
    where
        n       = length a
        es      = unit_matrix n
        a1      = take (n - 1 - k) a

        f bs b p l
            | length bs == 1    = p
            | otherwise         = f (tail bs)  b1 p1 (l+1)
            where
                b1      = [r v |v <- tail b]
                p1      = q' `mult` p
                q'      = [r e | e <- es]
                r       = reflection u (head bs)
                u       = replicate k 0 ++ drop l (head b)


reflection :: [Scalar] -> [Scalar] -> [Scalar] -> [Scalar]
reflection a e x
    --  A vector resulting from unitary complex
    --  Householder-like transformation of vector 'x'.
    --
    --  The operator of such transformation is defined
    --  by mapping vector 'a' to a multiple 'p' of vector 'e'
    --          U |a > = p | e >
    --  where scalar 'p' is chosen to guarantee unitarity
    --          < a | a > = < p e | p e>.
    --
    --  This transformation is not generally Hermitian, because
    --  the scalar 'p' might become complex - unless
    --          < a | e > = < e | a >,
    --  which is the case when both vectors are real, and
    --  when this transformation becomes a simple Hermitian
    --  reflection operation.
    --  See reference [1] for details.
    --
    | d == 0    = x
    | otherwise = [xk - z * yk |(xk, yk) <- zip x y]
    where
        z = s * bra_ket y x
        s = 2/h :+ (-2 * g)/h
        h = 1 + g^(2::Int)
        g = imagPart a_b / d
        d = a_a - realPart a_b
        y = normalized [ak - bk |(ak, bk) <- zip a b]
        p = a_a / (realPart (bra_ket e e))
        b = map ((sqrt p :+ 0) * ) e
        a_a = realPart (bra_ket a a)
        a_b = bra_ket a b

----------------------------------------------------------------------------
-- Category: Test data
--
----------------------------------------------------------------------------

-- matrixA :: [[Scalar]]
-- matrixA
--     --  Test matrix A represented as list of scalar columns.
--     =   [
--                 [1, 2, 4, 1, 5]
--         ,       [2, 3, 2, 6, 4]
--         ,       [4, 2, 5, 2, 3]
--         ,       [1, 6, 2, 7, 2]
--         ,       [5, 4, 3, 2, 9]
--         ]

----------------------------------------------------------------------------
-- Module documentation
-- ====================

-- Representation of vectors, matrices and scalars:
-- ------------------------------------------------
-- We have chosen to follow the same scheme as used in module Orthogonals:
-- vectors are represented here as lists of scalars, while matrices --
-- as lists of scalar columns (vectors). But while scalars over there are
-- generic and cover a range of types, the scalars of this module are
-- implemented as Complex Double. Although all algorithms here
-- operate on complex matrices and complex vectors, they will work
-- on real matrices without modifications. If however, the performance
-- is a premium it will be a trivial exercise to customize all these
-- algorithms to real domain. Perhaps the most important change should
-- be then made to a true workhorse of this module, the function 'reflection',
-- in order to convert it to a real reflection of a vector in a hyperplane
-- whose normal is another vector.
--
-- Schur triangularization of any matrix:
-- --------------------------------------
-- The Schur theorem states that there exists a unitary matrix Q such
-- that any nonsingular matrix A can be transformed to an upper triangular
-- matrix R via similarity transformation
--      R = Q' A Q
-- which preserves the eigenvalues. Here Q' stands for a Hermitian
-- conjugate of Q (adjoint, or Q-dagger).

-- Since the eigenvalues of a triangular matrix R are its diagonal
-- elements, finding such transformation solves the first part of
-- the eigenproblem. The second part, finding the eigenvectors of A,
-- is trivial since they can be computed from eigenvectors of R:
--      | x(A) > = Q | x(R) >
--
-- In particular, when matrix A is Hermitian, then the matrix R
-- becomes diagonal, and the eigenvectors of R are its normalized
-- columns; that is, the unit vectors. It follows that the eigenvectors
-- of A are then the columns of matrix Q.
-- But when A is not Hermitian one must first find the eigenvectors
-- of a triangular matrix R before applying the above transformation.
-- Fortunately, it is easier to find eigenvectors of a triangular matrix
-- R than those of the square matrix A.
--
-- Implementation of Schur triangularization via series of QR factorizations:
-- --------------------------------------------------------------------------
-- The methods known in literature as QR factorization (decomposition)
-- methods iteratively compose such unitary matrix Q from a series of
-- elementary unitary transformations, Q(1), Q(2)..:
--      Q = Q(1) Q(2) Q(3) ...
-- The most popular method of finding those elementary unitary
-- transformations relies on a reflection transformation, so selected as
-- to zero out all components of the matrix below its main diagonal. Our
-- implementation uses a complex variety of such a 'reflection', described
-- in the reference [1]. The columnar reduction of the lower portion of
-- the matrix to zeros is also known under the name of Householder
-- reduction, or Householder transformation. This is, however, not the
-- only possible choice for elementary transformations; see for example
-- our module Orthogonals, where such transformations are perfomed via
-- Gram-Schmidt orthogonalization procedure instead.
--
-- The iterative functions 'triangular' and 'triangular2' attempt to
-- triangularize any complex matrix A by a series of similarity
-- transformation, known in literature as QR decomposition.
-- Function 'triangular' does not deliver the transformation Q but
-- only a transformed matrix A, which should be close to triangular
-- form after a sufficient number of iterations. Use this function
-- if you are interested in eigenvalues only. But when you need
-- the eigenvectors as well, then use the function 'triangular2',
-- which also delivers the transformation Q, as shown below:
--   triangular         A => R  where R is upper triangular
--   triangular2        A => (R, Q)     such that R = Q' A Q
--
-- Tridiagonalization of Hermitian matrices:
-- -----------------------------------------
-- While the above functions are iterative and require a bit of
-- experimentation with a count of iterations to figure out whether
-- the required accuracy has yet been achieved, the tridiagonalization
-- methods transform any matrix A to a tridiagonal form in a finite
-- number of elementary transformations.
--
-- However, our implementation is not generic because it performs
-- tridiagonalization only on Hermitian matrices. It uses the same
-- unitary 'reflection', as the triangularization does.
--
-- Why would you care for such tridiagonalization at all? Many world
-- class algorithms use it as a first step to precondition the original
-- matrix A for faster convergence and for better stability and accuracy.
-- Its cost is small in comparison to the overall cost incurred during
-- the iterative stage. What's more, the triangularization iteration
-- does preserve the shape of tridiagonal matrix at each step - bringing
-- it only closer to the diagonal shape. So the tridiagonalization
-- is a recommended option to be executed before the iterative
-- triangulariation.
--
-- Again, we are offering here two versions of the tridiagonalization:
--
-- + tridiagonal        H -> T  where H is Hermitian and T is tridiagonal
-- + tridiagonal2       H -> (T, Q)     such that T = Q' H Q
--
-- Elementary transformations:
-- ---------------------------
-- All the above algorithms heavily rely on the function 'reflection'
-- which defines a complex reflection transformation of a vector. One use
-- of this function is to perform a Householder reduction of a column-vector,
-- to zero out all of its components but one. For example, the unitary
-- transformation 'subsAnnihilator 0' annihilates all subdiagonals lying
-- below the main diagonal. Similarly, 'subsAnnihilator 1' would zero out
-- all matrix components below its first subdiagonal - leading to a so-called
-- upper Hessenberg matrix.
--
-- + subsAnnihilator    A => Q  such that Q A has zeroed subdiagonals
-- + reflection         x => y  where y is a complex reflection of x
--
----------------------------------------------------------------------------
-- References:
-- [1]  Xiaobai Sun, On Elementary Unitary and Phi-unitary transformations,
--      Duke University, Department Of Computer Science, 1995,
--      http://citeseer.nj.nec.com/340881.html
---------------------------------------------------------------------------
--
-- Copyright:
--
--      (C) 2001 Numeric Quest, All rights reserved
--
--      Email: jans@numeric-quest.com
--
--      http://www.numeric-quest.com
--
-- License:
--
--      GNU General Public License, GPL
--
---------------------------------------------------------------------------