Mfasi

User:Mfasi (edit)

Revision as of 22:35, 16 March 2014

1,178 bytes added , 16 March 2014

→‎Your task

Mfasi

102

edits

@@ Line 62: / Line 62: @@
 * The project I would like to work on is '''Improve logm, sqrtm, funm'''.
-I would like to write such functions, using known algorithms for matrix functions. In particular, I would like to start implementing the Schur-Parlett recurrence, that should not be a hard task as syl () is already there, and then proceed implementing better - specific - algorithms for other relevant matrix functions. In particular, I think that the following list (ordered by relevance) should be respected:
+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.
-* funm
-* logm [just a check, see below]
+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 task will require to:
-* signm
-* sqrtm
+* modify the Schur function to deal with specific reordering of the eigenvalues;
-* rootm
+* implement a function that evaluates matrix polynomials using the Paterson-Stockmeyer scheme;
-* hyperbolic and trigonometric matrix functions
+* implement a function to find a block pattern (how to group the eigenvalues of the Schur form);
-Some other auxiliary functions should be required and implemented, possibly in C++ for performance's sake. It seems that something has already been done for the logm () function, see [http://octave.1599824.n4.nabble.com/GSoC2014-Looking-for-a-mentor-td4662648.html] and [http://octave.1599824.n4.nabble.com/logm-robustness-td2016796.html], while the implementation of funm () lacks of robustness as it uses a simple diagonalization scheme that should be avoided when dealing with general matrices. For the p-th root 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 but being rather involved could be tricky to implement in a robust way.
+* 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 with Parlett recurrence
+It seems that something has already been done for the logm () function, see [http://octave.1599824.n4.nabble.com/GSoC2014-Looking-for-a-mentor-td4662648.html] and [http://octave.1599824.n4.nabble.com/logm-robustness-td2016796.html], while the implementation of funm () lacks of robustness as it uses a simple diagonalization scheme that should be avoided when dealing with general matrices. For the p-th root 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 but being rather involved could be tricky to implement in a robust way.
 ==== Tentative timeline ====