[SciPy-User] raising a matrix to float power

[Charles R Harris]

On Sun, Jul 11, 2010 at 2:17 PM, David Goldsmith <d.l.goldsmith at gmail.com>wrote:

> On Sun, Jul 11, 2010 at 10:41 AM, Charles R Harris <
> charlesr.harris at gmail.com> wrote:
>
>>
>> On Sun, Jul 11, 2010 at 11:31 AM, Alexey Brazhe <brazhe at gmail.com> wrote:
>>
>>> Seems to be, but not for any matrix:
>>>
>>> #-----------------
>>>
>>> def mpower(M, p):
>>>     "Matrix exponentiation"
>>>     e,EV = linalg.eigh(M)
>>>     return dot(EV.T,
>>>                dot(diag((e)**p), EV))
>>>
>>> m = array([[1.0,2.0], [3.0,4.0]])
>>>
>>> then dot(m.T,m) does equal mpower(mpower(dot(m.T,m), 0.5), 2.0)
>>>
>>> But mpower(mpower(m,0.5),2) doesn't equal m!
>>>
>>>
>> For this algorithm the matrix needs to be Hermitean, which is the case for
>> W^T W. More generally, the matrix needs to be normal, i.e., commute with
>> it's transpose. A matrix can be diagonalized iff it is normal.
>>
>
> So it isn't just an algorithmic issue: the general matrix exponentiation
> only works for normal matrices (in theory as well as in computational
> practice)?  Or is it just the (M^a)^(1/a) = (M^(1/a))^a identity that fails
> if M isn't normal?
>

Well, you can square any matrix and the result has a square root. The series
expansion can even converge. For instance [[1,1],[0,1]] isn't normal but has
an n'th root [[1,1/n],[0,1]], i.e.,

In [6]: m = array([[1.,.2],[0.,1.]])

In [7]: dot(dot(dot(dot(m,m),m),m),m)
Out[7]:
array([[ 1.,  1.],
       [ 0.,  1.]])

In fact, all upper triangular matrices with ones on the diagonal have series
with a finite number of terms. So if the matrix is reduced to Jordan form
and none of the diagonal elements are zero, then it should be possible to
find the roots. OTOH,  [[0,1][0,0]] doesn't have a square root. The case of
normal matrices is easy to handle efficiently, however, while reduction to
Jordan form is often difficult and can be numerically tricky.

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100711/388f342a/attachment.html>