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

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

Hi Anders,

Excellent.  I think this gives us the info we need.  I hope to try this
out within the week.

Thanks so much,
-Chris


On Tue, 2011-04-05 at 20:55 +0200, Anders Peterson wrote:
> The class CholeskyDecomposition has three static factory methods: makeBig(), makeComplex(), makePrimitive() and makeJama(). You should use one of those, and not worry about which subclass is actually returned.
> 
> The pattern is the same for the other matrix decompositions: EigenvalueDecomposition, BidiagonalDecomposition, HessenbergDecomposition, LUDecomposition, QRDecomposition, SchurDecomposition, TridiagonalDecomposition and SingularValueDecomposition.
> 
> Don't think there is anything wrong with your matrices. It's just that JAMA wont allow even the smallest difference between transposed elements - they have to be EXACTLY the same. The "native" ojAlgo Cholesky decomposition wont check for symmetry unless you specify that, and when/if it does check it allows for small differences. I guess Parallel  does something similar.
> 
> /Anders
> 
> 
> On 5 apr 2011, at 07.49, Chris Forden wrote:
> 
> > 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
> > 
> > 
> > 
> > ------------------------------------------------------------------------------
> > 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
> 
> 
> ------------------------------------------------------------------------------
> 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