Decompositions

Back to: Main contents.

Contents

The LU decomposition (Zludpp)
The Cholesky decomposition (Zchol)
The QR decomposition (Zqrd, Zhqrd)
The spectral decomposition of a Hermitian matrix (Zspec)
The singular value decomposition (Zsvd)
Hessenberg form (Zhess)
Schur form (Schur)
Eigenvalue-vector form (Eig)

A matrix decomposition is a factorization of a matrix into the product of simpler matrices. For example, Gaussian elimination with partial pivoting when applied to a matrix A produces a factorization of the form

A = PLU,
where P is a permutation matrix, L is a lower triangular matrix, and U is an upper triangular matrix. This decomposition is widely used in the solution of linear systems--especially by the Solve suite.

At the heart of most computations with dense matrices is a decomposition. Jampack provides classes implementing the more commonly used decompositions. The computation of the decomposition is done by the class constructor.

The uses of a particular decomposition are too varied to be exhausted by a fixed set of methods. Consequently, Jampack's decomposition classes (with the exception of Zhqrd) have no methods of their own, and its members are public for ready access.

Certain matrix decompositions (ludpp, chol, and hqrd) are distinguished by the fact that they can be saved for further use on the history list of the matrix. In particular the Solve suite makes heavy use of this feature. While a decomposition is on the history list, it is hidden from applications outside the Jampack package. For more see Parameters.History and Matrices.Zmat.

The LU decomposition (Zludpp)

Back to: Main contents, Top of section.

The LU decomposition with partial pivoting is the decomposition computed by Gaussian elimination with partial pivoting. Specifically, any matrix A can be factored in the form

A = PLU,
where

P is a permutation matrix,

L is a lower triangular matrix with ones on its diagonal and with its subdiagonal elements less than one in magnitude,

U is upper triangular.
This decomposition is widely used in computing solutions to linear systems.

The pivoted LU decomposition is implemented by the class Zludpp. The matrix can be a general mxn Zmat, although square matrices are the rule in most applications. If m>n, the L-factor is lower trapezoidal; if m<n, the U-factor is upper trapezoidal.

The action of the permutation matrix P is specified by an integer vector of pivot indices. For more see the Pivot suite.

The members of Zlupp are all public.

public int nrl
The number of rows in L

public int ncl
The number of columns in L

public int nru
The number of rows in U

public int ncu
The number of columns in U

public int pvt[]
The pivot array

public Zltmat L
The L factor

public Zutmat U
The U factor

Note that L and U are of class Zltmat and Zutmat, subclasses of Zmat that tag them as representing lower and upper triangular matrices.

The decomposition itself is computed in the constructor for Zlupp.

Zludpp(Zmat A)
Returns the Zludpp of A.

A caution. The pivoted LU decomposition is always well defined, hence Zludpp does not return and exception. However, if A is singular, the matrix U will be singular (or numerically will have small diagonal elements), which may cause problems later.

The Cholesky decomposition (Zchol)

Back to: Main contents, Top of section.

If A is a (Hermitian) positive semi-definite matrix, then A can be factored in the form

A = RHR,
where R is upper triangular. This factorization is called the Cholesky decomposition of A. The class Zchol implements the Cholesky decomposition.

The Cholesky decomposition is defined only if A is Hermitian and satisfies xTAx >= 0 for all x. The first condition is checked by Zchol before any significant calculations can begin. The second condition is signaled by a failure of the algorithm for computing the Cholesky decomposition. The failure of either condition causes a JampackExcpetion to be thrown.

Here are the specifications for the Cholesky decomposition.

public int n
The order of A and R

public Zutmat R
The Cholesky factor of A

public Zchol(Zmat A)
Returns the Cholesky decomposition of A.
Throws JampackException if A is not Hermitian.
Throws JampackExcpetion if the decomposition does not exits

The QR decomposition (Zqrd, Zhqrd)

Back to: Main contents, Top of section.

For any mxn matrix A there is a unitary matrix Q such

QHA = R,
where R is zero below its main diagonal. Jampack provides two versions of this QR decomposition. The first is an explicit version in which Q and R are returned as matrices. The second version represents Q as a product of Householder transformations.

The explicit QR decomposition of an matrix mxn A is computed by the class Zqrd. Its data fields and constructor are

public Zmat Q
The Q factor of A

public Zmat R
The R factor of A

public Zqrd(Zmat A)
Returns the QR decomposition of A.

The factored form of the QR decomposition of a mxn matrix A is computed by the class Zhqrd (the "h" stands for Householder). The matrix Q is represented as the product of min(m,n) Householder transformations, whose generating vectors are stored in an array of Z1s. Moreover, if m>n, the the matrix R is made square by stripping it of its trailing m-n rows. Here are the specifications for Zhqrd

public int nrow
The number of rows in A

public int ncol
The number of columns in A

public int ntran
the number of Householder transformations into which Q is factored

public int Z1[] U
An array containing the ntran vectors generating the Householder transformations.

public Zutmat R
The R-factor of A. If nrow>ncol then R is square of order ncol. Otherwise R has the same dimensions as A.

public Zhqrd(Zmat A)
Returns the Householder QR decomposition of A.

Because the Q-factor produced by Zhqrd is in factored form, the class provides four methods for multiplying Q and its conjugate transpose into a matrix.

public Zmat qb(Zmat B)
Returns QB.

public Zmat qhb(Zmat B)
Returns QHB.

public Zmat bq(Zmat B)
Returns BQ.

public Zmat bqh(Zmat B)
Returns BQH.

The spectral decomposition of a Hermitian matrix (Zspec )

Back to: Main contents, Top of section.

For any Hermitian matrix A, there is a unitary matrix U such that

UHAU = D,
where D is a real diagonal matrix. The columns of U are eigenvectors of A. The corresponding diagonal elements of D are their eigenvalues. The decomposition is called the spectral decomposition of A. This decomposition is implemented by the class Zspec.

public Zmat U
The unitary matrix of eigenvectors of A

public Zdiagmat D
The diagonal matrix of eigenvalues of A

public Zspec (Zmat A)
Returns the spectral decomposition of A.
Throws a JampackException if A is not Hermitian.

Note that when real matrices are implemented, D will become a Ddiagmat.

The singular value decomposition (Zsvd)

Back to: Main contents, Top of section.

Let A be an mxn matrix m>=n there are unitary matrices U and V such that 

UHAV = | S |
       | 0 |
where S = diag(s1,...,sn) with
s1 >= s2 >= ... >= sn >= 0.
If m<n the decomposition has the form 
UHAV = | S  0 |,
where S is diagonal of order n.

The diagonals of S are the singular values of A. The columns of U are the left singular vectors of A and the columns of V are the right singular vectors.

The singular value decomposition is implemented by the class Zsvd.

static int MAXITER = 30
Limits the number of iterations in the SVD algorithm.

public Zmat U
The matrix of left singular vectors of X

public Zmat V
The matrix of right singular vectors of X

public Zdiagmat S
The diagonal matrix of singular values of X

public Zsvd(Zmat A)
Returns the SVD of A.
Throws a JampackException if MAXITER is exceeded.

Hessenberg form (Zhess)

Back to: Main contents, Top of section.

If A is a matrix of order n, there is a a unitary matrix U such that

H = UHAU
is upper Hessenberg; that is, H has the form illustrated below for n = 6. 
    | X X X X X X |
    | X X X X X X |
H = | O X X X X X |.
    | O O X X X X |
    | O O O X X X |
    | O O O O X X |
This reduction to Hessenberg form is implemented by the class Zhess.

public Zmat H
The Hessenberg form of A

public Zmat U
The unitary transformation that reduces A to H

public Zhess(Zmat A)
Returns the Hessenberg form of A.
Throws a JampackException if A is not square.

Schur form (Schur)

Back to: Main contents, Top of section.

If A is a matrix of order n, there is a unitary matrix U such that

T = UHAU
is upper triangular. The matrix T is called a Schur form of A. The diagonals of T are eigenvalues of A and the corresponding columns of U are its Schur vectors. The decomposition is not unique--the order of the eigenvalues on the diagonal of T can vary. However, in most applications any Schur form is sufficient.

The class Schur implements this decomposition.

public Zmat T
A Schur form of A

public Zmat U
The unitary transformation that reduces A to T

public static MAXITER = 30
Limits the number of iterations in the QR algorithm

public Schur(Zmat A)
Returns a Schur form of A.
Throws a JampackException if A is not square.

Eigenvalue-vector form (Eig)

Back to: Main contents, Top of section.

Given a diagonalizable matrix A of order n, there is a nonsingular matrix X such that

D = X-1AX
is diagonal. The columns of X are eigenvectors of A corresponding to the diagonal elements of D. Eig implements this eigenvalue-vector decomposition.

public Zmat X
The matrix of eigenvectors of A

public Zdiagmat D
The matrix of eigenvalues of A

public Eig(Zmat A)
Returns an eigenvalue-vector decomposition of A.
Throws a JampackException if A is not square.
The user should be aware that the eigenvalue-vector decomposition is not necessarily stable. In particular, the matrix X of eigenvectors may be ill conditioned.