------------------------------------------------------------------------------
-- Haskell module: Eigensystem
-- Date: initialized 2001-03-25, last modified 2001-03-25
-- Author: Jan Skibinski, Numeric Quest Inc.
-- Location: http://www.numeric-quest.com/haskell/Eigensystem.hs
-- See also: http://www.numeric-quest.com/haskell/QuantumVector.html
-- See also: http://www.numeric-quest.com/haskell/Orthogonals.html
--
-- Description:
--
-- This module extends the QuantumVector module by providing functions
-- to calculate eigenvalues and eigenvectors of Hermitian operators.
-- Such toolkit is of primary importance due to pervasiveness of
-- eigenproblems in Quantum Mechanics.
--
-- This module is organized in three layers:
--
-- 1. Interface to module QuantumVector, where all function signatures
-- are expressed in terms of linear operators, Dirac vectors and scalars.
--
-- Here the operators are defined directly via maps from input to
-- output vectors. In many cases it is much easier to define the operators
-- directly rather than to rely on their matrix representation.
--
-- 2. Conversion layer between operators and their matrix representation.
--
-- Sometimes it is more convenient to start with an underlying matrix
-- representation of an operator. There are also cases where a direct
-- manipulation on operators is too difficult, while it is trivial
-- to obtain the corresponding results via matrices. One example is a
-- computation of a Hermitian conjugate of A:
-- < ei | A' | ej > = conjugate < ej | A | ej >
-- (Here ' stands for a dagger)
-- If however the operator A is made from a product or a sum of simpler
-- operators, whose Hermitian conjugates are known to us, then the
-- direct approach from the upper layer could be easier and perhaps more
-- efficient in some cases.
--
-- 3. Implementation layer is stored in a separate module LinearAlgorithms,
-- where matrices are represented as lists of columns of scalars, and
-- vectors -- as lists of scalars.
--
-- This layer is completely independendent of the other two and can be
-- reused separately for applications other than those caring for the
-- QuantumVector module and its notation. It can also be reimplemented
-- via Haskell arrays, or perhaps by some other means, such as trees
-- of nodes relating square blocks of data to support paralleism.
--
-- See also bottom of the page for references and license.
-----------------------------------------------------------------------------
module Eigensystem (eigenvalues, adjoint) where
import Complex
import QuantumVector
import LinearAlgorithms (triangular, tridiagonal, triangular2)
import List (findIndex)
----------------------------------------------------------------------------
-- Category: Eigensystem for QuantumVector
----------------------------------------------------------------------------
eigenvalues :: Ord a => Bool -> Int -> [Ket a] -> (Ket a -> Ket a) -> [Scalar]
eigenvalues doTri n es a
-- A list of eigenvalues of operator 'a'
-- obtained after 'n' triangularizations
-- of a matrix corresponding to operator 'a'
-- where
-- 'es' is a list of base vectors
-- 'doTri' declares whether or not we
-- want the initial tridiagonalization
-- (applies to Hermitian operators only)
| doTri == True = f b1
| otherwise = f b
where
f c = diagonals $ operator es $ triangular n c
diagonals us = [toBra e <> us e | e <- es]
b = matrix es a
b1 = tridiagonal b
eigenpairs :: Ord a => Int -> [Ket a] -> (Ket a -> Ket a) -> ([Scalar], [Ket a])
eigenpairs n es a
-- A pair of lists (eigenvalues, eigenvectors) of hermitian
-- operator 'a' obtained after 'n' triangularizations of 'a'
-- where
-- 'es' is a list of base vectors
-- Note: For a moment this applies only to Hermitian operators
-- until we decide what would be the best way to compute eigenvectors
-- of a triangular matrix: the method from module Orthogonal, power
-- iteration, etc.
= (ls, xs)
where
(t, q) = triangular2 n b
b = matrix es a
ls = [ tk!!k | (tk, k) <- zip t [0..length t - 1] ]
xs = [compose qk es | qk <- q]
adjoint :: Ord a => [Ket a] -> (Ket a -> Ket a) -> (Ket a -> Ket a)
adjoint es a
-- A Hermitian conjugate of operator a,
-- (or a-dagger, or adjoint to a)
-- where 'es' is a list of base vectors
= operator es ms
where
ms = [[ conjugate (toBra ei <> vj) | vj <- v] | ei <- es]
v = [a ej | ej <- es]
----------------------------------------------------------------------------
-- Category: Conversion from operators to matrices and vice versa
----------------------------------------------------------------------------
operator :: Ord a => [Ket a] -> [[Scalar]] -> Ket a -> Ket a
operator basis ms x
-- Definition of an operator corresponding
-- to a matrix 'ms' given as a list of scalar
-- columns
-- where
-- 'basis' is a complete list of base vectors
-- 'x' is any ket vector from this space
= a >< x
where
a u = case (findIndex (u == ) basis) of
Just k -> compose (ms !! k) basis
Nothing -> error "Out of bounds"
matrix :: Ord a => [Ket a] -> (Ket a -> Ket a) -> [[Scalar]]
matrix basis a
-- List of scalar columns representing
-- the operator 'a' in a given 'basis'
= [[ei' <> vj | ei' <- e'] | vj <- v]
where
v = [a ej | ej <- basis]
e' = [toBra ei | ei <- basis]
----------------------------------------------------------------------------
-- 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]
]
opA = operator basisA matrixA
basisA = map Ket [1..5::Int] -- or: map Ket "abcde", etc.
---------------------------------------------------------------------------
-- Copyright:
--
-- (C) 2001 Numeric Quest, All rights reserved
--
-- Email: jans@numeric-quest.com
--
-- http://www.numeric-quest.com
--
-- License:
--
-- GNU General Public License, GPL
--
---------------------------------------------------------------------------