Quantcast

matrix functions

classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

matrix functions

CdeMills
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 element-by-element 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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: matrix functions

Bård Skaflestad
On Tue, 2010-12-21 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 'linear-algebra' add-on 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

_______________________________________________
Help-octave mailing list
[hidden email]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: matrix functions

CdeMills
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 'linear-algebra' add-on 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 linear-algebra package contains a 'thfm' function, which computes trigonometric funcs by using logm and expm, where the latter can be applied to non-diagonalisable matrices.

Pascal
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: matrix functions

Jordi Gutiérrez Hermoso
In reply to this post by Bård Skaflestad
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 '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. 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.

_______________________________________________
Help-octave mailing list
[hidden email]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

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 <[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 '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 <[hidden email]>
SINTEF ICT, Applied Mathematics

_______________________________________________
Help-octave mailing list
[hidden email]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: matrix functions

Bård Skaflestad
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 <[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 '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).


Regards,
--
Bård Skaflestad <[hidden email]>
SINTEF ICT, Applied Mathematics

_______________________________________________
Help-octave mailing list
[hidden email]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

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 <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 '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
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: matrix functions

Jordi Gutiérrez Hermoso
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.tar

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

_______________________________________________
Help-octave mailing list
[hidden email]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: matrix functions

Philip Nienhuis
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 help-octave April 22, 2010)
who seems to have held contacts with Higham and who got his permission
to use logm.m code for Octave.

>                    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.tar
>
> I didn't find a copyright statement there. Do you have a suggestion
> about how to approach the author over making this code free?

Just email him - his address is on the web page.

> Thanks,
> - Jordi G. H.
>

Philip
_______________________________________________
Help-octave mailing list
[hidden email]
https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: matrix functions

Miroslaw Kwasniak
In reply to this post by Jordi Gutiérrez Hermoso
On Wed, Dec 22, 2010 at 08:02:23AM -0600, Jordi Gutiérrez Hermoso wrote:
> I suppose you mean this one:
>
>      http://www.maths.manchester.ac.uk/~higham/NAMF/Funm_files.tar
>
> I didn't find a copyright statement there. Do you have a suggestion
> about how to approach the author over making this code free?

??
What it's unclear from this page?
http://www.maths.manchester.ac.uk/~higham/NAMF/

  Matlab M-files for the Schur-Parlett 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.
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: matrix functions

Jordi Gutiérrez Hermoso
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 M-files for the Schur-Parlett 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
Octave-Forge 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 Schur-Parlett
algorithm is fully described in the union of my various papers - and
in particular in my book

     http://www.maths.manchester.ac.uk/~higham/fm

It'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/060659909

The 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 M-file 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 = 1e-15;  % 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(F-FX,1)/norm(FX,1)
_______________________________________________
Help-octave mailing list
[hidden email]
https://mailman.cae.wisc.edu/listinfo/help-octave
Loading...