

Hello,
one of the basis of linear algebra is the decomposition of a square matrix into eigenvalues and eigenvectors, such that
[V, D] = eig(A)
A = V*D*inv(V)
This way, each Taylor series expansion can be written as
f(A) = f(A0) + c1 A + c2 A^2 + c3 A^3 + ...
= f(A0) + c1 [V D inv(V)] + c2[ V D inv(V) V D inv(V) ] + c3 [V D inv(V)] ^3
= V[ c0 + c1*D + c2*D^2 + c3*D^3]*inv(V)
= V[f(D)]inv(V)
This, provided that the factorisation on the first line exists, every function can be expanded in a matrix form by applying it elementbyelement on the diagonal matrix D. I've seen there exists an expm and logm functions, is there some generic function accepting a function handle and a square matrix as inputs, performing the factorisation and applying the function handle on the eigenvalues ?
Regards
Pascal


On Tue, 20101221 at 10:34 +0100, CdeMills wrote:
> Hello,
>
> one of the basis of linear algebra is the decomposition of a square matrix
> into eigenvalues and eigenvectors, such that
> [V, D] = eig(A)
> A = V*D*inv(V)
[leads to]
>
> f(A) = V[f(D)]inv(V)
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 'linearalgebra' addon package from OctaveForge does
contain a 'funm' implementation. I haven't used it (or any other such
packages for that matter), but you may want to start from there.
Do note, however, that the eigenvalue expansion might not have the best
cost/benefit ratio. Other implementations I've seen are based on the
Schur decomposition of the matrix rather than the eigenvalue
decomposition. Also, if you can bound the norm(A) you may want to look
into Padé approximations.
Regards,

Bård Skaflestad < [hidden email]>
SINTEF ICT, Applied Mathematics
_______________________________________________
Helpoctave mailing list
[hidden email]
https://wwwold.cae.wisc.edu/mailman/listinfo/helpoctave


Bård Skaflestad wrote
> f(A) = V[f(D)]inv(V)
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 'linearalgebra' addon package from OctaveForge does
contain a 'funm' implementation. I haven't used it (or any other such
packages for that matter), but you may want to start from there.
Do note, however, that the eigenvalue expansion might not have the best
cost/benefit ratio. Other implementations I've seen are based on the
Schur decomposition of the matrix rather than the eigenvalue
decomposition. Also, if you can bound the norm(A) you may want to look
into Padé approximations.
Nice. I also discovered that the linearalgebra package contains a 'thfm' function, which computes trigonometric funcs by using logm and expm, where the latter can be applied to nondiagonalisable matrices.
Pascal


On 21 December 2010 04:05, Bård Skaflestad < [hidden email]> 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 'linearalgebra' addon 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. 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?
 Jordi G. H.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://wwwold.cae.wisc.edu/mailman/listinfo/helpoctave


On Tue, 20101221 at 15:02 +0100, Jordi Gutiérrez Hermoso wrote:
> On 21 December 2010 04:05, Bård Skaflestad < [hidden email]> 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 'linearalgebra' addon 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 SchurParlett algorithm for
computing matrix functions," SIAM J. Matrix Anal. Appl., Vol.
25, Number 2, pp. 464485, 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,
TwentyFive 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 < [hidden email]>
SINTEF ICT, Applied Mathematics
_______________________________________________
Helpoctave mailing list
[hidden email]
https://wwwold.cae.wisc.edu/mailman/listinfo/helpoctave


On Tue, 20101221 at 16:14 +0100, Bård Skaflestad wrote:
> On Tue, 20101221 at 15:02 +0100, Jordi Gutiérrez Hermoso wrote:
> > On 21 December 2010 04:05, Bård Skaflestad < [hidden email]> 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 'linearalgebra' addon 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 followon threads).
Regards,

Bård Skaflestad < [hidden email]>
SINTEF ICT, Applied Mathematics
_______________________________________________
Helpoctave mailing list
[hidden email]
https://wwwold.cae.wisc.edu/mailman/listinfo/helpoctave


Bård Skaflestad wrote
On Tue, 20101221 at 16:14 +0100, Bård Skaflestad wrote:
> On Tue, 20101221 at 15:02 +0100, Jordi Gutiérrez Hermoso wrote:
> > On 21 December 2010 04:05, Bård Skaflestad <Bard.Skaflestad@sintef.no> 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 'linearalgebra' addon 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 followon 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 semidefinite 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 octaveforge 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 email archive.
Several times after that, funm was criticised (for the same good reasons), yet no one came up with a better version or howto 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


On 21 December 2010 16:17, Philip Nienhuis < [hidden email]> wrote:
> 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.
Are we sure N. Higham doesn't retain some copyright to the code? You
seem to know him. Could you ask him if he indeed cannot give us the
actual source code?
> 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 suppose you mean this one:
http://www.maths.manchester.ac.uk/~higham/NAMF/Funm_files.tarI didn't find a copyright statement there. Do you have a suggestion
about how to approach the author over making this code free?
Thanks,
 Jordi G. H.
_______________________________________________
Helpoctave mailing list
[hidden email]
https://wwwold.cae.wisc.edu/mailman/listinfo/helpoctave


Jordi Gutiérrez Hermoso wrote:
> On 21 December 2010 16:17, Philip Nienhuis< [hidden email]> wrote:
>> 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.
>
> Are we sure N. Higham doesn't retain some copyright to the code? You
> seem to know him.
No I don't (apart from knowing his name & some of his reputation).
Actually it was "Tommy Guy" (see thread on helpoctave April 22, 2010)
who seems to have held contacts with Higham and who got his permission
to use logm.m code for Octave.
Just email him  his address is on the web page.
> Thanks,
>  Jordi G. H.
>
Philip
_______________________________________________
Helpoctave mailing list
[hidden email]
https://wwwold.cae.wisc.edu/mailman/listinfo/helpoctave


On 28 December 2010 07:09, Miroslaw Kwasniak
< [hidden email]> wrote:
> What it's unclear from this page?
> http://www.maths.manchester.ac.uk/~higham/NAMF/>
> Matlab Mfiles for the SchurParlett algorithm for computing matrix
> functions are available in the tar file Funm_files.tar. The toolbox
> is distributed under the terms of the GNU General Public License
> (version 3 of the License, or any later version) as published by the
> Free Software Foundation.
That got added yesterday. N Higham responded privately to me about my
request to put up that license text. I'm reproducing his response
below. It looks like we should replace the funm function in
OctaveForge with the implementation in HIgham's page.
 Jordi G. H.
 Forwarded message 
From: Nick Higham < [hidden email]>
Date: 2010/12/27
Subject: Re: Using your (?) funm in Octave
To: Jordi Gutiérrez Hermoso < [hidden email]>
Cc: Edvin Deadman < [hidden email]>
Dear Jordi,
> currently available in Octave. We see that you offer some code on
> your website which you probably mean to distribute freely, but we
> cannot use it unless you give the code an explicit free license
> (Octave favours the GPL version 3, but any other free license of
> your choice is acceptable). Would it be too much of a hassle to
> explicitly put a free license on this?
As you thought, we intended to make the code freely available. I have
updated the Web page with an appropriate licensing statement  is this
OK?
> Additionally, we were wondering if there was any chance of getting
> the actual finished funm code.
The code on the web page was what we produced during the project. The
code in MATLAB is based on refinements of that original code,
including new research results, but I don't have anything
distributable due to funm being copyrighted by The MathWorks. However,
I think it's true to say that the state of the art SchurParlett
algorithm is fully described in the union of my various papers  and
in particular in my book
http://www.maths.manchester.ac.uk/~higham/fmIt's not an easy task to produce an independent high quality funm. One
possibility you can consider is to use what you have currently for
Hermitian matrices and otherwise invoke the "poor man's" version
below, which is based on
http://dx.doi.org/10.1137/060659909The behaviour of this method is not fully understood, but the method
can be justified based on backward error considerations. (Note I have
quickly coded up the Mfile below and not thoroughly tested it.)
Regards,
Nick Higham

function F = funm_davi07(A,fun)
%FUNM_DAVI07 Poor man's version of FUNM based on Davis (2007).
tol = 1e15; % Adjust to taste.
[V,D] = eig(A + tol*randn(size(A)));
F = V*diag(feval(fun,diag(D)))/V;
% For quick text:
% fun = @exp;A = rand(6)^4; F = funm_davi07(A,fun); FX=funm(A,fun);
norm(FFX,1)/norm(FX,1)
_______________________________________________
Helpoctave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/helpoctave

