[Octave-bug-tracker] [bug #36437] chol2inv can't handle sparse matrix

[Barbara Lócsi]

Follow-up Comment #1, bug #36437 (project octave):

I can see it is still not working, (only the error is bit different)


error: chol2inv: U must be a triangular matrix


The chol2inv function in sparse-chol.cc:526, requires the type of the matrix
to be determined. If it is an upper or lower matrix than it determines the
retval, but otherwise it gives this error.

This Y matrix is an upper matrix so i think the type of the matrix is
determined incorrectly.
MatrixType.cc lines 340-364. If it is a banded matrix then the function won't
check whether it is lower or upper matrix.
But in this case, this matrix (Y) is both banded and upper.


   1.41421   0.70711   0.00000   0.00000   0.00000   0.00000   0.00000  
0.00000   0.00000   0.00000
   0.00000   1.22474   0.81650   0.00000   0.00000   0.00000   0.00000  
0.00000   0.00000   0.00000
   0.00000   0.00000   1.15470   0.86603   0.00000   0.00000   0.00000  
0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   1.11803   0.89443   0.00000   0.00000  
0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   1.09545   0.91287   0.00000  
0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   1.08012   0.92582  
0.00000   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   1.06904  
0.93541   0.00000   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000  
1.06066   0.94281   0.00000
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000  
0.00000   1.05409   0.94868
   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000  
0.00000   0.00000   1.04881


Temporarily I deleted the else from else if, so even thougth the matrix is
both banded and upper/lower it will say it is upper/lower (before it said it
is banded).
In this case there is no more error from chol2inv but it calculates wrong
result.
The result of 


Z = chol2inv(full(Y)) 


is:


Z =

   0.909091  -0.818182   0.727273  -0.636364   0.545455  -0.454545   0.363636 
-0.272727   0.181818  -0.090909
  -0.818182   1.636364  -1.454545   1.272727  -1.090909   0.909091  -0.727273 
 0.545455  -0.363636   0.181818
   0.727273  -1.454545   2.181818  -1.909091   1.636364  -1.363636   1.090909 
-0.818182   0.545455  -0.272727
  -0.636364   1.272727  -1.909091   2.545455  -2.181818   1.818182  -1.454545 
 1.090909  -0.727273   0.363636
   0.545455  -1.090909   1.636364  -2.181818   2.727273  -2.272727   1.818182 
-1.363636   0.909091  -0.454545
  -0.454545   0.909091  -1.363636   1.818182  -2.272727   2.727273  -2.181818 
 1.636364  -1.090909   0.545455
   0.363636  -0.727273   1.090909  -1.454545   1.818182  -2.181818   2.545455 
-1.909091   1.272727  -0.636364
  -0.272727   0.545455  -0.818182   1.090909  -1.363636   1.636364  -1.909091 
 2.181818  -1.454545   0.727273
   0.181818  -0.363636   0.545455  -0.727273   0.909091  -1.090909   1.272727 
-1.454545   1.636364  -0.818182
  -0.090909   0.181818  -0.272727   0.363636  -0.454545   0.545455  -0.636364 
 0.727273  -0.818182   0.909091


,while the result of


Z = chol2inv(Y);
full(Z)


is:


ans =

   0.500000  -0.288675   0.204124  -0.158114   0.129099  -0.109109   0.094491 
-0.083333   0.074536  -0.067420
  -0.288675   0.833333  -0.589256   0.456435  -0.372678   0.314970  -0.272772 
 0.240563  -0.215166   0.194625
   0.204124  -0.589256   1.166667  -0.903696   0.737865  -0.623610   0.540062 
-0.476290   0.426006  -0.385337
  -0.158114   0.456435  -0.903696   1.500000  -1.224745   1.035098  -0.896421 
 0.790569  -0.707107   0.639602
   0.129099  -0.372678   0.737865  -1.224745   1.833333  -1.549449   1.341863 
-1.183412   1.058475  -0.957427
  -0.109109   0.314970  -0.623610   1.035098  -1.549449   2.166667  -1.876388 
 1.654819  -1.480115   1.338814
   0.094491  -0.272772   0.540062  -0.896421   1.341863  -1.876388   2.500000 
-2.204793   1.972027  -1.783765
  -0.083333   0.240563  -0.476290   0.790569  -1.183412   1.654819  -2.204793 
 2.833333  -2.534210   2.292280
   0.074536  -0.215166   0.426006  -0.707107   1.058475  -1.480115   1.972027 
-2.534210   3.166667  -2.864358
  -0.067420   0.194625  -0.385337   0.639602  -0.957427   1.338814  -1.783765 
 2.292280  -2.864358   3.500000


I think the chol2inv calculates it wrongly.


if (typ == MatrixType::Upper)
    {
      rinv = r.inverse (mattype, info, rcond, true, false);
      retval = rinv.transpose () * rinv;
    }


It is: (inv(Y))' * inv(Y)
While it should be inv (Y'*Y)







(file #36613)
    _______________________________________________________

Additional Item Attachment:

File name: bug_36437.diff                 Size:2 KB


    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?36437>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/