bugGNU Octave - Bugs: bug #34850, Behaviour of chol(a,'lower')

 
 

bug #34850: Behaviour of chol(a,'lower')

Submitted by:  None
Submitted on:  Thu 17 Nov 2011 02:28:28 PM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Performance
Status: FixedAssigned to: Carlo de Falco <cdf>
Originator Name: Salvatore FilipponeOriginator Email: -unavailable-
Open/Closed: ClosedRelease: dev
Operating System: Any

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 Jump to the original submission)

Tue 25 Aug 2015 02:20:02 AM UTC, comment #30:

Thanks Carlo for accepting the patch.

PrasannaKumar Muralidharan <prasannatsm>
Mon 24 Aug 2015 02:49:19 PM UTC, comment #29:

Hi,

I finally got around to test the patch,
it works correctly and provides a speed up
of about 7% for a 2e3 x 2e3 complex matrix
on my system.

I pushed the changeset to the default branch with this changeset:
http://hg.savannah.gnu.org/hgweb/octave/rev/5ce959c55cc0

thanks,
c.

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Mon 24 Aug 2015 05:59:08 AM UTC, comment #28:

Any update on this patch?

PrasannaKumar Muralidharan <prasannatsm>
Thu 22 Jan 2015 06:08:48 AM UTC, comment #27:

@Carlo: Did you get a chance to test the patch?

PrasannaKumar Muralidharan <prasannatsm>
Sun 28 Dec 2014 09:56:12 AM UTC, comment #26:

HI,

The patch seems OK but I cannot test at the moment as I am in the process of upgrading my OS so I won't be able to try and compile Octave again until tomorrow.

I'll report back when done.

c.

Anonymous
Sat 27 Dec 2014 05:37:57 PM UTC, comment #25:

@Carlo: Have you had a chance to review the attached patch? It seems like a worthwhile change to rely on LAPACK rather than doing transposes ourselves in Octave.

Rik <rik5>
Project Administrator
Wed 01 Oct 2014 10:07:05 PM UTC, comment #24:

I have been too busy with other stuff for a while.
I'll look into this ASAP

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Sun 28 Sep 2014 09:10:42 AM UTC, comment #23:

Any plans to merge this?

PrasannaKumar Muralidharan <prasannatsm>
Sun 24 Aug 2014 03:50:23 PM UTC, comment #22:

I have implemented the changes for double, float, complex types. Attached patch with this comment.

chol test passes without issues.

(file #31952)

PrasannaKumar Muralidharan <prasannatsm>
Fri 11 Jul 2014 01:17:38 PM UTC, comment #21:

DPOTRI needs to come after a call to DPOTRF and needs exactly the same value for UPPER.

Salvatore

Salvatore Filippone <sfilippone>
Fri 11 Jul 2014 11:42:37 AM UTC, comment #20:

is_upper is a member of class CHOL. Value of 'upper' is stored in 'is_upper'. It is later used in CHOL::inverse.

Is it not necessary while computing inverse (using dpotri lapack function)? If not necessary I can make necessary changes.

Thanks for the review.

PrasannaKumar Muralidharan <prasannatsm>
Fri 11 Jul 2014 10:33:32 AM UTC, comment #19:

Looks OK to me at a first look, except I don't understand what's the point in:

why not

?

I'm now testing the patch,
thanks for this contribution.
c.

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Tue 08 Jul 2014 04:32:51 PM UTC, comment #18:

I have made changes to send "L" to BLAS functions (for real matrices only). The output code seems to be working.

Please review / test the change and if everything works fine I will try to make changes to the other types (complex, float complex etc).

make check does not report any new error after this change. Validated sample input available in the description of this bug.

(file #31691)

PrasannaKumar Muralidharan <prasannatsm>
Mon 20 Jan 2014 04:57:55 PM UTC, comment #17:

Thanks for relaying that information Carlo. I'll set status to "in progress" for now? Also setting as a performance item, since it now works correctly but needs work to call the correct numerical library calls.

Mike Miller <mtmiller>
Project Administrator
Mon 20 Jan 2014 04:52:34 PM UTC, comment #16:

I'm reposting here the email reply by Salvatore as
the tracker does not automatically record replies sent via email.
c.

No, I did not have time then, and I did not work on the codes. I can
confirm that the workaround still works (as of the version I have
installed right now).
Looking at the source files the fix requires changing
dbleCHOL::init (const Matrix& a, bool calc_cond)

and its friends to accept the lower/upper flag and pass it to the
inner call. to the LAPACK routine.
Something along the lines of the code I am attaching should work, but
I did not test it, and I cannot do it any time soon. But if Carlo
has time to test I can look at the results.
NOTE: for complex data, in cholinv_internal you'll need to add std::conj .

BTW, I am not sure what happens in the sparse case...that's a whole
different game, and we'd need to look at the UMFPACK docs.
Salvatore

(file #30318, file #30319)

Anonymous
Mon 20 Jan 2014 09:27:14 AM UTC, comment #15:

The "bug" has not been really resolved
but merely "worked around" by applying
the patch in file #24462 .

A "clean" fix would require using the correct
LAPACK call rather than transposing the matrix
twice as is currently done.

Salvatore did you have time for looking into this?

c.

Anonymous
Sun 19 Jan 2014 10:51:57 PM UTC, comment #14:

Has this bug been resolved? The test cases provided for real and complex inputs to the chol function seem to return the results listed as correct in this report, tested in the 3.8.0 release.

Mike Miller <mtmiller>
Project Administrator
Mon 28 Nov 2011 01:21:40 PM UTC, comment #13:

Ok, time is a limited resource but this shouldn't be too hard.
Regards
Salvatore

Salvatore Filippone <sfilippone>
Mon 28 Nov 2011 10:18:27 AM UTC, comment #12:

I Agree with both parts of this message.

So I'll commit this temporary change in order
to try have it included in version 3.6
which is soon to be released and start working
on a new patch to better exploit LAPACK functionality.

Would you be willing to help with this? The files that would need to be modified are the the following:

chol.cc : main driver and user interface
http://hg.savannah.gnu.org/hgweb/octave/file/acaf33ccc04f/src/DLD-FUNCTIONS/chol.cc

*CHOL.[cc|h] : wrappers for LAPACK functions
FloatCHOL.[cc|h]
FloatCmplexCHOL.[cc|h]
CHOL.[cc|h]
CmplexCHOL.[cc|h]

http://hg.savannah.gnu.org/hgweb/octave/file/acaf33ccc04f/liboctave/FloatCHOL.cc

http://hg.savannah.gnu.org/hgweb/octave/file/acaf33ccc04f/liboctave/FloatCmplexCHOL.cc

http://hg.savannah.gnu.org/hgweb/octave/file/acaf33ccc04f/liboctave/CHOL.cc

http://hg.savannah.gnu.org/hgweb/octave/file/acaf33ccc04f/liboctave/CmplxCHOL.cc

Similar functions for Sparse matrices also exist but do not use LAPACK, so we should test they work correctly but probably they won't need fixing.

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Mon 28 Nov 2011 09:59:33 AM UTC, comment #11:

Well, that applying transpose() twice is safe should be clear: if A is Hermitian positive definite, so is A'=conjg(A), therefore for the purpose of having a "quick and dirty" fix this is OK.
It is not quite satisfactory in that we are looping around a problem of the interface, when the underlying engine (LAPACK) is perfectly capable of handling the requirement.

Salvatore

Salvatore Filippone <sfilippone>
Sat 26 Nov 2011 01:54:43 PM UTC, comment #10:

After a bit more testing it appears that it is OK
to use transpose() rather than hermitian() when the
operation is applied twice:

A = randn (4) + i * randn (4);
As = A' * A;

all (all (chol (triu (As'), 'upper')' == ...
chol (triu (As.'), 'upper').'))

ans = 1

so the attached patch works.
Salvatore, would you be able to test it too?

(file #24462)

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Mon 21 Nov 2011 11:13:12 AM UTC, comment #9:

Salvatore,

Thanks for providing a test.
So the behavior of chol is indeed incorrect
for complex data as I suspected, I'll try to fix this ASAP.

c.

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Mon 21 Nov 2011 09:49:35 AM UTC, comment #8:

Ok, for complex data the situation is even worse: the current code returns the transpose of the upper factor instead of the conjugate transpose. Test case:
b=[0,1,2,3;

> -1,0,1,2;
> -2,-1,0,1;
> -3,-2,-1,0]

b =

0 1 2 3
-1 0 1 2
-2 -1 0 1
-3 -2 -1 0

octave:5> ca=a+i*b
ca =

12 + 0i 2 + 1i 3 + 2i 4 + 3i
2 - 1i 14 + 0i 5 + 1i 3 + 2i
3 - 2i 5 - 1i 16 + 0i 6 + 1i
4 - 3i 3 - 2i 6 - 1i 16 + 0i

octave:8> chol(ca)
ans =

3.46410 + 0.00000i 0.57735 + 0.28868i 0.86603 + 0.57735i 1.15470 + 0.86603i
0.00000 + 0.00000i 3.68556 + 0.00000i 1.17576 + 0.24872i 0.56527 + 0.49744i
0.00000 + 0.00000i 0.00000 + 0.00000i 3.67048 + 0.00000i 1.01122 + 0.12870i
0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i 3.50864 + 0.00000i

octave:9> chol(ca,'lower')
ans =

3.46410 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i
0.57735 + 0.28868i 3.68556 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i
0.86603 + 0.57735i 1.17576 + 0.24872i 3.67048 + 0.00000i 0.00000 + 0.00000i
1.15470 + 0.86603i 0.56527 + 0.49744i 1.01122 + 0.12870i 3.50864 + 0.00000i

whereas Matlab correctly returns
ca=a+i*b

ca =

12.0000 2.0000 + 1.0000i 3.0000 + 2.0000i 4.0000 + 3.0000i
2.0000 - 1.0000i 14.0000 5.0000 + 1.0000i 3.0000 + 2.0000i
3.0000 - 2.0000i 5.0000 - 1.0000i 16.0000 6.0000 + 1.0000i
4.0000 - 3.0000i 3.0000 - 2.0000i 6.0000 - 1.0000i 16.0000

>> chol(ca)


ans =

3.4641 0.5774 + 0.2887i 0.8660 + 0.5774i 1.1547 + 0.8660i
0 3.6856 1.1758 + 0.2487i 0.5653 + 0.4974i
0 0 3.6705 1.0112 + 0.1287i
0 0 0 3.5086

>> chol(ca,'lower')


ans =

3.4641 0 0 0
0.5774 - 0.2887i 3.6856 0 0
0.8660 - 0.5774i 1.1758 - 0.2487i 3.6705 0
1.1547 - 0.8660i 0.5653 - 0.4974i 1.0112 - 0.1287i 3.5086

Salvatore Filippone <sfilippone>
Sun 20 Nov 2011 06:40:59 PM UTC, comment #7:

Sorry, posted too soon, I will check Matlab not so much to verify that hermitian is needed (that's certain), but to provide a viable test case.

Salvatore

Salvatore Filippone <sfilippone>
Sun 20 Nov 2011 06:39:46 PM UTC, comment #6:

Well, I will check Matlab as soon as possible, but to duplicate what LAPACK would be doing on 'L', x.hermitian() is definitely needed on both input and output.

Salvatore

Salvatore Filippone <sfilippone>
Sun 20 Nov 2011 06:33:42 PM UTC, comment #5:

Hi,
thanks for your feedback

>> I know nothing of Octave internals, but it seems to me that the
>> most logical way would be to propagate the "lower" choice down
>> to LAPACK; this also saves the transpose (which takes time).


Indeed that would be cleaner, but it requires changing quite
a large number of files, accounting for all possible base types,
e.g.:

dbleCHOL.[h|cc]
floatCHOL.[h|cc]
CmplxCHOL.[h|cc]
fCmplxCHOL.[h|cc]

and so on ...

As I said I have to think a little bit myself about how to do that
without too much code duplication and it will take time, that's why I'm
proposing a temporary patch.

>> I have trouble now looking at the patch something is worng with
>> the connection, but about a hour ago I got a glimpse, and I
>> remember you are calling X.transpose() for complex data;
>> in that case you want to have the conjugate transpose to feed ZPOTRF.


You are probably right on this one, X.hermitian () would be correct here.

But I think this shows another problem with the current implementation:
what chol.cc currently does is call xPOTRF with UPLO = 'U' and then
return the transpose of the factor,
I think for x = [c|z] it should actually also return the hermitian,
is this correct? could you check if chol ( ... ,'lower') currently
works as you expect for complex types?

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Sun 20 Nov 2011 05:44:33 PM UTC, comment #4:

Hi,
I know nothing of Octave internals, but it seems to me that the most logical way would be to propagate the "lower" choice down to LAPACK; this also saves the transpose (which takes time).
I have trouble now looking at the patch something is worng with the connection, but about a hour ago I got a glimpse, and I remember you are calling X.transpose() for complex data; in that case you want to have the conjugate transpose to feed ZPOTRF.

Salvatore

Salvatore Filippone <sfilippone>
Sun 20 Nov 2011 04:54:05 PM UTC, comment #3:

A proposed patch is attached.

(file #24403)

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Sun 20 Nov 2011 03:40:28 PM UTC, comment #2:

More info on this issue.

It seems that the lapack function DPOTRF is always invoked with
the input parameter UPLO set to 'U':

http://hg.savannah.gnu.org/hgweb/octave/file/5b1ddcf5ede1/liboctave/dbleCHOL.cc#l119

then, if the input argument 'lower' is passed to chol in Octave,
the factor returned by lapack is transposed

http://hg.savannah.gnu.org/hgweb/octave/file/5b1ddcf5ede1/src/DLD-FUNCTIONS/chol.cc#l55

So it seems the behavior observed by Salvatore is a feature rather than a bug.

One quick solution to get more compatible behavior would be to transpose the input to chol before factoring it, but I'll have to think a bit to find a cleaner way to fix this.

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Sun 20 Nov 2011 02:53:08 PM UTC, comment #1:

The same behavior is still present in the development release.

Carlo de Falco <cdf>
Project MemberIn charge of this item.
Thu 17 Nov 2011 02:28:28 PM UTC, original submission:

Hello,
Following the discussion https://mailman.cae.wisc.edu/pipermail/help-octave/2011-November/048569.html

To reproduce the problem: Build an SPD matrix
octave:1> a=[12,2,3,4;

> 2,14,5,3;
> 3,5,16,6;
> 4,3,6,16]

a =

12 2 3 4
2 14 5 3
3 5 16 6
4 3 6 16

Ask for the Cholesky factorization: in LAPACK, if you ask for the 'Upper' factor, only the upper triangle counts and is looked at/modified, and Octave agrees as in the following:

> chol(triu(a))

ans =

3.46410 0.57735 0.86603 1.15470
0.00000 3.69685 1.21725 0.63117
0.00000 0.00000 3.71057 1.14045
0.00000 0.00000 0.00000 3.60107

The symmetric choice of 'lower' in LAPACK only looks at the lower triangle; this is also true of MATLAB, but Octave gives the following:
octave:5> chol(tril(a),'lower')
ans =

3.46410 0.00000 0.00000 0.00000
0.00000 3.74166 0.00000 0.00000
0.00000 0.00000 4.00000 0.00000
0.00000 0.00000 0.00000 4.00000

i.e. it behaves as if factoring the upper triangle (which is just the main diagonal) and then transposing it, whereas according to the LAPACK mode of operation, I should get exactly* the same output I would get by asking for
(chol(triu(a)))'

The documentation of CHOL does not explicitly specify this behaviour, whereas LAPACK is very explicit about this.
As a minimum, this should be clarified in the docs (I'd prefer to see complete compatibility with MATLAB and LAPACK though).

Thanks
Salvatore

  • Note: up to working precision, as going from Upper to Lower might change the order of floating point operations, and thus there may be differences in the last (few) significant bit(s), also depending on the matrix condition. In the simple example above there is no difference.
Anonymous

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #31691:  dble_chol_propogate_lower.diff added by prasannatsm (8KiB - text/x-patch - Patch propogating lower to underlying library function)
file #30318:  dbleCHOL.cc added by None (11KiB - application/octet-stream)
file #30319:  dbleCHOL.h added by None (2KiB - application/octet-stream)
file #24462:  open_8Hshmam0.txt added by cdf (5KiB - text/plain)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by prasannatsm (Updated the item)
  • -unavailable- added by mtmiller (Posted a comment)
  • -unavailable- added by sfilippone (Posted a comment)
  • -unavailable- added by cdf (Posted a comment)
  • -unavailable- added by None (Submitted the item)
  •  

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

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 22 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon 24 Aug 2015 02:49:19 PM UTCcdfStatusIn Progress=>Fixed
      Open/ClosedOpen=>Closed
    Sun 24 Aug 2014 03:50:23 PM UTCprasannatsmAttached File-=>Added propagate_lower_chol_dbl_flt_cmplx.diff, #31952
    Tue 08 Jul 2014 04:32:51 PM UTCprasannatsmAttached File-=>Added dble_chol_propogate_lower.diff, #31691
    Mon 20 Jan 2014 04:57:55 PM UTCmtmillerItem GroupMatlab Compatibility=>Performance
      StatusNeed Info=>In Progress
    Mon 20 Jan 2014 04:52:34 PM UTCNoneAttached File-=>Added dbleCHOL.cc, #30318
      Attached File-=>Added dbleCHOL.h, #30319
    Sun 19 Jan 2014 10:51:57 PM UTCmtmillerCategoryNone=>Octave Function
      StatusPatch Submitted=>Need Info
      Operating SystemGNU/Linux=>Any
    Sat 26 Nov 2011 02:03:41 PM UTCcdfStatusIn Progress=>Patch Submitted
    Sat 26 Nov 2011 01:58:37 PM UTCcdfAttached File#24404=>Removed
    Sat 26 Nov 2011 01:54:43 PM UTCcdfAttached File-=>Added open_8Hshmam0.txt, #24462
    Mon 21 Nov 2011 11:13:12 AM UTCcdfStatusPatch Submitted=>In Progress
    Sun 20 Nov 2011 05:14:12 PM UTCcdfAttached File-=>Added open_cQbdpVI2.txt, #24404
    Sun 20 Nov 2011 05:13:48 PM UTCcdfAttached File#24403=>Removed
    Sun 20 Nov 2011 04:55:09 PM UTCcdfStatusConfirmed=>Patch Submitted
    Sun 20 Nov 2011 04:54:05 PM UTCcdfAttached File-=>Added open_iTFnjem2.txt, #24403
    Sun 20 Nov 2011 02:53:08 PM UTCcdfStatusNone=>Confirmed
      Assigned toNone=>cdf
      Release3.4.2=>dev

    Back to the top


    Powered by Savane 3.1-cleanup1