Tue 25 Aug 2015 02:20:02 AM UTC, comment #30:
Thanks Carlo for accepting the patch.
|
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.
|
Mon 24 Aug 2015 05:59:08 AM UTC, comment #28:
Any update on this patch?
|
Thu 22 Jan 2015 06:08:48 AM UTC, comment #27:
@Carlo: Did you get a chance to test the patch?
|
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.
|
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.
|
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
|
Sun 28 Sep 2014 09:10:42 AM UTC, comment #23:
Any plans to merge this?
|
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)
|
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
|
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.
|
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.
|
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)
|
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.
|
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)
|
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.
|
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.
|
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
|
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.
|
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
|
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)
|
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.
|
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
|
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
|
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
|
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?
|
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
|
Sun 20 Nov 2011 04:54:05 PM UTC, comment #3:
A proposed patch is attached.
(file #24403)
|
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.
|
Sun 20 Nov 2011 02:53:08 PM UTC, comment #1:
The same behavior is still present in the development release.
|
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.
|