# matrix functions

11 messages
Open this post in threaded view
|
Report Content as Inappropriate

## matrix functions

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

## Re: matrix functions

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

## Re: matrix functions

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

## Re: matrix functions

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

## Re: matrix functions

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

## Re: matrix functions

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

## Re: matrix functions

Open this post in threaded view
|
Report Content as Inappropriate

## Re: matrix functions

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

## Re: matrix functions

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