bugGNU Octave - Bugs: bug #53390, lu factorization of large sparse...

 
 

bug #53390: lu factorization of large sparse matrices often fails

Submitter:  Marco Caliari <caliari>
Submitted:  Tue 20 Mar 2018 03:27:30 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  In Progress Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 09 Apr 2018 12:39:16 AM UTC, comment #20: 

Oops.  Re-opening because the long term issue is not resolved.

Rik <rik5>
Group administrator
Mon 09 Apr 2018 12:38:14 AM UTC, comment #19: 

I checked your changeset in here (http://hg.savannah.gnu.org/hgweb/octave/rev/b8279dd83664).  The only addition I made was to add "local" to the call to warning.  This makes the change only for the enclosing block.  Afterwards, it reverts to whatever value it had before which in this case will be the user's own setting.

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Sun 08 Apr 2018 08:16:20 AM UTC, comment #18: 

I have chosen the warning off alternative. I hoped it was the easiest to revert when lu would have been fixed. But I discovered that something like


warning ("off", "not exsisting anymore id")


does not produce an error, as it does in matlab.

(file #43838)

Marco Caliari <caliari>
Group Member
Sat 07 Apr 2018 08:59:14 PM UTC, comment #17: 

@Marco, @Cristiano: Using 'test bicg' I now get a warning due to this bug report.


test bicg
warning: lu: function may fail when called with less than 4 output arguments and a sparse input
PASSES 13 out of 13 tests


The problem seems to be this BIST test


## Preconditioned technique
%!testif HAVE_UMFPACK
%! A = sparse (toeplitz ([2, 1, 0, 0, 0], [2, -1, 0, 0, 0]));
%! b = sum (A, 2);
%! [M1, M2] = lu (A + eye (5));
%! [x, flag] = bicg (A, b, [], 1, M1, M2);
%! ## b has two columns!
%! [y, flag]  = bicg (M1 \ A / M2, [M1 \ b, M2' \ b], [], 1);
%! assert (x, M2 \ y, 8 * eps);


Could you take a look at how to resolve this?  I have a couple of ideas on how to fix it.  First, if this matrix is simple enough then lu() will succeed.  In that case you could turn off the warning locally.  Alternatively, you could convert the matrix to full, calculate lu, and then use sparse() to turn M1 and M2 to sparse matrices.  As a third possibility, use the 4-output form of lu() and use P and Q appropriately to transform M1 and M2 before usage.

Rik <rik5>
Group administrator
Wed 28 Mar 2018 06:04:13 PM UTC, comment #16: 

I added a new warning ID and message to lu in this cset (http://hg.savannah.gnu.org/hgweb/octave/rev/1c077d652c57).

Rik <rik5>
Group administrator
Wed 28 Mar 2018 04:18:13 PM UTC, comment #15: 

I think that a warning, with the hint to use the 4 outputs version, is enough. In fact, if it works, it works. Otherwise, you get anyway an error.

Marco Caliari <caliari>
Group Member
Wed 28 Mar 2018 03:58:48 PM UTC, comment #14: 

I added just the note to the documentation here (http://hg.savannah.gnu.org/hgweb/octave/rev/fa2f8ffd088e).  Waiting to hear from Marco about whether sparse input with less than four outputs should be an error.

Rik <rik5>
Group administrator
Wed 28 Mar 2018 03:54:54 PM UTC, comment #13: 

Or warning?  It doesn't have to be an error if there is true use case for calling lu on sparse matrices with fewer than 4 outputs.

Rik <rik5>
Group administrator
Wed 28 Mar 2018 03:27:08 PM UTC, comment #12: 

@Marco: I can add that documentation fix, but is a more long term fix just to improve the input validation?  We could issue an error if the input is found to be sparse and < 4 output arguments are requested.  That is very easy to detect.

Rik <rik5>
Group administrator
Wed 28 Mar 2018 06:49:06 AM UTC, comment #11: 

@Rik: this bug has always been in octave and it is due, in my opinion, to a misleading documentation in umfpack (the author told me he will add a clarifying note to the documentation). If it came up only now, it is because people eihter do not use lu of sparse matrices or use it in the right way ([L,U,P,Q] = lu (S), maybe after noting that otherwise it fails). The proposal in comment #9 will take some time and it is not feasible for 4.4.0. What about adding sentence in the documentation? Something like


When called with two or three output arguments and a sparse input
matrix, 'lu' does not attempt to perform sparsity preserving column permutations.  Called with a fourth output argument, the sparsity preserving column transformation Q is returned, such that 'P * A * Q = L * U'. This is the preferred way to call lu with sparse input matrices.


Marco Caliari <caliari>
Group Member
Tue 27 Mar 2018 05:47:01 PM UTC, comment #10: 

@Marco: Does something need to happen on this bug before the 4.4.0 release?

Rik <rik5>
Group administrator
Fri 23 Mar 2018 03:13:07 PM UTC, comment #9: 

I made few tests and Option 3) in comment #8 seems to work (in particular, opts.type="ilutp", opts.droptol=0, opts.thresh=1). I can prepare a file _lu_.cc which is a simplified version of _ilu_.cc in which only the ilutp variant is kept (and forced to continue even if the pivot is zero). How to integrate this into sparse-lu.cc is beyond my possibilities and I need help.

Marco Caliari <caliari>
Group Member
Fri 23 Mar 2018 09:37:26 AM UTC, comment #8: 

I think we have to disable for the moment the calls [L,U]=lu(S) and [L,U,P]=lu(S) for sparse matrices. Even if umfpack does not fail, the result is unpredictable. I do not think it is a big issue: I cannot imagine how to use an LU factorization of a sparse matrix without pivoting and reduction of fill-in ([L,U,P,Q]=lu(S) syntax). Maybe for didactic purposes. Then, there are some options:

1) we implement it from scratch
2) GSOC
3) we try to adapt ilu (as far as I know, ilu is lu with an if clause)

Marco Caliari <caliari>
Group Member
Fri 23 Mar 2018 07:41:10 AM UTC, comment #7: 

It is confirmed that [L,U]=lu(A) or [L,U,P]=lu(A) does not use UMFPACK and UMFPACK cannot work without a minimal column ordering. We have to implement the algorithm from scratch. I will think about it.

Marco Caliari <caliari>
Group Member
Thu 22 Mar 2018 04:33:13 PM UTC, comment #6: 

@Marco: Okay, waiting information from umfpack author.  I've changed the status to "Need Info" until we get an answer.

Rik <rik5>
Group administrator
Thu 22 Mar 2018 04:27:59 PM UTC, comment #5: 

I can reproduce the problem with the mex interface of umfpack in matlab. Therefore, either there is a bug in umfpack or umfpack is not meant to perform an LU factorization without column pivoting (in fact, [L, U, P, Q] = lu (S) does not fail). Togheter with this sentence "With three or fewer outputs, lu uses its own sparse matrix routines" (available in old matlab documentation) I think that [L, U] = lu (S) in matlab does not use umfpack. I will contact umfpack's author.

Marco Caliari <caliari>
Group Member
Wed 21 Mar 2018 09:30:57 PM UTC, comment #4: 

@Rik: is it not a problem that Int (generic integer) is defined as SuiteSparse_long (8 bytes) and size of BLAS integer is 4?

Marco Caliari <caliari>
Group Member
Wed 21 Mar 2018 08:11:15 PM UTC, comment #3: 

It's probably nothing.  I wrote a small C++ program to check the sizes of long, and unsigned long and they are both 8 bytes on my machine.  If there was something like a SuiteSparse_int and it had a different value then 4 bytes then that would be a problem.

Rik <rik5>
Group administrator
Wed 21 Mar 2018 09:17:31 AM UTC, comment #2: 

If you set


Control (UMFPACK_PRL) = 3;


before umfpack_report_*, you get some statistics. This one in case it works


Symbolic object: OK
Numeric object:  OK

column-form matrix, n_row 1000 n_col 1000, nz = 1253. OK

column-form matrix, n_row 1000 n_col 1000, nz = 800. OK

permutation vector, n = 1000. OK

permutation vector, n = 1000. OK

UMFPACK V5.6.2 (Apr 25, 2013), Info:
    matrix entry defined as:          double
    Int (generic integer) defined as: SuiteSparse_long
    BLAS library used: Fortran BLAS.  size of BLAS integer: 4
    MATLAB:                           no.
    CPU timer:                        POSIX C clock_getttime ( ) routine.
    number of rows in matrix A:       1000
    number of columns in matrix A:    1000
    entries in matrix A:              1000
    memory usage reported in:         16-byte Units
    size of int:                      4 bytes
    size of SuiteSparse_long:         8 bytes
    size of pointer:                  8 bytes
    size of numerical entry:          8 bytes
[skipped]


This one in case of failure


Symbolic object: OK
UMFPACK V5.6.2 (Apr 25, 2013): ERROR: pattern of matrix (Ap and/or Ai) has changed

UMFPACK V5.6.2 (Apr 25, 2013), Info:
    matrix entry defined as:          double
    Int (generic integer) defined as: SuiteSparse_long
    BLAS library used: Fortran BLAS.  size of BLAS integer: 4
    MATLAB:                           no.
    CPU timer:                        POSIX C clock_getttime ( ) routine.
    number of rows in matrix A:       1000
    number of columns in matrix A:    1000
    entries in matrix A:              1000
    memory usage reported in:         16-byte Units
    size of int:                      4 bytes
    size of SuiteSparse_long:         8 bytes
    size of pointer:                  8 bytes
    size of numerical entry:          8 bytes
[skipped]


What do you think about the different sizes of SuiteSparse_long and int?

Marco Caliari <caliari>
Group Member
Wed 21 Mar 2018 12:31:18 AM UTC, comment #1: 

Confirmed.  I used the following test script (attached).


for i = 1:100
  S = sprand  (1000, 1000, 1 / 1000);
  [L,U] = lu (S);
  fprintf (stderr, "%d\n", i);
endfor


Octave rarely gets to more than 3 iterations before lu fails.


(file #43610)

Rik <rik5>
Group administrator
Tue 20 Mar 2018 03:27:30 PM UTC, original submission:  

If you try


S = sprand  (1000, 1000, 1 / 1000);
[L, U] = lu (S);


quite often you get


error: lu: sparse_lu: numeric factorization failed


If you check the value of status around line 730 of sparse-lu.cc, you find it is -11 which means (umfpack_numeric.h, umfpack.h)


The pattern (Ap and/or Ai) has changed since the call to
umfpack_*_*symbolic which produced the Symbolic object.


I don't see why it should have changed.

Marco Caliari <caliari>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #43838:  bicg.diff added by caliari (889B - application/x-tex)
file #43610:  tst_lu.m added by rik5 (108B - text/x-matlab)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by rik5
  • -email is unavailable- added by rik5 (Updated the item)
  • -email is unavailable- added by caliari (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 10 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-04-09 rik5 StatusFixed In Progress
        Open/ClosedClosed Open
    2018-04-09 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2018-04-08 caliari Attached File- Added bicg.diff, #43838
    2018-04-07 rik5 Carbon-Copy- Added cristianodorigo
    2018-03-23 caliari StatusNeed Info In Progress
    2018-03-22 rik5 StatusConfirmed Need Info
    2018-03-21 rik5 Attached File- Added tst_lu.m, #43610
        StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code