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

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

Hi Anders,

> Switch to using ojAlgo's own/native Cholesky decomposition and the
problem should go away.

According to http://ojalgo.org/generated/index.html,
CholeskyDecomposition is abstract, and no derived class is listed.  (The
generated Javadocs do not list any "All Known Implementing Classes:"
section for CholeskyDecomposition.  Should I be generating my own
Javadocs from the source, to overcome such missing refs to derived
classes?)



>  if the matrix does not pass as being symmetric you will get that
exception if you call solve(). ... I haven't actually looked at the code
you sent me yet, but I'm sure this is the problem. 

I have added a symmetry test:

private void printSymmetry(JamaMatrix in) {
	JamaMatrix diff = in.subtract(in.transpose());
	double assymNormRatio = diff.getFrobeniusNorm().doubleValue() 
				/ in.getFrobeniusNorm().doubleValue();
	System.out.println(diff.toString());
	System.out.println();
	System.out.println("assymNormRatio= " + assymNormRatio);
}

All the values of diff are small (0 or less than 1e-20), and
assymNormRatio = 5.3710432377270805E-17	
I would think that such small assymetries should be deemed
insignificant, right?

I then tested Positive Definiteness by printing the Eigen values with
this code:

private void printEigens(JamaMatrix in) {
	JamaEigenvalue eigMat = new JamaEigenvalue(); 
	eigMat.compute(in);
	Array1D<ComplexNumber> eigVals = eigMat.getEigenvalues();
	System.out.println();
	System.out.println("eigen Vals= " + eigVals);
	double minEigValR =99999.9;
	for (ComplexNumber c : eigVals) {
		if (minEigValR > c.getReal())
			minEigValR = c.getReal();
	}
	System.out.println("min eigVal real = " + minEigValR);
}


The output was:

eigen Vals= [(4.498282249210615E-4 + 0.0i), (1.5887605890203683E-4 +
0.0i), (1.1332549725716817E-4 + 0.0i), (1.0609982549417765E-4 + 0.0i),
(5.2699055744352376E-5 + 0.0i), (4.8733492649113855E-5 + 0.0i),
(4.2166615429069484E-5 + 0.0i), (3.8203524668310696E-5 + 0.0i),
(3.629597352711602E-5 + 0.0i), (2.6518439788904983E-5 + 0.0i),
(1.5907235051737824E-5 + 0.0i), (1.2258822768511711E-5 + 0.0i),
(1.1744258961933221E-5 + 0.0i), (1.0690845009451796E-5 + 0.0i),
(8.786473150043942E-6 + 0.0i), (5.600551936654057E-6 + 0.0i),
(4.221645596951178E-6 + 0.0i), (4.109779674404523E-6 + 0.0i)]
min eigVal real = 4.109779674404523E-6

The fact that all Eigenvalues have real parts greater than zero, means
the matrix is positive definite, right?

Thanks for your help,
--Chris


On Mon, 2011-04-04 at 22:05 +0200, Anders Peterson wrote:
> "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
> 
> 
> ------------------------------------------------------------------------------
> Xperia(TM) PLAY
> It's a major breakthrough. An authentic gaming
> smartphone on the nation's most reliable network.
> And it wants your games.
> http://p.sf.net/sfu/verizon-sfdev
> _______________________________________________
> ojAlgo-user mailing list
> ojA...@li...
> https://lists.sourceforge.net/lists/listinfo/ojalgo-user