bugGNU Octave - Bugs: bug #29768, logm fix

 
 

bug #29768: logm fix

Submitter:  Jaroslav Hajek <highegg>
Submitted:  Tue 04 May 2010 08:52:07 AM UTC
   
 
Category:  Libraries Severity:  2 - Minor
Priority:  5 - Normal Item Group:  Feature Request
Status:  Fixed Assigned to:  highegg
Originator Name:  Tommy Guy Open/Closed:  * Closed
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

Wed 28 Jul 2010 07:01:26 PM UTC, comment #31: 

I added the GPL license block to the head of the file since Prof. Higham agreed to anything covered by GPL version 3.  I also made a few changes to the documentation string.

I'd say this feature has been implemented and I'm going to close this report.

Rik <rik5>
Group administrator
Wed 28 Jul 2010 07:24:46 AM UTC, comment #30: 

I polished the code somewhat, taking into account Rik's remarks, Octave's style as well as my own ideas; then pushed a version to Savannah:
http://hg.savannah.gnu.org/hgweb/octave/rev/cace99cb01ab

Please review and comment.

Jaroslav Hajek <highegg>
Wed 28 Jul 2010 07:02:27 AM UTC, comment #29: 

Dear Rik,

I learnt something new about "find", thanks!
I agree to set m=7 as default value, too.

Marco

Marco Caliari <caliari>
Group Member
Tue 27 Jul 2010 08:05:42 PM UTC, comment #28: 

Marco,

I had two thoughts on the code.

Where you have written

min (find (tau <=  theta))

you could use

find (tau <= theta, 1)

Which finds the index of the first occurrence that matches.  This might be clearer than using min in this case.

Second, a value for m is assigned only if the tau is less than theta(7) and other conditions are met.  It would probably be a good idea to assign a value to it before the loop is entered in case opt_iters has been set very low and the loop exits before assigning a value to m.  What is the proper value in this case?  7?

Rik <rik5>
Group administrator
Tue 27 Jul 2010 07:52:01 PM UTC, comment #27: 

The recoding to use m < 16 for the Pade section really speeded things up.  On the same random test matrix the benchmark timings were:

Elapsed time is 1.554 seconds.
Elapsed time is 0.009721 seconds.
Elapsed time is 2.249 seconds.
Elapsed time is 1.706 seconds.
Elapsed time is 0.7743 seconds.

The Pade section changed from 4.4 seconds to 1.7 seconds.
Overall the code is now just 2.4X longer than the old eigenvalue-based code.

Rik <rik5>
Group administrator
Tue 27 Jul 2010 01:17:32 PM UTC, comment #26: 

As said, I implemented Algorithm 11.9 in Higham's book. The values for theta(3:6) are taken from Table 11.1 (which reports them with three digits, so I lowered the third by one, to be safe).

Marco

(file #21084)

Marco Caliari <caliari>
Group Member
Tue 27 Jul 2010 11:38:04 AM UTC, comment #25: 

I think prof. Higham was already contacted (see earlier posts) and he reportedly agreed to anything within the bounds of GPL. So there is no need, I suppose.

Jaroslav Hajek <highegg>
Tue 27 Jul 2010 11:32:12 AM UTC, comment #24: 

I agree that rsf2csf should only check that UR and TR are square and of the same dimension.

Coming back to logm, I see that m==16 for the Pade' approximation. In algorithm 11.9 of the book "Function of matrices" by N. Higham there is a simple strategy to reduce m to 7 at most. This would reduce a lot the cost of Pade' approximation.

Do we need a written permission to use it?

Marco

Marco Caliari <caliari>
Group Member
Tue 27 Jul 2010 10:37:59 AM UTC, comment #23: 

Aha, I see. So the formulae can be slightly simpler.
Do you think rsf2csf should perform some checks whether the matrices are of the correct form?
Usually where a valid factorization is expected (e.g. qrupdate) we only check dimensions. Verifying that a matrix is unitary would be very costly anyway, so probably no checks are needed.

Jaroslav Hajek <highegg>
Tue 27 Jul 2010 09:52:00 AM UTC, comment #22: 

The diagonal blocks of real Schur are of the form

[alpha, beta;gamma, alpha]

with beta * gamma < 0 (otherwise there are real eigenvalues). Then, beta / (beta - gamma) is always greater than zero.

Marco

Marco Caliari <caliari>
Group Member
Tue 27 Jul 2010 09:37:29 AM UTC, comment #21: 

Hi Marco,

I found out I have already played with this a couple of months before and had an incomplete Fortran implementation, so I completed and pushed it. This code stores the rotations in helper vectors so that it can process both matrices strictly column-wise to be cache-friendly.

I invite you to verify the new function. One thing I would like to note: you used

      s = sqrt(T(m,m-1)/(T(m-1,m)-T(m,m-1)));
      c = sqrt(T(m-1,m)/(T(m-1,m)-T(m,m-1)));

I do not think, however, that T(m-1,m) > 0 is guaranteed as an output of real schur, is it? In that case, with c complex,
      G = [c, s;-s', c];
would not be a rotation...

I used instead the formulas

      c = sqrt(abs(T(m-1,m)/(T(m-1,m)-T(m,m-1))));
      s = sqrt(abs(T(m,m-1)/(T(m-1,m)-T(m,m-1)))) * sign(T(m-1,m));

the sign of s should ensure that the conjugate with positive imaginary part comes first.

I will let schur() and sqrtm() take advantage of the new feature.

Jaroslav Hajek <highegg>
Tue 27 Jul 2010 07:43:19 AM UTC, comment #20: 

Dear Jaroslav,

since you are going to translate to C++, then it could be faster to locally apply the Givens rotations. I attach the global transformation approach (rsf2csf_g.m) and the local one (rsf2cst_l.m). I also changed the rotation to be in the form [c,s;-s',c] and not [c,s;-s,c'].

(file #21080, file #21081)

Marco Caliari <caliari>
Group Member
Mon 26 Jul 2010 06:05:47 PM UTC, comment #19: 

Thanks Marco, I'll try to convert this into a C++ version (to be used in the ComplexSCHUR classes). The schur (A, "complex") can take this approach.

Jaroslav Hajek <highegg>
Mon 26 Jul 2010 05:35:21 PM UTC, comment #18: 

I added the rsf2csf function call into the logm script and annotated it with tic/toc pairs to benchmark the various sections.  The modified logm is attached.

Benchmarking on a 500x500 random matrix
Schur decomposition (real)
  Elapsed time is 1.505 seconds.
rsf2csf
  Elapsed time is 0.04366 seconds.
sqrtm
  Elapsed time is 2.221 seconds.
Pade approximation
  Elapsed time is 4.461 seconds.
Recovery
  Elapsed time is 0.7695 seconds.

Overall the time on my machine went from 13 seconds to 9 seconds or about a -30% reduction.  The results also show that the next most time consuming block which needs work is the Pade approximation.

Using the old eigenvector code the time was 2.665 seconds.  I still like my suggestion in comment #11 to split the calculation based on tolerance.  If the user desires maximal accuracy then they can wait for the Pade approximation.  Otherwise, give them something fast that they can work with.

(file #21075)

Rik <rik5>
Group administrator
Mon 26 Jul 2010 07:30:44 AM UTC, comment #17: 

Here it is my rsf2csf.m. I found a little bit faster to collect all the complex Givens rotations into a single G matrix and then transform than to transform after each rotation, as in the code you pointed out.

Cheers,

Marco

(file #21067)

Marco Caliari <caliari>
Group Member
Fri 23 Jul 2010 11:35:21 PM UTC, comment #16: 

As a rule, yes, we're always interested in new code.  I don't know why the original rsf2csf function was removed from Octave-Forge.  You would probably need to query someone over at the Octave Forge site. 

You might want to take a look at this posting from 2001 which has an implementation including some benchmarking:
http://velveeta.che.wisc.edu/octave/lists/bug-octave/2001/80.

Rik <rik5>
Group administrator
Fri 23 Jul 2010 07:48:22 AM UTC, comment #15: 

I have seen that rsf2csf.m was part of Octave, but removed for "license issues" http://www-old.cae.wisc.edu/pipermail/help-octave/2007-September/005923.html. What was the problem? I think I can write such a function from scratch. Are you interested?

Marco

Marco Caliari <caliari>
Group Member
Mon 21 Jun 2010 08:34:08 AM UTC, comment #14: 

Notice the huge increase of time consumed by th eschur decomposition itself. This is because schur(..., "complex") is performed completely in complex arithmetics, which is quite wasteful. A much better approach would be to start with the real factorization and convert it to the complex one by factoring 2x2 diagonal blocks. Matlab even has a function for it, rsf2csf.

Jaroslav Hajek <highegg>
Thu 17 Jun 2010 10:13:44 PM UTC, comment #13: 

This feels like squeezing a balloon.  See the following timed runs for not using/using the "complex" argument to schur.

Elapsed time is 1.629 seconds.
warning: Matrix has nonegative eigenvalues.  Principal matrix logarithm is not defined.  Computing non-principal logarithm.
Elapsed time is 6.846 seconds.
Elapsed time is 5.34 seconds.
Elapsed time is 0.5854 seconds.

Elapsed time is 4.165 seconds.
Elapsed time is 2.662 seconds.
Elapsed time is 4.447 seconds.
Elapsed time is 0.7709 seconds.

In the first case, the schur decomposition is very fast but the sqrtm is exceedingly slow.  For the second case, schur is very slow and sqrtm is 3x faster.  Overall the "complex" argument helped but only reduced total running time by 11% on an average of 10 random 500x500 matrices. 

I guess we wait for comment on accuracy vs. speed.  Another thought is to implement a tolerance or have a way for the user to select the algorithm.

Rik <rik5>
Group administrator
Thu 17 Jun 2010 06:29:02 AM UTC, comment #12: 

Ugh, yes. One thing is that we need a complex schur instead of the real (quasi-triangular) one; without that, the triangular sqrtm optimizations won't take place. The newest Octave recognizes the Matlab-compatible syntax schur (a, "complex"); that alone should give you a much better performance. On a random 500x500 matrix, it cuts down time from 14s to 4.8s for me.

It still lags behind the original logm, however. Measuring the individual stages, I get:

octave:4> logm (a);
schur decomposition
Elapsed time is 1.62743 seconds.
square roots
Elapsed time is 1.76724 seconds.
pade approximation
Elapsed time is 1.13181 seconds.
recovery
Elapsed time is 0.195883 seconds.

we could save some time and memory in stages 2 and 3 by rewriting code into C++, so that it can avoid temporary matrices.

There is still a problem with the costly complex schur decomposition. This is currently done simply by converting into complex. By contrast, real schur decomposition is 2x faster, but does not give us a triangular factor.

I don't reckon we can beat the original logm in terms of speed, probably not even if we code it in C++ (although that can get us closer). So it all comes down to speed vs. accuracy issues. Probably those who really use logm may comment.

Jaroslav Hajek <highegg>
Thu 17 Jun 2010 05:16:31 AM UTC, comment #11: 

The proposed new algorithm is awfully slow compared to the old one.  I ran a test with the new algorithm, including new sqrtm, on a random 500x500 matrix.  The results were 13.9 seconds for the new calculation versus 2.8 seconds for the old calculation (5x slower).  Numerically, the two results were within rounding errors of each other.

If the only problem with the old calculation is the occasional problem of inverting the eigenvector matrix, wouldn't it make sense to use the Pade approach only for these matrices and use the fast eigenvalue decomposition for the rest?

In pseudo-code

[V, D] = eig(a);
if (condest (V) < SINGULAR_MATRIX_CUTOFF)
  Eigenvalue approach
else 
  Pade approach
endif

Rik <rik5>
Group administrator
Thu 10 Jun 2010 02:11:45 AM UTC, comment #10: 

I have submitted a revised version of the patch.  The only change I've made is to style in logm.m.

If any style guidelines, or anything else stand in the way of clearing this bug, let me know.

Richard T. Guy

(file #20715)

Anonymous
Wed 09 Jun 2010 07:46:00 AM UTC, comment #9: 

The patch should first be adapted to Octave's coding standards (use specific block ends, 2 spaces indent, parenthesize if-conditions, use # for comments, etc).
Then I think it can be applied and we may test further what needs to be done.

Jaroslav Hajek <highegg>
Sat 05 Jun 2010 01:53:45 PM UTC, comment #8: 

Does anyone know where this bug stands?  It seems like there was fruitful work then things just stagnated.

Is the logm fix going to be pushed into the main build, or are there plans to speed things up farther any maybe add more ancillary features?

Richard T. Guy

Anonymous
Fri 07 May 2010 10:19:46 AM UTC, comment #7: 

I removed the old unconditional error checks because Matlab doesn't seem to do them either:

>> sqrtm([1 -1 -1; 0 1 -1; 0 1 -1])


ans =

   1.0e+08 *

    0.0000    1.0183   -1.0183
         0    0.5092   -0.5092
         0    0.5092   -0.5092


I believe the residual argument is intended to verify the result:

octave:1> [x,err] = sqrtm([1 -1 -1; 0 1 -1; 0 1 -1])
warning: sqrtm: matrix is singular, may not have a square root
x =

  1.0000e+00 + 0.0000e+00i  -3.0000e+00 - 1.3422e+08i  1.0000e+00 + 1.3422e+08i
  0.0000e+00 + 0.0000e+00i  0.0000e+00 - 6.7109e+07i  0.0000e+00 + 6.7109e+07i
  0.0000e+00 + 0.0000e+00i  0.0000e+00 - 6.7109e+07i  0.0000e+00 + 6.7109e+07i

err =  1.9640

isntead of computing the exact residual (which can be costly), Matlab also allows estimating it using a special power method technique. Octave currently can't.


Jaroslav Hajek <highegg>
Thu 06 May 2010 06:31:33 PM UTC, comment #6: 

I just tested the new sqrtm algorithm against my logm code and found that it does help speed.

Old sqrtm:
octave:23> n = rand(500,500);
octave:24> tic, logm(n); toc
warning: Matrix has nonegative eigenvalues.  Principal matrix logarithm is not defined.  Computing non-principal logarithm.
Elapsed time is 94 seconds.


New sqrtm:
octave:29> n = rand(500,500);
octave:30> tic, logm(n); toc
warning: Matrix has nonegative eigenvalues.  Principal matrix logarithm is not defined.  Computing non-principal logarithm.
Elapsed time is 74.21 seconds.

The warning on singular matrices is helpful, but in logm, which uses the sqrtm function potentially many times, the warning message goes to the screen 100 times.  Is there a quick way to fix this or suppress the warning?

Another potential problem with some singular matrices is that sometimes a result is returned with no warning the result is incorrect:
OLD:
octave:32> sqrtm([1 -1 -1; 0 1 -1; 0 1 -1])
error: sqrtm: failed to find sqrt
NEW:

octave:36> sqrtm([1 -1 -1 ; 0 1 -1; 0 1 -1])
ans =

  1.0000e+00 - 0.0000e+00i  -3.0000e+00 - 4.2513e+08i  1.0000e+00 + 4.2513e+08i
  0.0000e+00 - 0.0000e+00i  0.0000e+00 - 2.1256e+08i  0.0000e+00 + 2.1256e+08i
  0.0000e+00 - 0.0000e+00i  0.0000e+00 - 2.1256e+08i  0.0000e+00 + 2.1256e+08i

octave:37> ans*ans
ans =

    1.00000 +  0.00000i   17.39062 -  0.00000i   -6.71875 -  0.00000i
    0.00000 +  0.00000i    6.33594 +  0.00000i   -6.33203 +  0.00000i
    0.00000 +  0.00000i    6.33594 +  0.00000i   -6.33594 +  0.00000i


However, the speedup is good, and this appears to be the right direction.

TG

Anonymous
Thu 06 May 2010 11:55:23 AM UTC, comment #5: 

I checked in the following patch:
http://hg.savannah.gnu.org/hgweb/octave/rev/f9860b622680

Summary of changes:

I wrote an in-place column-major (i.e. cache-aligned) version of the triple loop from Higham's paper for computing the upper triangular sqrtm. The in-place version reduces memory usage when applied to temporaries (such as the result of a real->complex conversion).

For diagonal matrices, sqrtm will simply reduce to sqrt.
For upper triangular matrices, the triple loop will be used directly. Lower triangular matrices are handled via double transpose. If the matrix has a positive diagonal, the sqrtm will be done in real arithmetic.

Here goes a benchmark:

n = 500;

disp ("diagonal matrix");
x = diag (rand (n, 1)); x = x*x;
tic; sqrtm (x); toc

disp ("upper triangular matrix");
x = triu (rand (n) + eye (n)); x = x * x;
tic; sqrtm (x); toc

disp ("lower triangular matrix");
x = tril (rand (n) + eye (n)); x = x * x;
tic; sqrtm (x); toc

disp ("full matrix");
x = rand (n) + eye (n); x = x * x;
tic; sqrtm (x); toc

previously, I got:

diagonal matrix
Elapsed time is 0.562037 seconds.
upper triangular matrix
Elapsed time is 0.554102 seconds.
lower triangular matrix
Elapsed time is 0.577609 seconds.
full matrix
Elapsed time is 1.78512 seconds.

with the new patch, I get:

diagonal matrix
Elapsed time is 0.000648975 seconds.
upper triangular matrix
Elapsed time is 0.0239501 seconds.
lower triangular matrix
Elapsed time is 0.026382 seconds.
full matrix
Elapsed time is 1.54298 seconds.

we still don't have the error estimating code (Matlab has it).

Jaroslav Hajek <highegg>
Thu 06 May 2010 01:00:00 AM UTC, comment #4: 

I haven't had a chance to test the patch put out by David yet, but it appears to be in the right track.

I have two suggestions:
1) Yes, it is good to fix the Frob norm, but this is ancillary to the real problem.
2) The best solution is to have Schur recognize a matrix as triangular, avoid copies, and return the matrix itself as well as eye(n) for the unitary matrix.

Also, for the record, I'm including the text of the email from Dr. Higham:

>Dear Richard,


>Thanks for your email and explanation of what you are proposing >to do.
>Yes - anything allowed by the GPL is fine.


>Regards,
>---Nick Higham


RTG

Anonymous
Wed 05 May 2010 09:46:14 PM UTC, comment #3: 

Here is a test patch that implements some of these ideas with sqrtm.. It hasn't really been tested, but I see speedups only for some matrices (eg randn(n,n)) and not others (eg eye(n))..

Any other ideas?

D.


(file #20455)

David Bateman <dbateman>
Group Member
Tue 04 May 2010 12:17:15 PM UTC, comment #2: 

For 4) there seems to be quite a few things we might do to speed up sqrtm. The frobnorm function seems inefficient and might be rewritten as a template like

Template <class T>
static T
frobnorm (const Array<T>& A)
{
  T sum = T();
  T *fv = A.fortran_vec ();
  for (octave_idx_type i = 0; i < A.numel (); i++)
    sum += real(fv[i] * conj(fv[i]));
  return sqrt (sum);
}

using fortan_vec to avoid the overhead of the checking of the indexing. Similarly the sqrtm_from_schur function might be rewritten to use fortan_vec.

The biggest gain however will come from writing the methods

ComplexMatrix schur (Matrix&, std::string&)

and

FloatComplexMatrix schur (FloatMatrix&, std::string&)

Currently in Fsqrtm. A copy of the real matrix is to a complex matrix is made, and then in the method ComplexSCHUR::init essentially another copy of the complex matrix is made by the copy of write to the variable schur_mat. If the above methods are written the second of these might be avoided.

D.
 

David Bateman <dbateman>
Group Member
Tue 04 May 2010 08:56:53 AM UTC, comment #1: 

Regarding 4), I reckon it will be more Octavish if we just make sqrtm work efficiently for triangular matrices. schur can even utilize the MatrixType system and flag the matrix as triangular already when formed.

More opinions? Do we want the permission from prof. Higham sent to Octave, or is this enough? (I suppose it is).

Jaroslav Hajek <highegg>
Tue 04 May 2010 08:52:07 AM UTC, original submission:  

Octave dev team,

I am submitting a patch to fix the numerically unstable matrix
logarithm function logm.m

1) The original implementation used an eigenvalue decomposition and a
direct inversion of the eigenvector matrix.  This was unstable.  For
instance, the matrix a = [1 -1 -1;0 1 -1; 0 0 1] gave poor results.
2) I utilized a better algorithm based on the Pade approximation.
Code from the mftoolbox by D. Higham was adapted for use.  That code
was published under GPL 3.0 and I contacted him for permission to use
the code.  He agreed to anything within bounds of GPL.
3) I built the difference file with

>> hg diff -U 3

If you require something different, just let me know.  Three files
change: PROJECTS, Changelog, and logm.m.  The last file was completely
rewritten.
4) Only downside to this algorithm is that speed has suffered, but
accuracy trumps speed.  I am thinking about ways to speed it up,
particularly by exploiting a cheaper algorithm to compute sqrtm of an
upper triangular matrix, but the cheap algorithm fails on singular
matrices.

Please let me know if you will accept the patch!

Richard T. Guy
Dept. of Computer Science
Wake Forest University
(soon to be University of Toronto)

Jaroslav Hajek <highegg>

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #21084:  logm.m added by caliari (4KiB - text/x-objcsrc)
file #21080:  rsf2csf_l.m added by caliari (2KiB - text/x-objcsrc)
file #21081:  rsf2csf_g.m added by caliari (2KiB - text/x-objcsrc)
file #21075:  logm.m added by rik5 (3KiB - application/octet-stream)
file #21067:  rsf2csf.m added by caliari (2KiB - text/x-objcsrc)
file #20715:  updated_logm_diff.txt added by None (5KiB - text/plain - updated change for logm.)
file #20455:  changeset.sqrtm added by dbateman (12KiB - application/octet-stream)
file #20419:  octave_logm_patch.txt added by highegg (5KiB - text/plain)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by dbateman (Posted a comment)
  • -email is unavailable- added by highegg
  • -email is unavailable- added by highegg (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 11 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2010-07-28 rik5 StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2010-07-27 caliari Attached File- Added logm.m, #21084
    2010-07-27 caliari Attached File- Added rsf2csf_l.m, #21080
        Attached File- Added rsf2csf_g.m, #21081
    2010-07-26 rik5 Attached File- Added logm.m, #21075
    2010-07-26 caliari Attached File- Added rsf2csf.m, #21067
    2010-06-10 None Attached File- Added updated_logm_diff.txt, #20715
    2010-05-05 dbateman Attached File- Added changeset.sqrtm, #20455
    2010-05-04 highegg Attached File- Added octave_logm_patch.txt, #20419
        Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code