commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From "Brent Worden" <br...@worden.org>
Subject RE: [math] Greetings from a newcomer (but not to numerics)
Date Mon, 26 May 2003 03:21:05 GMT
> I would prefer to implement something that is reasonable for a large
> class of matrices (not just pos definite).  The RealMatrixImpl class
> requires in-memory storage, so we are not going to be working on *large*
> matrices in that implementation.
I have no qualms about that at all.

> Do any of you have implementations that you would like to contribute?

Here is a Cholesky Decomposition class usable by commons-math which can be
used or modified as you see fit:

package org.apache.commons.math;

/**
 * Solve linear systems of equations using Cholesky decomposition.
 *
 * See
http://ikpe1101.ikp.kfa-juelich.de/briefbook_data_analysis/node33.html
 * for further information.
 *
 * @author Brent Worden
 */
public class CholeskyDecomposition {
	/**
	 * Default constructor. Does nothing as this is a utility class.
	 */
	public CholeskyDecomposition(){
		super();
	}

    /**
     * Decompose a symmetric, positive definite matrix into using Cholesky
     * factorization.  Using this technique, mat is decomposed into LL',
where
     * L is a lower triangular matrix.
     * @param mat the matrix to be decomposed.
     * @return the lower triangular matrix of the Cholesky factorization.
     * @throws IllegalArgumentException if mat is not square or is singular.
     */
    public RealMatrix decompose(RealMatrix mat){
        if((mat == null) ||
           (mat.getColumnDimension() <= 0) ||
           (mat.getColumnDimension() != mat.getRowDimension()))
        {
            throw new IllegalArgumentException("matrix must be square.");
        }

        int n = mat.getColumnDimension(); // size of square matrix
        double[][] a = mat.getData();     // original matrix

        // decomposition
        RealMatrix decomposition = new RealMatrixImpl(n, n);
        double[][] data = decomposition.getData();

        data[0][0] = Math.sqrt(a[0][0]);
        for (int i = 1; i < n; i++) {
            data[i][0] = a[i][0] / data[0][0];
        }
        for (int i = 0; i < n; i++) {
            double sum = a[i][i];
            for (int j = 0; j < i; j++) {
                sum = sum - (data[i][j] * data[i][j]);
            }
            if(sum <= 0.0){
                // TODO: should this be a specialized exception?
                throw new IllegalArgumentException("can not decompose
singular matrix.");
            }
            data[i][i] = Math.sqrt(sum);
            for (int j = n - 1; j > i; j--) {
                sum = a[j][i];
                for (int k = 0; k < i; k++) {
                    sum = sum - (data[j][k] * data[i][k]);
                }
                data[j][i] = sum / data[i][i];
            }
        }

        return decomposition;
    }

    /**
     * Take the output of {@link
CholeskyDecomposition#decompose(RealMatrix)}
     * and solves a linear set of equations, LL'x = b.
     * @param mat LL' matrix generated from Cholesky factorization.
     * @param b ??? what is the common name for this vector ???
     * @return the solution vector, x
     * @throws IllegalArgumentException if mat is not square or is not the
same
     *         size as b.
     */
    public double[] solveUsingDecomposition(RealMatrix mat, double[] b){
        if((mat == null) ||
           (mat.getColumnDimension() <= 0) ||
           (mat.getColumnDimension() != mat.getRowDimension()) ||
           (b == null) ||
           (b.length != mat.getColumnDimension()))
        {
            throw new IllegalArgumentException(
                "matrix must be square and the same dimension as b
vector.");
        }

        int n = b.length;
        double[] x = new double[n];
        double[][] data = mat.getData();

        for (int i = 0; i < n; i++) {
		double sum = b[i];
            for (int j = i - 1; j >= 0; j--) {
			sum = sum - (data[i][j] * x[j]);
		}
            x[i] = sum / data[i][i];
    	  }
        for (int i = n - 1; i >= 0; i--) {
            double sum = x[i];
            for (int j = i + 1; j < n; j++) {
                sum = sum - (data[j][i] * x[j]);
            }
            x[i] = sum / data[i][i];
        }

        return x;
    }

    /**
     * Solves a linear set of equations, Ax = b. A must be of a symmetric,
     * positive definite matrix.
     * @param a the symmetric, positive definite matrix
     * @param b ??? what is the common name for this vector ???
     * @return the solution vector, x
     */
    public double[] solve(RealMatrix a, double[] b){
        RealMatrix decomposition = decompose(a);
        return solveUsingDecomposition(decomposition, b);
    }
}

And here is a simple unit test:

package org.apache.commons.math;

import junit.framework.TestCase;

/**
 * @author Brent Worden
 */
public class CholeskyDecompositionTest extends TestCase {

	/**
	 * Constructor for CholeskyDecompositionTest.
	 * @param name
	 */
	public CholeskyDecompositionTest(String name) {
		super(name);
	}

	public void testSolve(){
        double[][] data = {
            { 1,  -1,   1},
            {-1,  10, -10},
            { 1, -10,  14}
        };
        double[] b = {2, -2, 6};
        double[] expected = {2, 1, 1};

        CholeskyDecomposition cd = new CholeskyDecomposition();
        double[] x = cd.solve(new RealMatrixImpl(data), b);

        for (int i = 0; i < x.length; i++) {
            assertEquals(expected[i], x[i], 1e-2);
        }
    }
}

Brent Worden
http://www.brent.worden.org


---------------------------------------------------------------------
To unsubscribe, e-mail: commons-dev-unsubscribe@jakarta.apache.org
For additional commands, e-mail: commons-dev-help@jakarta.apache.org


Mime
View raw message