com.imsl.math

Class Cholesky

java.lang.Object

com.imsl.math.Cholesky

All Implemented Interfaces:

Serializable, Cloneable

public class Cholesky
extends Object
implements Serializable, Cloneable

Cholesky factorization of a matrix of type double.

Class Cholesky uses the Cholesky-Banachiewicz algortithm to factor the matrix A.

The Cholesky factorization of a matrix is \(A = RR^T\), where R is a lower triangular matrix. Thus, $$A = RR^T = \begin{pmatrix} R_{11} & 0 & 0\\ R_{21} & R_{22} & 0\\ R_{31} & R_{32} & R_{33} \end{pmatrix} \begin{pmatrix} R_{11} & R_{21} & R_{31}\\ 0 & R_{22} & R_{32}\\ 0 & 0 & R_{33} \end{pmatrix} = \begin{pmatrix} R^2_{11} & & (Symmetric) \\ R_{21}R_{11} & R^2_{22}+R^2_{22} & \\ R_{31}R_{11} & R_{31}R_{21}+R_{32}R_{22} & R^2_{31}+R^2_{32}+R^2_{33} \end{pmatrix} $$ which leads to the following for the entries of the lower triangular marix R: $$ R_{i,j} = \frac{1}{R_{j,j}}\left ( A_{i,j}-\sum_{k=1}^{j-1}R_{i,k}R_{j,k} \right )\mbox{,} \,\,\,\, \mbox{for}\,\,\, i>j $$ and $$ R_{i,i} = \sqrt{A_{i,i} - \sum_{k=1}^{i-1}R_{i,k}^2} \,\,\,\mbox{.} $$

The method update is based on the LINPACK routine SCHUD; see Dongarra et al. (1979) and updates the \(RR^T\) Cholesky factorization of the real symmetric positive definite matrix A after a rank-one matrix is added. Given this factorization, update computes the factorization \(\tilde R\) such that $$A + xx^T = \tilde R \tilde R^T \,\, \mbox{.}$$

Similarly, the method downdate, based on the LINPACK routine SCHDD; see Dongarra et al. (1979), downdates the \(RR^T\) Cholesky factorization of the real symmetric positive definite matrix A after a rank-one matrix is subtracted. downdate computes the factorization \(\tilde R\) such that $$A - xx^T = \tilde R \tilde R^T \,\, \mbox{.}$$ This is not always possible, since \(A - xx^T\) may not be positive definite.

downdate determines an orthogonal matrix U as the product \(G_N\ldots G_1\) of Givens rotations, such that $$U \left[ \begin{array}{l} R^T \\ 0 \\ \end{array} \right] = \left[ \begin{array}{l} {\tilde R}^T \\ x^T \\ \end{array} \right] $$

By multiplying this equation by its transpose and noting that \(U^TU = I\), the desired result $$RR^T - xx^T = \tilde R \tilde R^T $$ is obtained.

Let a be the solution of the linear system \(Ra = x\) and let $$\alpha = \sqrt {1 - \left\| a \right\|_2^2 }$$

The Givens rotations, \(G_i\), are chosen such that $$G_1 \cdots G_N \left[ \begin{array}{l}a \\ \alpha \end{array} \right] = \left[ \begin{array}{l} 0 \\ 1 \end{array} \right] $$

The \(G_i\), are (N + 1) by (N + 1) matrices of the form $$G_i = \left[ {\begin{array}{*{20}c}{ I_{i - 1} } & 0 & 0 & 0 \\ 0 & {c_i } & 0 & { - s_i } \\ 0 & 0 & {I_{N - i} } & 0 \\ 0 & {s_i } & 0 & {c_i } \\ \end{array}} \right] $$ where \(I_k\) is the identity matrix of order k; and \(c_i= \cos\theta _i, s_i= \sin\theta_i\) for some \(\theta_i\).

The Givens rotations are then used to form $$ {\tilde R}^T,\,\,G_1 \cdots \,G_N \left[ \begin{array}{l} R^T \\ 0 \\ \end{array} \right] = \left[ \begin{array}{l} {\tilde R}^T \\ {\tilde x^T} \\ \end{array} \right] $$

The matrix $$\tilde R$$ is lower triangular and $$\tilde x = x$$ because $$ x = \left( {R \,\,\, 0} \right) \left[ \begin{array}{l}a \\ \alpha \\ \end{array} \right] = \left( R \,\,\, 0 \right) U^TU \left[ \begin{array}{l}a \\ \alpha \\ \end{array} \right] = \left( {\tilde R \,\,\, \tilde x} \right) \left[ \begin{array}{l}0 \\ 1 \\ \end{array} \right] = \tilde x $$.

See Also:

Example, Serialized Form

Nested Class Summary

Nested Classes

Modifier and Type	Class and Description
static class	Cholesky.NotSPDException The matrix is not symmetric, positive definite.

Constructor Summary

Constructors

Constructor and Description
Cholesky(double[][] a) Create the Cholesky factorization of a symmetric positive definite matrix of type double.

Method Summary

Modifier and Type	Method and Description
void	downdate(double[] x) Downdates the factorization by subtracting a rank-1 matrix.
double[][]	getR() Returns the R matrix that results from the Cholesky factorization.
double[][]	inverse() Returns the inverse of this matrix
double[]	solve(double[] b) Solve Ax = b where A is a positive definite matrix with elements of type double.
void	update(double[] x) Updates the factorization by adding a rank-1 matrix.

Methods inherited from class java.lang.Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Constructor Detail

Cholesky

public Cholesky(double[][] a)
         throws SingularMatrixException,
                Cholesky.NotSPDException

Create the Cholesky factorization of a symmetric positive definite matrix of type double.

Parameters:

a - a double square matrix to be factored

Throws:

IllegalArgumentException - Thrown when the row lengths of matrix a are not equal (for example, the matrix edges are "jagged".)

SingularMatrixException - Thrown when the input matrix A is singular.

Cholesky.NotSPDException - Thrown when the input matrix is not symmetric, positive definite.

Method Detail

getR

public double[][] getR()

Returns the R matrix that results from the Cholesky factorization.

Returns:

a double matrix which contains the lower triangular R matrix that results from the Cholesky factorization such that \(A = RR^T\)

solve

public double[] solve(double[] b)

Solve Ax = b where A is a positive definite matrix with elements of type double.

Parameters:

b - a double array containing the right-hand side of the linear system

Returns:

a double array containing the solution to the system of linear equations

inverse

public double[][] inverse()

Returns the inverse of this matrix

Returns:

a double matrix containing the inverse

update

public void update(double[] x)

Updates the factorization by adding a rank-1 matrix. The object will contain the Cholesky factorization of \(A + xx^T\), where A is the previously factored matrix.

Parameters:

x - A double array which specifies the rank-1 matrix. x is not modified by this function.

downdate

public void downdate(double[] x)
              throws Cholesky.NotSPDException

Downdates the factorization by subtracting a rank-1 matrix. The object will contain the Cholesky factorization of \(A - xx^T\), where A is the previously factored matrix.

Parameters:

x - A double array which specifies the rank-1 matrix. x is not modified by this function.

Throws:

Cholesky.NotSPDException - if \(A - xx^T\) is not symmetric positive-definite.