[Scilab-Dev] Problem with matrix logarithms, in the complex case

[Serge Steer]

joao moreira a écrit :
> Hello,
>
> I tried posting this on the users list, but got no answers, so I'll try
> here. Basically, I'm trying to understand why the matrix logarithm is
> implemented the way it is (in macros/elem/logm.sci, in scilab 4.1.2),
> and if
> there is some reason for it not calculating in the case described below.
>
> --------
>
> Hi all,
>
> I have an issue with matrix logarithms, in the complex case. Consider
> the following matrix :
>
> A=[
>     -0.01+0*%i 0.13+0*%i;
>      0.00+0*%i 0.08+0*%i;
>   ];
>
> Coefficients are real, but I force them to be complex, hoping that
> scilab will use a complex logarithm matrix function. But scilab says
> "unable to diagonalize!"... I understand that as a real matrix, I
> can't calculate its log, but as a complex matrix, why doesn't this
> work ?
>
> Using algorithms from Golub, van Loan, "Matrix Computations", I coded
> a complex matrix log function, and found logm(A) = L where
>
> L=[
>     -4.6051702+3.1415927*%i  3.0036378-4.5378561*%i;
>      0.0000000+0.0000000*%i -2.5257286+0.0000000*%i;
>   ];
>
> and it's easy to check that exp(log(A)) == log(exp(A)) == A.
>
> I looked in the scilab source, and found that logm is actually a .sci
> file that does
> .........................
> a=a+0*%i;   //Set complex
> [s,u,bs]=bdiag(a);
>   if maxi(bs)>1 then
>     error('logm: unable to diagonalize!');return
>   end
> .........................
>
> bdiag does block diagonalisation, finds a 2x2 block, and refuses to
> continue !... Is this normal ? is there some reason for scilab not
> calculating this ?
>
> Maybe I'm missing something obvious, thanks for helping.

No you are right, there is a problem with your matrix. We uses bdiag in
logm to ensure that the a matrix is diagonalizable. But with your matrix
bdiag with default rmax find that a is not diagonalizable. The default
rmax is too strong for this use.
It should be better to do something like

rmax=max(norm(a,1),1/sqrt(%eps))
[s,u,bs]=bdiag(a,rmax);

I modified the Scilab code this way. But I think a more clever algorithm
as to developped to be able to compute logm for quasi non diagonalizable
matrices...



Thanks for the report
Serge
>
> Joao
>
>
> -- 
> Joao Moreira (math student)
> Paris, France
> joao.mrr1 at gmail.com <mailto:joao.mrr1 at gmail.com>
>