Hipparchus-Math / hipparchus

An efficient, general-purpose mathematics components library in the Java programming language
Apache License 2.0
142 stars 41 forks source link

SemiDefinitePositiveCholeskyDecomposition wrong output? #255

Open JMmontilla opened 1 year ago

JMmontilla commented 1 year ago

I was comparing the decomposition of sigma points of MerweUnscentedTransform.java with my own implementation using Matlab, and the sigma points 7 and 13 were exactly the same for the Hipparchus implementation. I have tracked it down to the cholesky descomposition in the unscentedTransform method of AbstractUnscentedTransform.java, where the last column is all zeros.

I leave here the example of matrix that I am decomposing:

example_matrix = equicCov.getMatrix.getData

example_matrix =

   0.002617121548164  -0.000000001368900   0.000000006200114  -0.000000000938743   0.000000001649951  -0.000000046381356
  -0.000000001368900   0.000000000000035   0.000000000000008   0.000000000000003  -0.000000000000012  -0.000000000000023
   0.000000006200114   0.000000000000008   0.000000000000033   0.000000000000007  -0.000000000000003  -0.000000000000120
  -0.000000000938743   0.000000000000003   0.000000000000007   0.000000000000223  -0.000000000000116   0.000000000000256
   0.000000001649951  -0.000000000000012  -0.000000000000003  -0.000000000000116   0.000000000000122  -0.000000000000165
  -0.000000046381356  -0.000000000000023  -0.000000000000120   0.000000000000256  -0.000000000000165   0.000000000001166

And the code in matlab to test Hipparchus decomposition against Matlab's:

% Set tuning parameters
alpha   = 1/sqrt(6);
kappa   = 0;
beta    = 2;

n_x = 6;
% Obtain scaling factors and weights for means
lambda = alpha^2*(n_x+kappa) - n_x;                    
factor = n_x + lambda;

cov                     = equicCov.getMatrix.scalarMultiply(factor);
descp_orekit      = SemiDefinitePositiveCholeskyDecomposition(cov).getL().getData()

descp_matlab    = chol(factor*equicCov.getMatrix.getData(),'lower')

These are my outputs:

descp_orekit =

   0.051157810236213                   0                   0                   0                   0                   0
  -0.000000026758383   0.000000183878808                   0                   0                   0                   0
   0.000000121195853   0.000000058920215   0.000000122227742                   0                   0                   0
  -0.000000018349943   0.000000012429491   0.000000065605215   0.000000467137920                   0                   0
   0.000000032252183  -0.000000062846295  -0.000000025892626  -0.000000241376570   0.000000241208261                   0
  -0.000000906632940  -0.000000259075206   0.000000042993360   0.000000512286272  -0.000000111560863                   0

descp_matlab =

   0.051157810236213                   0                   0                   0                   0                   0
  -0.000000026758383   0.000000183878808                   0                   0                   0                   0
   0.000000121195853   0.000000058920215   0.000000122227742                   0                   0                   0
  -0.000000018349943   0.000000012429491   0.000000065605215   0.000000467137920                   0                   0
   0.000000032252183  -0.000000062846295  -0.000000025892626  -0.000000241376570   0.000000241208261                   0
  -0.000000906632940  -0.000000259075206   0.000000042993360   0.000000512286272  -0.000000111560863   0.000000005912609

I can work with my implementation for the UT, but I open this issue just in case there is a bug. It could also be a numerical limitation for this particular matrix.

BryanCazabonne commented 1 year ago

Hi!

I've some difficulties to reproduce your issue. Could you give us the exact value of cov

Thank you, Bryan

BryanCazabonne commented 1 year ago

Indeed, using you current input matrix I have a value for the (6,6) element of the output matrix