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 |
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 |
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 |
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 |
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 |
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 |
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 |
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 |
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 |
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 |
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 |
Free forum by Nabble | Edit this page |