bugGNU Octave - Bugs: bug #36437, chol2inv can't handle sparse matrix

 
 

bug #36437: chol2inv can't handle sparse matrix

Submitter:  Nir Krakauer <nir_krakauer>
Submitted:  Wed 09 May 2012 01:02:15 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Nir Krakauer Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Thu 11 Aug 2016 05:47:13 PM UTC, comment #5: 

FYI - tests that rely on CHOLMOD should be marked with "testif HAVE_CHOLMOD", applied that change here

http://hg.savannah.gnu.org/hgweb/octave/rev/2c7acd39e9b7

Mike Miller <mtmiller>
Group Member
Tue 09 Aug 2016 11:10:03 PM UTC, comment #4: 

I have pushed the fix.  Some notes of this fix:

  • sparse matrices that were previously identified as banded, but were banded AND lower/upper, are now upper or lower and not banded.
  • I noticed that non-sparse matrices do not check for banded matrices.


http://hg.savannah.gnu.org/hgweb/octave/rev/5417dad26a25

http://hg.savannah.gnu.org/hgweb/octave/rev/71d86e88589f

Note, I had to adjust the precision in one of your tests.

Carnë Draug <carandraug>
Group Member
Sat 06 Aug 2016 05:48:35 PM UTC, comment #3: 

I added some tests, you can look at it at my repository:
https://bitbucket.org/BarbaraLocsi/octave/commits/branch/barbara-chol2inv

Barbara Lócsi <bumi>
Tue 26 Jul 2016 12:53:33 PM UTC, comment #2: 

I couldn't find any tests for chol2inv.  Could you write some before I push this?

Carnë Draug <carandraug>
Group Member
Sat 12 Mar 2016 10:03:35 PM UTC, comment #1: 

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)

Barbara Lócsi <bumi>
Wed 09 May 2012 01:02:15 AM UTC, original submission:  


chol2inv seems not to work when the input is a sparse upper triangular matrix. For example:

n = 10;
X = spdiags(ones(n, 1) * [1 2 1], [-1 0 1], n, n);
Y = chol(X);
Z = chol2inv(Y)
error: chol2inv: spchol2inv requires triangular matrix

On the other hand,
Z = chol2inv(full(Y))
returns the expected result

Nir Krakauer <nir_krakauer>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #36613:  bug_36437.diff added by bumi (2KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by bumi (Updated the item)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by nir_krakauer (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 7 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-08-09 carandraug StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2016-03-12 bumi Attached File- Added bug_36437.diff, #36613
    2013-10-28 mtmiller CategoryNone Octave Function
        Item GroupNone Incorrect Result
        StatusNone Confirmed
        Release3.6.1 dev

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code