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

[Anders P. <an...@op...>]

"java.lang.RuntimeException: Matrix is not symmetric positive definite" is a message to you - not (necessarily) a bug.

You've chosen to use the JAMA Cholesky decomposition (part of ojAlgo). JAMA checks for symmetry using A[k][j] == A[j][k] and if the matrix does not pass as being symmetric you will get that exception if you call solve(). I have no plans to change JAMA's behavior in this respect.

I haven't actually looked at the code you sent me yet, but I'm sure this is the problem. Switch to using ojAlgo's own/native Cholesky decomposition and the problem should go away.

/Anders


On 4 apr 2011, at 05.52, Chris Forden wrote:

> 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
>>> 
>>> 
> 
> <Inverted covar matrix that fails oj Cholesky compute.txt><CholeskyBugDemo.zip>------------------------------------------------------------------------------
> Create and publish websites with WebMatrix
> Use the most popular FREE web apps or write code yourself; 
> WebMatrix provides all the features you need to develop and 
> publish your website. http://p.sf.net/sfu/ms-webmatrix-sf
> _______________________________________________
> ojAlgo-user mailing list
> ojA...@li...
> https://lists.sourceforge.net/lists/listinfo/ojalgo-user