Re: [GSoC2014] Looking for a mentor

[Mogrob Sanit]

On 11/03/2014 22:40, Philip Nienhuis wrote:
> mfasi wrote
>> Hi all,
>>
>> I am a little bit sorry for spamming the maintainers list, but I was
>> suggested to do so on IRC - and I am in good company. As the subject says,
>> I am looking for a mentor to apply for the Google Summer of Code 2014. The
>> topic I am interested in is the implementation of some matrix functions in
>> Octave (project 
> *
>> Improve logm, sqrtm, funm
> *
>> ), as I already have a little experience on the domain. Further
>> information on my 
>> public application draft <http://wiki.octave.org/User:Mfasi>  
>> . Thank you!
> 
> Have you found a mentor yet?
> I regret I am too little proficient in linear algebra to help you out,
> although I wrote the existing funm.m in the OF linear-algebra package. While
> that works OK for my aims (relatively small (n < 100), well conditioned,
> positive semi-definite symmetric matrices) I'd love to see a more robust
> funm.m implementation. But be warned: as Higham wrote in his email to our
> maintainers list, writing a robust funm.m from scratch is no easy task.
> 
> Just some remarks:
> 
> As to logm(), AFAIK that one has already been optimized by "Tommy Guy" some
> years back, along the lines set out by Higham (IIRC Tommy Guy implemented
> some of Higham's methods). I am not aware of things left out (maybe
> complex?).
> 
> As you can see in the OF linear-algebra package, some matrix functions (trig
> and hyperbolic) have been implemented, by invoking expm(), in the thfm.m
> function in that same package. funm.m delegates those trig and hyperbolic
> matrix functions to thfm.m. 
> It would be wise to keep that structure intact.
> 
> Philip
> 
> 
> 
> --
> View this message in context: 
> http://octave.1599824.n4.nabble.com/GSoC2014-Looking-for-a-mentor-tp4662648p4663004.html
> Sent from the Octave - Maintainers mailing list archive at Nabble.com.
> 

Hi Philip,

I have not found a mentor yet, and the fact that you are the only one
that has answered on that topic so far let me think you are the only one
that could be interested. Moreover, the fact of being the author of the
actual funm () makes you one of the most recommended people to mentor
this project. In conclusion, I hope you will change your mind about it. :-D

Talking about serious stuff, I think you are absolutely right about the
behaviour of the current funm () function. The basic weakness, as far as
I can see, is that it uses a special case of the definition by means of
the Jordan decomposition, and I would expect it to be not accurate when
the eigenvectors are very ill-conditioned and to have sometimes a
non-predictable outcome when the Jordan blocks have size bigger than
one, as you pointed out in the source comments. Maybe it will be a
little slow for very big matrices, requiring the computation of the
diagonal decomposition, but I am not sure that this can be completely
eliminated. AFAIK, the best approach for general functions should be the
Schur-Parlett recurrence, that basically computes the matrix function on
the diagonal block of the Schur decomposition and then builds back the
matrix using Sylvester equations. It is robust in general but there are
some issues that have to be addressed to get a good implementations.

On the other hand, according to the documentation, the logm () function
is implemented using some Padé approximants, so I guess it should work
for complex in general, but you are probably right when dealing with
branch points and branch cut. I will try to look at the logm () code to
shed some light on the question.

Best,
Massimiliano