Getting a more functional Style for Cholesky

[mikera]

Hi @prasant94 - The initial Cholesky stuff looks good, but I think we should try to convert this to a more "functional" interface.

By which I mean:

• There should be no public mutator functions (e.g. setExpectedMaxSize) - any such values should be set as final fields when objects are constructed. It's OK to mutate as an implementation detail of course
• Immutable result objects should be returned from decompositions
• There is a pure functional interface to make the decomposition happen.
• API users should never need to use anything from the "*.impl" packages - this should be purely for implementation details

As an example, I've converted the old Cholesky code to return a CholeskyResult in this style

[mikera]

Note: This should hopefully make the Cholesky decomposition easier to understand and use for users.... stateful interfaces are always trickier to use because it isn't obvious what order you need to set values etc.

[prasant94]

If the result is returned as a CholeskyResult, why not remove the get methods from Cholesky completely?
How about having the decompose method in ICholesky interface but not any of the get methods, and the get methods in CholeskyResult?

[prasant94]

I mean,

• ICholesky interface that has a decompose method, and is implemented by CholeskyCommon and CholeskyLDU.
• No get methods in CholeskyCommon and CholeskyLDU.
• CholeskyResult class that has getL, getD and getU, and is the return type for the decompose methods.

or maybe a CholeskyLDUResult that extends CholeskyResult. Then we can add getD to this and remove it form CholeskyResult

[prasant94]

I also removed the setExpectedMaxSize function from CholeskyLDU. Now it just assigns new arrays when the decompose method is called.
Take a look, if this is ok, I'll do this for the other classes as well.

public class CholeskyLDU implements ICholeskyLDU {

    // it can decompose a matrix up to this width
    private int maxWidth;
    // width and height of the matrix
    private int n;

    // the decomposed matrix
    private Matrix L;
    private double[] el;

    // the D vector
    private double[] d;

    // tempoary variable used by various functions
    double vv[];

    public CholeskyResult decompose( AMatrix mat ) {
        if( mat.rowCount() != mat.columnCount() ) {
            throw new RuntimeException("Can only decompose square matrices");
        }
        n = mat.rowCount();
        this.vv = new double[n];
        this.d = new double[n];
        L = mat.toMatrix();
        this.el = L.data;

        double d_inv=0;
        for( int i = 0; i < n; i++ ) {
            for( int j = i; j < n; j++ ) {
                double sum = el[i*n+j];

                for( int k = 0; k < i; k++ ) {
                    sum -= el[i*n+k]*el[j*n+k]*d[k];
                }

                if( i == j ) {
                    // is it positive-definate?
                    if( sum <= 0.0 )
                        return null;

                    d[i] = sum;
                    d_inv = 1.0/sum;
                    el[i*n+i] = 1;
                } else {
                    el[j*n+i] = sum*d_inv;
                }
            }
        }
        // zero the top right corner.
        for( int i = 0; i < n; i++ ) {
            for( int j = i+1; j < n; j++ ) {
                el[i*n+j] = 0.0;
            }
        }

        return new CholeskyResult(L, DiagonalMatrix.create(d), L.getTranspose());
    }

    public ADiagonalMatrix getD() {
        return DiagonalMatrix.create(d);
    }

    public AMatrix getL() {
        return L;
    }

    public double[] _getVV() {
        return vv;
    }

    @Override
    public AMatrix getU() {
        return L.getTranspose();
    }
}

[mikera]

This looks like it is going in the right direction, though I think we can simplify it further.

The object that represents the result (ICholeskyLDU) should ideally not be the same object as the one that implements the decompose methods.

This is an example of the "Single Responsibility Principle" - each object / function should have one and only one reason for existing. So we should ideally separate:

• The object/function that performs the decomposition (decompose)
• The object that represents the result (anything implementing ICholeskyLDU)

Does that make sense? There's some more reading on the SRP here if you are interested, in general I find it to be a very useful technique: http://en.wikipedia.org/wiki/Single_responsibility_principle

Originally the SRP applied to mutable objects, but it also works well in functional programming - each entity should have one function / reason for existing.