This project has moved and is read-only. For the latest updates, please go here.

Decomposition Design

Apr 9, 2010 at 12:15 AM

Just wanted to confirm our design decision for the matrix decomposition classes: let's work with the Cholesky factor class. The first design follows dnA:

  • Cholesky: defines the public API for cholesky decompositions. Depending on the type of the input matrix, it creates a subclass of AbstractCholesky
  • AbstractCholesky: defines all the glue for implementing a cholesky decomposition (keeping track of whether the decomposition was already computed etc)
  • DenseCholesky: subclass of AbstractCholesky which now also routes the actually computation to the correct algorithm (either native or managed)

Another design would be to put all the routing into one Cholesky class. This means merging AbstractCholesky and Cholesky into one class and adding the necessary code to route to the right algorithm inside this class. Pros/Cons?

Another decision we want to make is whether to compute decompositions at constructor time or to do lazy evaluation. I personally prefer that the default is at constructor time as it will be very clear when the computation is done.

Apr 9, 2010 at 5:01 AM

>Another design would be to put all the routing into one Cholesky class. This means merging AbstractCholesky and Cholesky into one class and adding the necessary code to >route to the right algorithm inside this class. Pros/Cons?

I'll think about this weekend - but it seems odd to have a the base class constructor construct one of its sub-classes. One option would be to get rid of the routing, merge UserCholesky and AbstractCholesky into a non-abstract class (Cholesky) and then have DenseCholesky sub-class it. The problem then is if a user does new Cholesky(someDenseMatrix) instead of new DenseCholesky(someDenseMatrix) - less optimized code would be called.

>Another decision we want to make is whether to compute decompositions at constructor time or to do lazy evaluation.

>I personally prefer that the default is at constructor time as it will be very clear when the computation is done.

I'm trying to think of a case where someone would want the decomp to occur on another thread. If they do, it would be cleaner to create the decomp classes and then call Compute.

Apr 9, 2010 at 9:12 AM

Then maybe a hybrid approach? Maybe that matches what you had in mind with your second design, Jurgen.

  • Cholesky with a static method as the primary api (maybe with overloads) to do the basic routing to the correct concrete subclass and then return a decomposition object of type Cholesky; Cholesky to have a protected ctor.
  • DenseCholesky (inheriting Cholesky) or other concrete implementations to route to the correct algorithm; or maybe ManagedRealDenseCholesky plus native versions. Public with public ctor for special scenarios.

The concrete implementations/subclasses could still do the actual computation in a Compute method after construction instead of in the constructor directly (so it could be called separately in special scenarios), but the Cholesky class would not have any such methods, it just holds the results; its static method could automatically call the compute method after constructing one of the subclasses so that in the normal case the user only touches the Cholesky class and doesn't have to call Compute on the result manually.

How would we deal with the second part though: Solve, to actually solve Ax=b or AX=B efficiently based on A's decomposition? Maybe an abstract method in Cholesky to be implemented by the subclasses?

Apr 10, 2010 at 10:56 AM
Edited Apr 10, 2010 at 10:57 AM

OK - I think I see what you are saying. A Cholesky abstract class with static "factory" methods, so something like (based of dnA and only for discussion):

    public abstract class Cholesky
    {
        public static Cholesky Create(Matrix matrix);      

        public bool IsPositiveDefinite();

        public Matrix Factor();

        public double Determinant();

        public abstract Matrix Solve(Matrix input);

        public abstract void Solve(Matrix input, Matrix result);

        public abstract Vector Solve(Vector input);

        public abstract void Solve(Vector input, Vector result);

        public abstract void Compute();

    }

or maybe

    public abstract class Cholesky
   {
public static Cholesky Create(Matrix matrix);

        public bool IsPositiveDefinite();

        public Matrix Factor();

        public double Determinant();

        public Matrix Solve(Matrix input);

        public void Solve(Matrix input, Matrix result);

        public Vector Solve(Vector input);

        public void Solve(Vector input, Vector result);

        protected abstract void Compute();

        protected abstract void DoSolve(Matrix factor, Matrix result);

        protected abstract void DoSolve(Matrix factor, Vector result);

    }

Then we'd have internal DenseCholesky and SparseCholesky sub-classes. We don't need to have versions for managed or native, since the one class would call the LAPACK provider (that would be native or managed).

 

 

 

 

 

 

Apr 10, 2010 at 11:41 PM

So the use case would be

var A = new DenseMatrix(10,10);
var x = new DenseMatrix(10);
// code to add stuff to A

var chol = Cholesky.Create(A);

var y = chol.Solve(x);
var det = chol.Determinant();

I kind of like that. The only weird thing is that for most other things in Math.Net Numerics, you always create an object using a constructor and it might be weird (from the user's perspective) why we suddenly chose to go for a static factory method?

J

 

Apr 11, 2010 at 12:23 AM
Edited Apr 11, 2010 at 12:30 AM

Exactly (both posts).

We could even add extension methods like public static Cholesky Cholesky(this Matrix a) { return Cholesky.Create(a); } for further simplification and discovery, i.e.

var chol = A.Cholesky();   or maybe: var chol = A.CholeskyDecompose();

If it makes sense from the consistency point of view we could actually make DenseCholesky and SparseCholesky public as well, with a public constructor. After all, DenseMatrix and SparseMatrix are public. In the end a tradeoff between consistency and discoverability/simplicity. Not sure about that yet.


Slightly different approach, just for discussion:

  • less weight on the Cholesky base class (like you'd often work directly with DenseMatrix in code, instead of Matrix)
  • public DenseCholesky and SparseCholesky class with public constructors (taking DenseMatrix or SparseMatrix respectively)
  • extension method Cholesky, overloaded with specific return types, i.e. DenseMatrix -> DenseCholesky, SparseMatrix -> SparseCholesky
  • additional extension method overload for the base classes: Matrix -> Cholesky, simply calling Cholesky.Create as mentioned above
Apr 11, 2010 at 4:21 AM

Sounds good. But there is one case I forgot, UserCholesky - when users extend Matrix. I guess they would just use the factory method and it will return UserCholeskty if the Matrix is not a DenseMatrix or SparseMatrix.

 

May 13, 2010 at 7:37 PM
OK, I'm going to take a stab at this. Last question: shall I put this under MathNet.Numerics.LinearAlgebra.Double.Decompositions? J
Apr 5, 2011 at 5:51 PM

Hi, I am just migrating my code from MathNet Iridium to MathNet Numerics, and wanted to confirm that based on this post, this piece of code

Matrix C = Matrix.Transpose(S.CholeskyDecomposition.TriangularFactor);  

should become

var C = Cholesky.Create(S).Factor.Transpose();

Thanks !

Apr 5, 2011 at 6:22 PM

Hi dondestasmanu,

Yes, that sounds about right. Or you could use the extension methods:

var C = S.Cholesky().Factor.Transpose();

Thanks,
Chris

Apr 5, 2011 at 6:26 PM

Do I need to do something to "see" the extension methods ? maybe my dll version is outdated, but i believe I have the last stable release (v4.0.30319).

Thanks!

Apr 5, 2011 at 6:47 PM

Good point, it seems that you currently need to include the following namespace (we may want to change that):

MathNet.Numerics.LinearAlgebra.Generic.Factorization

Thanks,
Chris

Apr 5, 2011 at 7:01 PM

you're quick ! Thanks a lot for your help Chris !!

Jun 6, 2011 at 1:05 PM

Hi,

I try to do some Cholesky but I can't do the final step.

So I have a Matrix MKorr and a Vector U.

I want to do a Cholesky decomposition of MKorr and mult with U....

I did

var Chol = MKorr.Cholesky();

 

and I have the correct results (to tst I have 2*2 Matrix and can check it) but I can't mult Chol*U.

 Any idea why?

thanks

 

Jun 6, 2011 at 2:03 PM

Hi,

>var Chol = MKorr.Cholesky();

Chol is a Cholesky decomposition instance, not the results of the decomposition. You'll need to use Chol.Factor or you can do:

var Chol = MKorr.Cholesky().Factor;

Regards,
Marcus


Jun 8, 2011 at 2:40 PM

Thanks cuda.

This works but I face another prob.

I work with Var/Cov matrix and calculate this matrix out of a Correlation matrix.

This Correlation matrix can have negativ correlations but in this case I can't do Cholesky.

The reason is that the matrix is not positiv definit. Any workaround/solution for this?

 

Thanks,

TTR

Jun 8, 2011 at 4:02 PM

Cholesky just requires that the matrix be positive definite, not that all of the matrix's values be positive.  What is the exact error message you are getting? Can you post or e-mail me the matrix you are trying do decompose?