[ojAlgo-user] Files to demo Cholesky failure in versions 28 and 29.45

[Chris F. <cf...@co...>]

Hi Anders,

I eventually did find a case in which ojAlgo's latest snapshot (29.45)
failed to calculate an inverse, although it successfully decomposed the
input matrix.  I have attached a small Eclipse project in zipped format,
which I created by extracting 3 files from our project.  I added some
statements in CovarianceWeighter.calcInverse(JamaMatrix), declaring
convenient local vars for debugging.  The code requires ojAlgo 29.45,
and if you want for comparison, Parallel Colt 0.9.4.

You can reproduce the problem by debugging
TestCovarianceWeighter.testSelectRefStations() as a Junit test; an
exception gets thrown:

java.lang.RuntimeException: Matrix is not symmetric positive definite.
at
org.ojalgo.matrix.jama.CholeskyDecomposition.solve(CholeskyDecomposition.java:163)
at org.ojalgo.matrix.jama.JamaCholesky.solve(JamaCholesky.java:223)
at
org.ojalgo.matrix.jama.JamaAbstractDecomposition.solve(JamaAbstractDecomposition.java:63)
at org.ojalgo.matrix.jama.JamaCholesky.solve(JamaCholesky.java:1)
at org.ojalgo.matrix.jama.JamaCholesky.getInverse(JamaCholesky.java:90)
at
gov.usgs.dismodel.calc.inversion.CovarianceWeighter.calcInverse(CovarianceWeighter.java:212)
at
gov.usgs.dismodel.calc.inversion.CovarianceWeighter.calcRefSubtractedInverse(CovarianceWeighter.java:424)
at
gov.usgs.dismodel.calc.inversion.CovarianceWeighter.selectReferenceStation(CovarianceWeighter.java:401)
at
gov.usgs.dismodel.calc.inversion.TestCovarianceWeighter.testSelectRefStations(TestCovarianceWeighter.java:275)

Parallel Colt performed the decomposition and the inversion without
throwing any exceptions.  OjAlgo (since v28) succeeds in decomposing the
input matrix into an L matrix, but throws a non SPD exception when
getInverse() is then called.  I suspect that the problem might be a
tolerance for SPDness; if so, some public API for adjusting the
tolerance, might be all that is required and should be mentioned in the
Javadocs for JamaCholesky.

The problem described above did not force us to replace ojAlgo with P.
Colt, however.  Instead I merely used a general matrix inversion,
supplied by ojAlgo, in CovarianceWeighter.calcInverse(JamaMatrix).
However, a different matrix, in the attached file “Inverted covar matrix
that fails oj Cholesky compute.txt”, caused JamaCholesky.compute() to
return false, and did not yield L and R matrices that, when multiplied,
reconstituted the original input matrix.  For that calculation, we went
to Parallel Colt which succeeds in decomposing the input.  I have tried
ojAlgo v28 and v29.45 with this matrix file.  I have not ruled out the
possibility that in either of these submitted examples, that the inputs
were really not SPD, although they looked symmetric and can be
decomposed and reconstituted (L*L') by P. Colt.

Below is some code I used on “Inverted covar matrix that fails oj
Cholesky compute.txt”, adapted from
http://code.google.com/p/dismodel/source/browse/ehz_gps/trunk/modeling/dismodel/src/main/java/gov/usgs/dismodel/calc/inversion/CovarianceWeighter.java


	/**
	 * Calculate what our MATLAB scripts call the W matrix, for solving our
systems
	 * of linear equations.  This weights equations by the inverse of the
data
	 * covariance matrix.
	 */
	protected void calcWeighterMatrix(){
	    try {
			JamaCholesky cholTemp = new JamaCholesky();
			if (cholTemp.compute(invcov)){
			    weighterMatrix = cholTemp.getR();
			}
			else {
				JamaMatrix R  = cholTemp.getR();
				JamaMatrix L  = cholTemp.getL();
				JamaMatrix RT  =R.transpose();
				JamaMatrix product =  L.multiplyRight((BasicMatrix)R);

		        /* Parallel Colt's Cholesky Decomposition: */
		        DenseDoubleMatrix2D in = new
DenseDoubleMatrix2D(ArrayUtils.toRawCopyOf(invcov));
		        DenseDoubleCholeskyDecomposition coltCho = new 
		                DenseDoubleCholeskyDecomposition(in);
		        DoubleMatrix2D cL = coltCho.getL();
		        DoubleMatrix2D cLt = coltCho.getLtranspose();
		        DoubleMatrix2D cReconstituted = new
DenseDoubleMatrix2D(in.rows(), in.rows());
		        cL.zMult(cLt, cReconstituted);

				int dummy  = 4; // Set breakpoint here;
			}
		} catch (Exception e) {
			e.printStackTrace();
		}
	}
	

Thanks,
--Chris

On Fri, 2011-04-01 at 09:31 +0200, Anders Peterson wrote: 
> Did you encounter any problems switching from v29 to v30?
> 
> Is the bug still there or not?
> 
> /Anders
> 
> 
> On 28 mar 2011, at 22.25, Chris Forden wrote:
> 
> > Hi,
> > 
> > Thanks!  I will try the recent snapshot.
> > 
> > --Chris
> > 
> >