Re: matrix functions

[Philip Nienhuis]

Bård Skaflestad wrote:
> 
> On Tue, 2010-12-21 at 16:14 +0100, Bård Skaflestad wrote:
>> On Tue, 2010-12-21 at 15:02 +0100, Jordi Gutiérrez Hermoso wrote:
>> > On 21 December 2010 04:05, Bård Skaflestad <address@hidden>
>> wrote:
>> > 
>> > > I was about to suggest you use the 'funm' function, but now I see
>> that
>> > > 'funm' is only available (by default) in That Other Software.  On the
>> > > other hand, the 'linear-algebra' add-on package from OctaveForge does
>> > > contain a 'funm' implementation.
>> > 
>> > Huh, that's awkward, because the docstring for expm mentions funm,
>> > which isn't included in core. I see the funm implementation is quite
>> > naïve and relies on the matrix being diagonalisable.
>> 
>> Naïve indeed.  The implementation uses exactly the algorithm alluded to
>> by the OP and, as you say, implements even this algorithm poorly.  I
>> almost feel like withdrawing my original reference to this particular
>> function...
> 
> Oh, lest I forget:  This has been extensively discussed in the mailing
> list before.  See, e.g., the thread entitled "logm robustness" from
> April of this year (and a few similarly entitled follow-on threads).
> 

Being the author of that -indeed- naive funm I may be entitled to say some
more about it. 
Not so much for defense BTW, I do agree that funm.m as it stands is nothing
more than a kludge.
But a kludge that works sufficiently reliable for my own needs (symmetric
positive semi-definite tridiagonal matrices of limited order (< 50), no
repeated eigenvalues, ...). I wrote it > 10 yrs ago for certain problems
where "naive was good enough".

At that time, Paul Kienzle helped me improve it quite a bit. Soon after
inclusion in octave-forge Rolf Fabian criticised the robustness and Paul
Kienzle made a few improved versions that he shared with me. Due to lack of
time he and I gave up soon, and none of PK's improved versions ever made it
to svn (or whatever was in place at the time). I've got them archived
somewhere in a stack of decommissioned hard disks, plus probably in an old
e-mail archive.

Several times after that, funm was criticised (for the same good reasons),
yet no one came up with a better version or how-to suggestions for
improvement. Oh well; as I said, it works OK for my needs...

Prof. N. Higham developed a much improved funm version for Matlab, I suppose
funded by TMW. That is, the header of ML's R2010b funm.m (I made sure I
didn't peek at the code proper) shows it to be (C) TMW.
So, anyone in need for a better funm.m probably needs to start from scratch.
Higham's papers and the theses of his PhD students should provide a good
start. Much of it seems to be available through his site (see below).
For a start, some more elaborate funm versions can be obtained from Higham's
site, AFAICS under GPL (one has to search a bit around there):
  http://www.maths.manchester.ac.uk/~higham/NAMF/

I do hope someone will take this up. Linear algebra is not quite my
specialty and Higham's work shows that developing a sufficiently reliable
funm is not an easy job (to wit: a quick scroll down shows that TMW's
version contains > 650 lines of code).

Philip

-- 
View this message in context: 
http://octave.1599824.n4.nabble.com/matrix-functions-tp3137935p3159793.html
Sent from the Octave - General mailing list archive at Nabble.com.