Re: matrix functions

[Bård Skaflestad]

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...

>  Even worse, it
> calls inv instead of using the backslash operator. Is this the best
> that can be done in general? Surely if expm works for defective
> matrices, something can be done in general if we assume analytic
> functions, which funm is assuming anyways?

There are at least two classical references in this area

      * Davies, P. I. and N. J. Higham, "A Schur-Parlett algorithm for
        computing matrix functions," SIAM J. Matrix Anal. Appl., Vol.
        25, Number 2, pp. 464-485, 2003.
      * Golub, G. H. and C. F. Van Loan, Matrix Computation, Third
        Edition, Johns Hopkins University Press, 1996, p. 384.

And, of course, I would be remiss if I failed to mention Moler and Van
Loan's "Nineteen Dubious Ways to Compute the Exponential of a Matrix"
and "Nineteen Dubious Ways to Compute the Exponential of a Matrix,
Twenty-Five Years Later", although these papers are mostly about the
exponential.

That said, implementing some of these more advanced algorithms is not
entirely trivial.  A bit of work and analysis is needed...


Regards,
-- 
Bård Skaflestad <address@hidden>
SINTEF ICT, Applied Mathematics