Package 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:

  • Nested Class Summary

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

  • Constructor Summary

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

  • Method Summary

    Modifier and Type
    Method
    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 Details

  • Method Details

    • 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.