bugGNU Octave - Bugs: bug #41116, Incorrect result for sparse LU...

 
 

bug #41116: Incorrect result for sparse LU factorization

Submitted by:  None
Submitted on:  Mon 06 Jan 2014 02:33:23 AM UTC  
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Michael C. Grant Originator Email:  -email is unavailable-
Open/Closed:  Closed Release:  3.8.0
Operating System:  Mac OS

Add a New Comment (Rich Markup)
   

You are not logged in

Please log in, so followups can be emailed to you.

 

( Jump to the original submission)

Mon 03 Feb 2014 05:02:43 PM UTC, comment #25: 

Fantastic, thanks. I have applied your patch together with a few stylistic fixes of my own:

http://hg.savannah.gnu.org/hgweb/octave/graph/93c019f00e59

I believe this closes this bug report. The results will be in the upcoming 3.8.1 release.

Jordi GutiƩrrez Hermoso <jordigh>
Project Administrator
Sun 02 Feb 2014 05:57:55 PM UTC, comment #24: 

My apologies, I'm finally getting the hang of Mercurial and your processes. (Yes, I'm slow; I'm just now getting CVX out of Subversion!) The file "lupatch.diff" here is a formal Mercurial changeset generated per section D.3 of the manual. Integrates the patch (built against the current tip) and a commit message with the bug number referenced.

(file #30440)

Michael C. Grant <mcgrant>
Wed 29 Jan 2014 12:16:11 AM UTC, comment #23: 

New patch created that adds a test with this data set. I've run 'make check' to confirm it passes.

Proposed commit message:

correct numerical errors in sparse LU factorization.

* lu.cc: modified to apply pivots as warranted to L and U.
* sparse-base-lu.cc: compute correct matrix size for single-output case.

(file #30400)

Michael C. Grant <mcgrant>
Tue 28 Jan 2014 09:32:07 PM UTC, comment #22: 

I will do both in about 20-30 minutes... hang on

Michael C. Grant <mcgrant>
Tue 28 Jan 2014 09:26:40 PM UTC, comment #21: 

Hey, I'm about to apply this.

One request: can you please turn the original report into a test case? Look at this for examples of how to do it:

http://hg.savannah.gnu.org/hgweb/octave/file/5b88b2cb431e/libinterp/corefcn/lu.cc#l768

Also, would it be too much trouble to provide a commit message?

http://wiki.octave.org/Commit_message_guidelines#Examples

Jordi GutiƩrrez Hermoso <jordigh>
Project Administrator
Tue 28 Jan 2014 01:34:57 AM UTC, comment #20: 

Ping... I would like to be able announce CVX works with 3.8.1... pretty please :-) (Yeah I do have a hack I can use if needed, but I'd really prefer not to!)

Michael C. Grant <mcgrant>
Fri 24 Jan 2014 10:44:17 PM UTC, comment #19: 

Lets do it so that it is in the coming releases.

Michael Godfrey <godfrey>
Project Member
Fri 24 Jan 2014 09:09:46 PM UTC, comment #18: 

Agreed. I am not convinced this exactly mirrors the approach Matlab uses; I am not sure that's possible without bringing in another LU library. But it does seem to coincide exactly most of the time; and at least when it doesn't, we still have that ever-so-desirable property of "correctness" :-)

Michael C. Grant <mcgrant>
Thu 23 Jan 2014 11:58:04 PM UTC, comment #17: 

It seems to me that this should go in an an error-correcting
patch for stable, and all. It may not be perfect, but it
is quite sure that it corrects problems with the current
code which are quite serious.

Michael

Michael Godfrey <godfrey>
Project Member
Thu 23 Jan 2014 11:10:47 PM UTC, comment #16: 

Ping... Sorry to be a pest (no, I'm not :-))... I've been using this patch for awhile without incident, I'd love to see this in 3.8.1. Can we get other testers if needed?

Michael C. Grant <mcgrant>
Thu 16 Jan 2014 05:14:25 PM UTC, comment #15: 

John, any chance you can verify this and commit this into the gui-default build? I'd like to remove the workaround code I have in CVX before I push out an Octave-compatible alpha package.

Michael C. Grant <mcgrant>
Wed 15 Jan 2014 05:54:12 AM UTC, comment #14: 

And I have to be honest, I have no clue why anyone would ever be interested in the nargout=1 version of lu(). It discards pivoting so when is it ever practical to use it?

Michael C. Grant <mcgrant>
Wed 15 Jan 2014 05:52:07 AM UTC, comment #13: 

I am uploading a patch for both lu.cc and sparse-base-lu.cc. It turns out that the nargout=1 mode wasn't being constructed properly, either: the size of the output matrix is supposed to be the same size as the original, but for matrices with more columns than rows, it was square.

I'm basically hoping against hope here that Marco is right, and that the column reorderings performed during UMFPACK's numerical factorization step will not destroy the upper triangular nature of U. When that assumption is correct, this should indeed produce results identical to MATLAB's.

If that assumption fails, then at least we still have x = LL * UU. And since Octave knows that UU is a permuted upper triangular matrix anyway, there shouldn't be a performance loss if someone is depending on triangularity simply to achieve high solve performance. The only people that might get tripped up are those that depend on UU's upper triangularity, and I don't actually know who such people might be.

(file #30277)

Michael C. Grant <mcgrant>
Tue 14 Jan 2014 02:03:55 PM UTC, comment #12: 

Hmm, you might be right. Here's the documentation: " [L,U,P] = lu(A) returns unit lower triangular matrix L, upper triangular matrix U, and permutation matrix P so that P*A = L*U."

For this particular value of x, you're right, even with the column permutations, the upper triangular property is preserved. My concern was that there might be scenarios where the column permutations would destroy the upper triangular property of the matrix.

Certainly, it is possible for an arbitrary upper triangular matrix U2 and an arbitrary column permutation Q2 that the matrix U=U2*Q2' would fail to be upper triangular. But perhaps the way UMFPACK chooses its numerical swaps, this would not happen.

Michael C. Grant <mcgrant>
Tue 14 Jan 2014 01:55:57 PM UTC, comment #11: 

What about the controls UMFPACK_ORDERING_GIVEN and UMFPACK_ORDERING_NONE?

Anyway, [LL,UU] = lu(x) with my patch (i.e., with column permutation) behaves like [LL,UU] = lu(x) in Matlab, which is supposed to work without column permutation. Are we sure that column permutation is completely disabled in Matlab?

Marco

Marco Caliari <caliari>
Project Member
Tue 14 Jan 2014 12:44:56 PM UTC, comment #10: 

Oh sure, I agree that LL*UU=x with your modification. And indeed that may be the best we can do. I think your path is probably what should be done for now.

But MATLAB ensures that UU is upper triangular when nargout<=3---it expects column permutation to be completely disabled. That behavior cannot be reproduced with UMFPACK.

So there will still remain MATLAB incompatibility unless we bring in a different sparse LU factorization.

Michael C. Grant <mcgrant>
Tue 14 Jan 2014 08:34:34 AM UTC, comment #9: 

Sorry, I rewrite

Dear Michael,

[LL2,UU2]=lu(full(x));

makes no column permutation: in fact LL2*UU2=x, without any Q matrix.
With the workaround I suggested, I get from

[LL,UU]=lu(x);

exactly the same (correct) factorization, that is no column permutation in the
final result. It seems to me that the column reordering we see between the
symbolic and the numerical factorization in UMFPACK is a sort of "inner"
reordering that we have to take into account, but without a real column
permutation (that is Q different from identity) in the result.

Marco

Marco Caliari <caliari>
Project Member
Tue 14 Jan 2014 08:31:22 AM UTC, comment #8: 

Dear Michael,

[LL2,UU2]=lu(full(x));

makes no column permutation: in fact LL2*UU2=x, without any Q matrix.
With the workaround I suggested, I get from

Marco Caliari <caliari>
Project Member
Mon 13 Jan 2014 08:26:53 PM UTC, comment #7: 

I now suspect that MATLAB does not use UMFPACK for LU factorization when nargout=1,2, or 3. I don't think you can use UMFPACK in those cases.

It looks like UMFPACK reserves the right to do column reordering during numerical factorization, even if you force a fixed column ordering during symbolic factorization. I see no way to disable that.

If I am right, it is not possible to duplicate MATLAB functionality for nargout=1,2,3 without replacing UMFPACK with something else; say, SuperLU.

Michael C. Grant <mcgrant>
Mon 13 Jan 2014 04:31:08 PM UTC, comment #6: 

Actually, I was wrong: apparently, when nargout==2 or nargout==3, there is not supposed to be any column permutation at all. I'll figure out what hast to be changed to prevent column pivoting.

Michael C. Grant <mcgrant>
Thu 09 Jan 2014 12:33:02 PM UTC, comment #5: 

I have not tested this yet, but I can confirm that you have made the correct interpretation. This is how the factorization is supposed to work.

Michael C. Grant <mcgrant>
Thu 09 Jan 2014 08:48:34 AM UTC, comment #4: 

Dear all,

I found that q after UMFPACK_DNAME (get_numeric) around l.416 of SparsedbleLU.cc is not the identity permutation. I don't know if this is the wanted behavior. In fact, UMFPACK_FIXQ = 1.0 only means that the symbolic column permutation is kept in the numeric factorization.
Anyway, since there is a column permutation, I had to rewrite the code around l.233 of lu.cc and permute U

PermMatrix Q = fact.Pc_mat ();
SparseMatrix U = fact.U () * Q.transpose ();
retval(1) = octave_value (U, MatrixType (MatrixType::Permuted_Upper, nr, fact.col_perm ()));

I'm not sure of the last line, I copied it from the row permutation of L and changed row_perm with col_perm (???)
The result is now correct. If this is a real fix and not a simple workaround, than it should be applied to the 3-argument version and to the float, complex versions.

I hope someone more confident with the code can use my investigation to provide a complete patch.

Marco

Marco Caliari <caliari>
Project Member
Mon 06 Jan 2014 03:25:15 PM UTC, comment #3: 

I have confirmed that the three-argument version is also incorrect, but the four-argument version is correct. I've studied the LU code a bit but I don't know enough about UMFPACK to figure it out, and I have another project calling. I'll get back to this when I can, unless someone else does before me. I would love to see CVX working with Octave 4.0 or earlier; we're very very close!

load lutest
[LL2,UU2]=lu(x);
x-LL2*UU2
[LL3,UU3,PP3]=lu(x);
x-PP3'*LL3*UU3
[LL4,UU4,PP4,QQ4]=lu(x);
x-PP4'*LL4*UU4*QQ4'

Michael C. Grant <mcgrant>
Mon 06 Jan 2014 04:15:54 AM UTC, comment #2: 

I'm going to hazard a guess that this occurs only when the matrix has more columns than rows. This produces a correct results:

load lutest
[UU,LL]=lu(x');
UU=UU';
LL=LL';
LL*UU-x

Michael C. Grant <mcgrant>
Mon 06 Jan 2014 03:33:41 AM UTC, comment #1: 

I see the same result going all the way back to version 3.2.4, so I guess this has always been done incorrectly.  I'm surprised no one noticed before now.

John W. Eaton <jwe>
Project Administrator
Mon 06 Jan 2014 02:33:23 AM UTC, original submission:  

load lutest
[LL,UU]=lu(x);
[LL2,UU2]=lu(full(x));
LL*UU-x % this should be zero, but is not
LL2*UU2-x % this is zero
LL-LL2 % this is zero
UU-UU2 % this should be zero, but is not

Inspection reveals that the columns of UU remain un-permuted. But using the three-argument version of LU does not seem to help:

[LL3,UU3,PP3]=lu(x);
LL3*UU3-PP3*x

PP3 is actually the identity matrix.
I have verified this both with a version of Octave 3.8 compiled with Homebrew, and the prebuilt binaries from Octave Forge.

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #30440:  lupatch.diff added by mcgrant (9KiB - application/octet-stream)
file #30400:  patch2.diff added by mcgrant (8KiB - application/octet-stream)
file #30277:  patch.diff added by mcgrant (8KiB - application/octet-stream)
file #30199:  lutest added by None (248B - 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 jordigh (Posted a comment)
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by mcgrant (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

     

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-02-03 jordigh StatusPatch Submitted => Fixed
        Open/ClosedOpen => Closed
    2014-02-02 mcgrant Attached File- => Added lupatch.diff, #30440
    2014-01-29 mcgrant Attached File- => Added patch2.diff, #30400
    2014-01-28 jordigh StatusConfirmed => Patch Submitted
    2014-01-15 mcgrant Attached File- => Added patch.diff, #30277
    2014-01-06 jwe StatusNone => Confirmed
    2014-01-06 None Attached File- => Added lutest, #30199

    Back to the top


    Powered by Savane 3.5