102
edits
Mfasi
Joined 28 February 2014
→Your task
| Line 64: | Line 64: | ||
Basically, I would like to provide Octave with an (almost) state of the art library to compute the most diffused matrix functions in a robust way. Octave's standard library as well as the linear algebra package (Octave Forge) already implement some of these functions, notably funm (), expm (), logm() and sqrtm(). Some of these implementations are good, while others should be improved or, in some cases, completely reimplemented. | Basically, I would like to provide Octave with an (almost) state of the art library to compute the most diffused matrix functions in a robust way. Octave's standard library as well as the linear algebra package (Octave Forge) already implement some of these functions, notably funm (), expm (), logm() and sqrtm(). Some of these implementations are good, while others should be improved or, in some cases, completely reimplemented. | ||
One of the functions that should be completely rewritten is probably funm (), that is currently implemented using a simple but in some cases weak approach inspired by the definition of matrix function via the Jordan canonical form. This approach fails when the matrix is not diagonalizable and can be very inaccurate when the eigenvectors are ill-conditioned. To avoid these issues I would like to use the Schur-Parlett recurrence with block reordering, a well known algorithm that though having some issues that should be addressed, is in theory very robust. This | One of the functions that should be completely rewritten is probably funm (), that is currently implemented using a simple but in some cases weak approach inspired by the definition of matrix function via the Jordan canonical form. This approach fails when the matrix is not diagonalizable and can be very inaccurate when the eigenvectors are ill-conditioned. To avoid these issues I would like to use the Schur-Parlett recurrence with block reordering, a well known algorithm that though having some issues that should be addressed, is in theory very robust. This implementation will require to: | ||
* modify the Schur function to deal with specific reordering of the eigenvalues; | * modify the Schur function to deal with a specific reordering of the eigenvalues; | ||
* implement a function that evaluates matrix polynomials using the Paterson-Stockmeyer scheme; | * implement a function that evaluates matrix polynomials using the Paterson-Stockmeyer scheme; | ||
* implement a function to find a block pattern (how to group the eigenvalues of the Schur form); | * implement a function to find a block pattern (how to group the eigenvalues of the Schur form); | ||
* implement a function to find a confluent permutation (of course, this is a very specific task). | * implement a function to find a confluent permutation (of course, this is a very specific task). | ||
The following functions that should be addressed are expm (), logm () and sqrtm (). The implementation of the first function is stable, using Padé approximants, but could be improved using an heuristic that switches to some other kind of evaluation when such rational functions are ill-conditioned. The logm () functions is a recent implementation and uses a scaling and squaring method | The following functions that should be addressed are expm (), logm () and sqrtm (). The implementation of the first function is stable, using Padé approximants, but could be improved using an heuristic that switches to some other kind of evaluation when such rational functions are ill-conditioned. The logm () functions is a recent implementation and uses a scaling and squaring method. Maybe some little improvements can be done, but I would not put them between the priorities. I would not touch sqrtm () either, since it seems to be a good C++ implementation of an algorithm by Higham that should be very similar to the one in Matlab. I would instead focus on developing the p-th root, as there is a [http://poisson.phc.unipi.it/~maxreen/bruno/pdf/B.%20Iannazzo%20and%20C.%20Manasse%20-%20A%20Schur%20logarithmic%20algorithm%20for%20fractional%20powers%20of%20matrices%20-%20SIMAX.pdf recent work] that seems promising, though being rather involved could be tricky to implement in a robust way. | ||
The hyperbolic and trigonometric functions have already been implemented straightforwardly in thfm (), from the linear-algebra package. As they rely on sqrtm (), expm () and logm (), they will benefit of all the improvements to such functions. | |||
==== Tentative timeline ==== | ==== Tentative timeline ==== | ||