Re: sqrt() of a Matrix in a DLD

[Paul Kienzle]

Use sqrtm, expm, logm, funm, thfm if you want matrix functions.

Note: sqrtm, expm are fast and stable. funm uses
the diagonalization you suggest, which is fast and
unstable.   thfm implements the trig and hyperbolic
functions using a combination of sqrtm, expm, inv
and logm, so will perform as well as they do on
ill-conditioned cases.

We have code which handles logm for ill-conditioned
cases as well which I believe has been posted but
not yet incorporated into octave-forge/octave.  Many
thanks to Ross Lippert for all the work he put into it.
I was waiting until I could improve the performance
of the well-conditioned cases, but clearly if I haven't
found the time in the last three years it is time to put
it in as is and let somebody else take over.

Paul Kienzle
address@hidden

On Mar 7, 2004, at 3:44 PM, Vic Norton wrote:

> Something is very peculiar about this, John. What does sqrt(m) mean 
> anyhow?
>
> For example suppose
>
>    m = [1 -3; 0 2];
>
> then, according to octave,
>
>    sm = sqrt(m) =
>
>    1.00000 + 0.00000i  0.00000 + 1.73205i
>    0.00000 + 0.00000i  1.41421 + 0.00000i
>
> The complex elements are there alright, but
>
>    sm^2 =
>
>    1.00000 + 0.00000i  0.00000 + 4.18154i
>    0.00000 + 0.00000i  2.00000 + 0.00000i
>
> is certainly not m.
>
> In fact m is diagonalizable with positive diagonal elements
>
>    m = inv(t) * d * t ,  t = [1 3; 0 1],  d = diag([1 2]).
>
> It follows that
>
>    rm = inv(t) * sqrt(d) * t =
>
>    1.00000  -1.24264
>    0.00000   1.41421
>
> does have the property that  rm^2 = m.
>
> IMHO, rm is the square root of m.
>
>
> Regards,
>
> Vic
>
>
> At 1:39 PM -0600 3/5/04, John W. Eaton wrote:
> > For example, if
> > you write
> >
> >   sm = sqrt (m);
> >
> > then sm will be complex if any element of m is negative.  Do you want
> > to preserve that behavior in your C++ function?
> >
> > jwe
>
>
>
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
>
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------