bugGNU Octave - Bugs: bug #32364, problem multiplying permutation...

 
 

bug #32364: problem multiplying permutation and sparse matrices

Submitter:  John W. Eaton <jwe>
Submitted:  Sat 05 Feb 2011 03:02:24 PM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  None Assigned to:  dbateman
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Sun 13 Feb 2011 05:36:51 PM UTC, comment #4: 

Well as you accidentally pushed my patch in

http://hg.savannah.gnu.org/hgweb/octave/rev/afb65a7bc065

I pushed a changelog and am closing this bug

D.

David Bateman <dbateman>
Group Member
Wed 09 Feb 2011 07:02:26 AM UTC, comment #3: 

It seems to work for me.  Would you please add a ChangeLog entry and commit the change?

Thanks.

John W. Eaton <jwe>
Group administrator
Tue 08 Feb 2011 11:53:48 PM UTC, comment #2: 

Strictly speaking the sparse code doesn't need to make the assumption that the row indexes are in order, however doing so makes some code simpler, while making other code more difficult. I made the choice when I started writing the sparse code to assume that the row indexes were ordered and all of the sparse code would have to be gone over to remove this assumption. So probably the best thing to do here is to modify PM SM so that it to orders the row indexes. Fortunately, by its nature SM PM already does so.

Here is a patch that I think does the right thing but its received little testing

D.


(file #22632)

David Bateman <dbateman>
Group Member
Tue 08 Feb 2011 07:57:56 AM UTC, comment #1: 

[Reposting because I screwed up a verbatim tag and the report on the web was incomplete; I wish this thing had a preview button...]

While investigating the following bug report

https://mailman.cae.wisc.edu/pipermail/bug-octave/2011-February/017300.html

I uncovered a problem with multiplying permutation and sparse matrices.

Given


A = sparse ([1   1   1   0   0   0   0   0   0
             0   0   0   1   1   1   0   0   0
             0   0   0   0   0   0   1   1   1
             1   0   0   1   0   0   1   0   0
             0   1   0   0   1   0   0   1   0
             0   0   1   0   0   1   0   0   1])
[L, U, P, Q, R] = lu (A);


I was trying to verify that


P*(R\A)*Q - L*U


was numerically zero, but found


octave:8> P*(R\A)*Q - L*U
ans =

Compressed Column Sparse (rows = 6, cols = 9, nnz = 4 [7.4%])

  (4, 5) -> -0.33333
  (4, 5) ->  0.33333
  (3, 7) -> -0.33333
  (3, 7) ->  0.33333


Note the strange result that some elements apparently have two values!

However,


full (P)*(R\A)*Q - L*U


returns a matrix of zeros.  Playing around a bit more, I noticed


octave:10> P*A
ans =

Compressed Column Sparse (rows = 6, cols = 9, nnz = 18 [33%])

  (1, 1) ->  1
  (3, 1) ->  1
  (1, 2) ->  1
  (4, 2) ->  1
  (1, 3) ->  1
  (6, 3) ->  1
  (2, 4) ->  1
  (3, 4) ->  1
  (2, 5) ->  1
  (4, 5) ->  1
  (2, 6) ->  1
  (6, 6) ->  1
  (5, 7) ->  1
  (3, 7) ->  1
  (5, 8) ->  1
  (4, 8) ->  1
  (5, 9) ->  1
  (6, 9) ->  1


Note the elements in columns 7 and 8 are out of order.  I assume that the sparse matrix multiplication algorithm assumes that the elements in each column are ordered by ascending row number.

I've looked at the code in Sparse-perm-op-defs.h, but I think it will take me a while to figure this out, so I would really appreciate it if someone who understands this code better could take a look at this problem.

I'm adding the following people to the CC list for this report: Jason Reidy since he is listed as the author of Sparse-perm-op-defs.h, Jaroslav Hajek because he is the author of the permutation matrix code, and David Bateman because he has good knowledge of the sparse matrix code.

John W. Eaton <jwe>
Group administrator
Sat 05 Feb 2011 03:02:24 PM UTC, original submission:  

While investigating the following bug report

https://mailman.cae.wisc.edu/pipermail/bug-octave/2011-February/017300.html

I uncovered a problem with multiplying permutation and sparse matrices.

Given


A = sparse ([1   1   1   0   0   0   0   0   0
             0   0   0   1   1   1   0   0   0
             0   0   0   0   0   0   1   1   1
             1   0   0   1   0   0   1   0   0
             0   1   0   0   1   0   0   1   0
             0   0   1   0   0   1   0   0   1])
[L, U, P, Q, R] = lu (A);


I was trying to verify that


P*(R\A)*Q - L*U
-vertabim-

was numerically zero, but found

+verbatim+
octave:8> P*(R\A)*Q - L*U
ans =

Compressed Column Sparse (rows = 6, cols = 9, nnz = 4 [7.4%])

  (4, 5) -> -0.33333
  (4, 5) ->  0.33333
  (3, 7) -> -0.33333
  (3, 7) ->  0.33333


Note the strange result that some elements apparently have two values!

However,


full (P)*(R\A)*Q - L*U


returns a matrix of zeros.  Playing around a bit more, I noticed


octave:10> P*A
ans =

Compressed Column Sparse (rows = 6, cols = 9, nnz = 18 [33%])

  (1, 1) ->  1
  (3, 1) ->  1
  (1, 2) ->  1
  (4, 2) ->  1
  (1, 3) ->  1
  (6, 3) ->  1
  (2, 4) ->  1
  (3, 4) ->  1
  (2, 5) ->  1
  (4, 5) ->  1
  (2, 6) ->  1
  (6, 6) ->  1
  (5, 7) ->  1
  (3, 7) ->  1
  (5, 8) ->  1
  (4, 8) ->  1
  (5, 9) ->  1
  (6, 9) ->  1


Note the elements in columns 7 and 8 are out of order.  I assume that the sparse matrix multiplication algorithm assumes that the elements in each column are ordered by ascending row number.

I've looked at the code in Sparse-perm-op-defs.h, but I think it will take me a while to figure this out, so I would really appreciate it if someone who understands this code better could take a look at this problem.

I'm adding the following people to the CC list for this report: Jason Reidy since he is listed as the author of Sparse-perm-op-defs.h, Jaroslav Hajek because he is the author of the permutation matrix code, and David Bateman because he has good knowledge of the sparse matrix code.


John W. Eaton <jwe>
Group administrator

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #22632:  patch.perm added by dbateman (1KiB - 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 dbateman (Updated the item)
  • -email is unavailable- added by jwe (Submitted the item)
  • -email is unavailable- added by jwe
  • -email is unavailable- added by jwe
  • -email is unavailable- added by jwe
  •  

    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 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2011-02-13 dbateman Assigned toNone dbateman
        Open/ClosedOpen Closed
    2011-02-08 dbateman Attached File- Added patch.perm, #22632
    2011-02-05 jwe Carbon-Copy- Added -email is unavailable-
        Carbon-Copy- Added highegg
        Carbon-Copy- Added dbateman

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code