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

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

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