Package stallone.algebra

Source Code of stallone.algebra.RealCholeskyDecomposition

/*
* To change this template, choose Tools | Templates
* and open the template in the editor.
*/
package stallone.algebra;

import java.security.InvalidParameterException;

import stallone.api.algebra.*;
import stallone.api.doubles.Doubles;
import stallone.api.doubles.IDoubleArray;


/**
* Cholesky decomposition.<br/>
*
* <p>For an n-by-n matrix A, the Cholesky decomposition is the n-by-n
* upper triangular matrix R so that</p>
*
* <pre>
* A = R^T*R is satisfied.
* </pre>
*
* <p>The Cholesky decomposition exists, iff the input matrix A is positive
* definite. Otherwise the algorithm will fail, so that the Cholesky
* decomposition can be used as a test for positive definiteness of A using
* {@link #hasDecomposition()} .
* </p>
*
* @author Axel Rack
*/
public class RealCholeskyDecomposition implements ICholeskyDecomposition {

    /**
     * The input matrix A for the decomposition.
     */
    private IDoubleArray matrixA;
    /**
     * The result matrix R of the decomposition.
     */
    private IDoubleArray matrixR;

    /**
     * Create a new instance of this which tries to compute the decomposition
     * matrix L from A. If the input A is not positive definite, the
     * decomposition will fail and L will be {@code null}.
     *
     * <p>
     * If input A is not purely real or non-symmetric, an unchecked
     * {@see InvalidParameterException} is thrown.
     * </p>
     *
     * @param input The input matrix, which must be purely real and symmetric.
     */
    public RealCholeskyDecomposition(IDoubleArray input) {
        setInputMatrix(input);
        setResultMatrix(computeDecomposition(input));
    }

    @Override
    public boolean hasDecomposition() {
        return getRMatrix() != null;
    }

    /**
     * Get the decomposition matrix L if it exists.
     *
     * @return The output matrix L or {@code null}, if the algorithm failed.
     */
    @Override
    public IDoubleArray getRMatrix() {
        return matrixR;
    }

    @Override
    public IDoubleArray getInputMatrix() {
        return matrixA;
    }

    /**
     * Set the input matrix A to given matrix.
     *
     * @param matrix The new input matrix.
     */
    protected void setInputMatrix(IDoubleArray matrix) {
        if (matrix == null) {
            throw new InvalidParameterException("Input matrix must not be null");
        }
        if (!Algebra.util.isSymmetric(matrix)) {
            throw new InvalidParameterException("Input matrix must be symmetric");
        }
        this.matrixA = matrix;
    }

    /**
     * Set the decomposition result matrix R to given {@code IMatrix}.
     *
     * @param matrix The new result.
     */
    protected void setResultMatrix(IDoubleArray matrix) {
        this.matrixR = matrix;
    }

    /**
     * Computes the Cholesky decomposition of given input matrix A, i.e.
     * a matrix L satisying A=R^T*R for positive definite A.<br>
     * If A is not positive definite, the algorithm fails.
     *
     * @param input The {@code IMatrix} to be decomposed.
     * @return The upper triangular decomposition matrix R, or {@code null} if
     * the algorithm fails.
     */
    private IDoubleArray computeDecomposition(IDoubleArray input) {
        int n = input.rows();
        IDoubleArray result = Doubles.create.array(n, n);
        result.zero();

        for (int i = 0; i < n; i++) {
            for (int j = 0; j <= i; j++) {
                double sum = 0.0;
                for (int k = 0; k < j; k++) {
                    sum += result.get(i, k) * result.get(j, k);
                }
                if (i == j) {
                    result.set(i, i, Math.sqrt(input.get(i, i) - sum));
                } else {
                    result.set(i, j, 1 / result.get(j, j) *
                            (input.get(i, j) - sum));
                }
            }
            if (result.get(i, i) <= 0) {
                result = null;
                break;
            }
        }

        return result;
    }
}
TOP

Related Classes of stallone.algebra.RealCholeskyDecomposition

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.