org/apache/commons/math/linear/CholeskySolver.java

/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.apache.commons.math.linear;


/**
 * Solves a linear equitation with symmetrical, positiv definit 
 * coefficient matrix by Cholesky decomposition.
 * <p>
 * For every symmetric, positiv definit matrix <code>M</code> there is a
 * lower triangular matrix <code>L</code> so that <code>L*L^T=M</code>. 
 * <code>L</code> is called the <i>Cholesky decomposition</i> of <code>M</code>.  
 * For any constant vector <code>c</code> it can be used to solve 
 * the linear equitation <code>M*x=L*(L^T*x)=c</code>.<br>
 * Compared to the LU-decompoistion the Cholesky methods requires only half 
 * the number of operations.
 * <p>
 * @author Stefan Koeberle, 11/2003
 */
public class CholeskySolver {
    
    private double numericalZero = 10E-12;
    
    /** The lower triangular matrix */
    private Array2DRowRealMatrix decompMatrix;
    
   
    /** 
     * Creates a new instance of CholeskySolver
     */
    public CholeskySolver() {
    }//constructor CholeskySolver
    
    
    /** 
     * Every double <code>d</code> satisfying 
     * <code>java.lang.Math.abs(d) <= numericalZero</code> 
     * is considered equal to <code>0.0d.</code>
     */
    public void setNumericalZero(double numericalZero) {
        this.numericalZero = numericalZero;
    }//setNumericalZero
    
    /**
     * See <code>setNumericalZero</code>
     */
    public double getNumericalZero() {
        return numericalZero;
    }//getNumericalZero
    
    
    /**
     * Calculates the Cholesky-decomposition of the symmetrical, positiv definit 
     * matrix <code>M</code>.
     * <p>
     * The decomposition matrix is internally stored.
     * <p>
     * @throws IllegalArgumentException   if <code>M</code> ist not square or 
     *                                    not positiv definit
     */
    public void decompose(RealMatrix m) 
    throws IllegalArgumentException {
       
       decompMatrix = null;
       double[][] mval = m.getData();
       int numRows = m.getRowDimension();
       int numCols = m.getColumnDimension();
       if (numRows != numCols) 
           throw new IllegalArgumentException("matrix is not square"); 
       double[][] decomp = new double[numRows][numCols];       
       double sum;
       
       //for all columns
       for (int col=0; col<numCols; col++) {
          
           //diagonal element
           sum = mval[col][col];
           for (int k=0; k<col; k++) 
               sum = sum - decomp[col][k]*decomp[col][k];
           if (sum <= numericalZero) {
               throw new IllegalArgumentException(
                             "Matrix is not positiv definit");
           }
           decomp[col][col] += Math.sqrt(sum);
           
           //column below diagonal
           for (int row=col+1; row<numRows; row++) {
               sum = mval[row][col];
               for (int k=0; k<col; k++) 
                   sum = sum - decomp[col][k]*decomp[row][k];
               decomp[row][col] = sum/decomp[col][col]; 
           }//for
           
       }//for all columns
       
       decompMatrix = new Array2DRowRealMatrix(decomp);
       
    }//decompose
    
    
    /**
     * Returns the last calculated decomposition matrix.
     * <p>
     * Caution: Every call of this Method will return the same object.
     * Decomposing another matrix will generate a new one.
     */
    public Array2DRowRealMatrix getDecomposition() {
        return decompMatrix;
    }//getDecomposition
    
    
    /**
     * Returns the solution for a linear system with constant vector <code>c</code>. 
     * <p>
     * This method solves a linear system <code>M*x=c</code> for a symmetrical,
     * positiv definit coefficient matrix <code>M</code>. Before using this 
     * method the matrix <code>M</code> must have been decomposed.
     * <p>
     * @throws IllegalStateException    if this methode is called before 
     *                                  a matrix was decomposed
     * @throws IllegalArgumentException if the dimension of <code>c</code> doesn't
     *                                  match the row dimension of <code>M</code>
     */    
    public double[] solve(double[] c) 
    throws IllegalStateException, IllegalArgumentException {
      
        if (decompMatrix == null) {
            throw new IllegalStateException("no decomposed matrix available");
        }//if
        if (decompMatrix.getColumnDimension() != c.length) 
           throw new IllegalArgumentException("matrix dimension mismatch"); 
       
        double[][] decomp = decompMatrix.getData();
        double[] x = new double[decomp.length];
        double sum;
        
        //forward elimination
        for (int i=0; i<x.length; i++) {
            sum = c[i];
            for (int k=0; k<i; k++) 
                sum = sum - decomp[i][k]*x[k];
            x[i] = sum / decomp[i][i];
        }//forward elimination
        
        //backward elimination
        for (int i=x.length-1; i>=0; i--) {
            sum = x[i];
            for (int k=i+1; k<x.length; k++) 
                sum = sum - decomp[k][i]*x[k];        
            x[i] = sum / decomp[i][i];
        }//backward elimination
        
        return x;
    }//solve
    
    
    /**
     * Returns the solution for a linear system with a symmetrical, 
     * positiv definit coefficient matrix <code>M</code> and 
     * constant vector <code>c</code>. 
     * <p>
     * As a side effect, the Cholesky-decomposition <code>L*L^T=M</code> is 
     * calculated and internally stored.
     * <p>
     * This is a convenience method for <code><pre>
     *   solver.decompose(m);
     *   solver.solve(c);
     * </pre></code>
     * @throws IllegalArgumentException if M ist not square, not positive definit
     *                                  or the dimensions of <code>M</code> and 
     *                                  <code>c</code> don't match.                     
     */
    public double[] solve(RealMatrix m, double[] c) 
    throws IllegalArgumentException {
        decompose(m);
        return solve(c);
    }//solve 
    
    
    /**
     * Returns the determinant of the a matrix <code>M</code>.
     * <p>
     * Before using this  method the matrix <code>M</code> must 
     * have been decomposed. 
     * <p>
     * @throws IllegalStateException  if this method is called before 
     *                                a matrix was decomposed
     */
    public double getDeterminant() {
        
        if (decompMatrix == null) {
            throw new IllegalStateException("no decomposed matrix available");
        }//if
        
        double[][] data = decompMatrix.getData();
        double res = 1.0d; 
        for (int i=0; i<data.length; i++) {
            res *= data[i][i];
        }//for
        res = res*res;
        
        return res;
    }//getDeterminant
    
}//class CholeskySolver