[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: sqrt() of a Matrix in a DLD
From: |
Paul Kienzle |
Subject: |
Re: sqrt() of a Matrix in a DLD |
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
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
-------------------------------------------------------------