From help-octave-request at bevo dot che dot wisc dot edu Sun Mar  7 15:17:19 2004
Subject: Re: sqrt() of a Matrix in a DLD
From: Paul Kienzle <pkienzle at users dot sourceforge dot net>
To: Octave_post Help <help-octave at bevo dot che dot wisc dot edu>
Date: Sun, 7 Mar 2004 16:17:05 -0500

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
pkienzle at users dot sf dot net

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
-------------------------------------------------------------