[SciPy-User] raising a matrix to float power

[David Goldsmith]

On Sun, Jul 11, 2010 at 3:02 PM, Charles R Harris <charlesr.harris at gmail.com
> wrote:

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

OK, so it is, at least in part, an algorithmic issue; which begs the
question: does scipy have a more generally applicable, though still robust
(i.e., behaves nicely when it fails) pow(matrix, float) function?

DG

>
> Chuck
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>


-- 
Mathematician: noun, someone who disavows certainty when their uncertainty
set is non-empty, even if that set has measure zero.

Hope: noun, that delusive spirit which escaped Pandora's jar and, with her
lies, prevents mankind from committing a general suicide.  (As interpreted
by Robert Graves)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20100711/0f9cb4d0/attachment.html>