ojalgo-user

[ojAlgo-user] Cholesky decomp broken for many matrices in versions 28 and 29

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

Hi,

Thanks for the wonderful software which our the main math workhorse lib
in http://code.google.com/p/dismodel/  !

We found that there was a bad bug in your Cholesky decomposition
whenever we used it to try to solve some matrices which had been
processed by what we call differencing out our reference-station data
form our covariance matrices.  I understand that probably is not
meaningful to you, but we solved our problem quickly by adding all the
Parallel Colt libs to our project just for Cholesky decomp.  That was an
inelegant solution but worked well.  The problem occurred in modest
sized matrices (such as 18 x 18).

I wonder if your recent (version 30, unreleased) work on your Cholesky
algorithm is likely to fix that bug.  Is it necessary for us to write a
bug report?  The current status of our project makes it not very easy to
extract example matrices to send you for analysis and testing, although
I might be able to do so if necessary.

Thanks so much again for ojAlgo! 

--Chris Forden

Re: [ojAlgo-user] Cholesky decomp broken for many matrices in versions 28 and 29

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

Sounds like you don't use any of the snapshot releases. There has been at least one bug fix related to solving equation systems and inverting matrices using the Cholesky and LU decompositions. Please, try the latest snapshot release. If the bug is still there I will need some sort of test case that demonstrates the problem.

It's been very long since v29 was released. Too long. A lot of things have changed with v30. The change log is probably not complete...

/Anders


On 28 mar 2011, at 18.44, Chris Forden wrote:

> Hi,
> 
> Thanks for the wonderful software which our the main math workhorse lib
> in http://code.google.com/p/dismodel/  !
> 
> We found that there was a bad bug in your Cholesky decomposition
> whenever we used it to try to solve some matrices which had been
> processed by what we call differencing out our reference-station data
> form our covariance matrices.  I understand that probably is not
> meaningful to you, but we solved our problem quickly by adding all the
> Parallel Colt libs to our project just for Cholesky decomp.  That was an
> inelegant solution but worked well.  The problem occurred in modest
> sized matrices (such as 18 x 18).
> 
> I wonder if your recent (version 30, unreleased) work on your Cholesky
> algorithm is likely to fix that bug.  Is it necessary for us to write a
> bug report?  The current status of our project makes it not very easy to
> extract example matrices to send you for analysis and testing, although
> I might be able to do so if necessary.
> 
> Thanks so much again for ojAlgo! 
> 
> --Chris Forden
> 
> 
> ------------------------------------------------------------------------------
> 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
> 
> 

Re: [ojAlgo-user] Cholesky decomp broken for many matrices in versions 28 and 29

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

Hi,

Thanks!  I will try the recent snapshot.

--Chris



On Mon, 2011-03-28 at 22:21 +0200, Anders Peterson wrote:
> Sounds like you don't use any of the snapshot releases. There has been at least one bug fix related to solving equation systems and inverting matrices using the Cholesky and LU decompositions. Please, try the latest snapshot release. If the bug is still there I will need some sort of test case that demonstrates the problem.
> 
> It's been very long since v29 was released. Too long. A lot of things have changed with v30. The change log is probably not complete...
> 
> /Anders
> 
> 
> On 28 mar 2011, at 18.44, Chris Forden wrote:
> 
> > Hi,
> > 
> > Thanks for the wonderful software which our the main math workhorse lib
> > in http://code.google.com/p/dismodel/  !
> > 
> > We found that there was a bad bug in your Cholesky decomposition
> > whenever we used it to try to solve some matrices which had been
> > processed by what we call differencing out our reference-station data
> > form our covariance matrices.  I understand that probably is not
> > meaningful to you, but we solved our problem quickly by adding all the
> > Parallel Colt libs to our project just for Cholesky decomp.  That was an
> > inelegant solution but worked well.  The problem occurred in modest
> > sized matrices (such as 18 x 18).
> > 
> > I wonder if your recent (version 30, unreleased) work on your Cholesky
> > algorithm is likely to fix that bug.  Is it necessary for us to write a
> > bug report?  The current status of our project makes it not very easy to
> > extract example matrices to send you for analysis and testing, although
> > I might be able to do so if necessary.
> > 
> > Thanks so much again for ojAlgo! 
> > 
> > --Chris Forden
> > 
> > 
> > ------------------------------------------------------------------------------
> > 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
> > 
> > 
> 
> 
> ------------------------------------------------------------------------------
> 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

Re: [ojAlgo-user] Cholesky decomp broken for many matrices in versions 28 and 29

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

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
> 
> 
> 
> On Mon, 2011-03-28 at 22:21 +0200, Anders Peterson wrote:
>> Sounds like you don't use any of the snapshot releases. There has been at least one bug fix related to solving equation systems and inverting matrices using the Cholesky and LU decompositions. Please, try the latest snapshot release. If the bug is still there I will need some sort of test case that demonstrates the problem.
>> 
>> It's been very long since v29 was released. Too long. A lot of things have changed with v30. The change log is probably not complete...
>> 
>> /Anders
>> 
>> 
>> On 28 mar 2011, at 18.44, Chris Forden wrote:
>> 
>>> Hi,
>>> 
>>> Thanks for the wonderful software which our the main math workhorse lib
>>> in http://code.google.com/p/dismodel/  !
>>> 
>>> We found that there was a bad bug in your Cholesky decomposition
>>> whenever we used it to try to solve some matrices which had been
>>> processed by what we call differencing out our reference-station data
>>> form our covariance matrices.  I understand that probably is not
>>> meaningful to you, but we solved our problem quickly by adding all the
>>> Parallel Colt libs to our project just for Cholesky decomp.  That was an
>>> inelegant solution but worked well.  The problem occurred in modest
>>> sized matrices (such as 18 x 18).
>>> 
>>> I wonder if your recent (version 30, unreleased) work on your Cholesky
>>> algorithm is likely to fix that bug.  Is it necessary for us to write a
>>> bug report?  The current status of our project makes it not very easy to
>>> extract example matrices to send you for analysis and testing, although
>>> I might be able to do so if necessary.
>>> 
>>> Thanks so much again for ojAlgo! 
>>> 
>>> --Chris Forden
>>> 
>>> 
>>> ------------------------------------------------------------------------------
>>> 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
>>> 
>>> 
>> 
>> 
>> ------------------------------------------------------------------------------
>> 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
> 
> 
> 
> ------------------------------------------------------------------------------
> 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
> 
> 

Re: [ojAlgo-user] Cholesky decomp broken for many matrices in versions 28 and 29

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

Hi Anders,

I will try testing v30 in the next day or two.  Meanwhile a colleague has
reported other problems, unrelated to Cholesky decomp, but I do not know the
details.

--Chris

-----Original Message-----
From: Anders Peterson [mailto:an...@op...] 
Sent: Friday, April 01, 2011 12:32 AM
To: ojAlgo ojAlgo
Subject: Re: [ojAlgo-user] Cholesky decomp broken for many matrices in
versions 28 and 29

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
> 
> 
> 
> On Mon, 2011-03-28 at 22:21 +0200, Anders Peterson wrote:
>> Sounds like you don't use any of the snapshot releases. There has been at
least one bug fix related to solving equation systems and inverting matrices
using the Cholesky and LU decompositions. Please, try the latest snapshot
release. If the bug is still there I will need some sort of test case that
demonstrates the problem.
>> 
>> It's been very long since v29 was released. Too long. A lot of things
have changed with v30. The change log is probably not complete...
>> 
>> /Anders
>> 
>> 
>> On 28 mar 2011, at 18.44, Chris Forden wrote:
>> 
>>> Hi,
>>> 
>>> Thanks for the wonderful software which our the main math workhorse lib
>>> in http://code.google.com/p/dismodel/  !
>>> 
>>> We found that there was a bad bug in your Cholesky decomposition
>>> whenever we used it to try to solve some matrices which had been
>>> processed by what we call differencing out our reference-station data
>>> form our covariance matrices.  I understand that probably is not
>>> meaningful to you, but we solved our problem quickly by adding all the
>>> Parallel Colt libs to our project just for Cholesky decomp.  That was an
>>> inelegant solution but worked well.  The problem occurred in modest
>>> sized matrices (such as 18 x 18).
>>> 
>>> I wonder if your recent (version 30, unreleased) work on your Cholesky
>>> algorithm is likely to fix that bug.  Is it necessary for us to write a
>>> bug report?  The current status of our project makes it not very easy to
>>> extract example matrices to send you for analysis and testing, although
>>> I might be able to do so if necessary.
>>> 
>>> Thanks so much again for ojAlgo! 
>>> 
>>> --Chris Forden
>>> 
>>> 
>>>
----------------------------------------------------------------------------
--
>>> 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
>>> 
>>> 
>> 
>> 
>>
----------------------------------------------------------------------------
--
>> 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
> 
> 
> 
>
----------------------------------------------------------------------------
--
> 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
> 
> 


----------------------------------------------------------------------------
--
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

Re: [ojAlgo-user] Cholesky decomp broken for many matrices in versions 28 and 29

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

Hi Anders,

For better or worse, I cannot reproduce the bug even in version 28, which is
where I first noticed it.  I have not yet tried 29 or any snapshot of 30
recently.  Perhaps I am not just remembering the steps-to-reproduce
correctly.  Or perhaps the JVM has changed.  I went back to our code just
before swapping out ojAlgo that still has some comments about the problem:

		cho = new JamaCholesky();
		cho.compute(in);
		invcov = null; // in case the inverse (below) throws
		invcov = cho.getInverse();

However, I cannot get the exception to throw (unless in is diagonal with
some zeroes on the diagonal, in which case it should throw).

I have tried 32-bit Linux and 64-bit Vista, both running 32-bit JVMs and
Eclipse.

My colleague recently mentioned he had some troubles, and I now see that the
code he mentioned also uses Cholesky matrices.  I will ask him if he still
has the trouble.

Sorry for the bother (but in September, there really was a problem).  :)

--Chris


-----Original Message-----
From: Anders Peterson [mailto:an...@op...] 
Sent: Friday, April 01, 2011 12:32 AM
To: ojAlgo ojAlgo
Subject: Re: [ojAlgo-user] Cholesky decomp broken for many matrices in
versions 28 and 29

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
> 
> 
> 
> On Mon, 2011-03-28 at 22:21 +0200, Anders Peterson wrote:
>> Sounds like you don't use any of the snapshot releases. There has been at
least one bug fix related to solving equation systems and inverting matrices
using the Cholesky and LU decompositions. Please, try the latest snapshot
release. If the bug is still there I will need some sort of test case that
demonstrates the problem.
>> 
>> It's been very long since v29 was released. Too long. A lot of things
have changed with v30. The change log is probably not complete...
>> 
>> /Anders
>> 
>> 
>> On 28 mar 2011, at 18.44, Chris Forden wrote:
>> 
>>> Hi,
>>> 
>>> Thanks for the wonderful software which our the main math workhorse lib
>>> in http://code.google.com/p/dismodel/  !
>>> 
>>> We found that there was a bad bug in your Cholesky decomposition
>>> whenever we used it to try to solve some matrices which had been
>>> processed by what we call differencing out our reference-station data
>>> form our covariance matrices.  I understand that probably is not
>>> meaningful to you, but we solved our problem quickly by adding all the
>>> Parallel Colt libs to our project just for Cholesky decomp.  That was an
>>> inelegant solution but worked well.  The problem occurred in modest
>>> sized matrices (such as 18 x 18).
>>> 
>>> I wonder if your recent (version 30, unreleased) work on your Cholesky
>>> algorithm is likely to fix that bug.  Is it necessary for us to write a
>>> bug report?  The current status of our project makes it not very easy to
>>> extract example matrices to send you for analysis and testing, although
>>> I might be able to do so if necessary.
>>> 
>>> Thanks so much again for ojAlgo! 
>>> 
>>> --Chris Forden
>>> 
>>> 
>>>
----------------------------------------------------------------------------
--
>>> 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
>>> 
>>> 
>> 
>> 
>>
----------------------------------------------------------------------------
--
>> 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
> 
> 
> 
>
----------------------------------------------------------------------------
--
> 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
> 
> 


----------------------------------------------------------------------------
--
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

[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
> > 
> > 

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

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

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

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