[SciPy-User] raising a matrix to float power
josef.pktd@gmai...
josef.pktd@gmai...
Sat Jul 10 19:24:39 CDT 2010
On Sat, Jul 10, 2010 at 8:15 PM, Charles R Harris
<charlesr.harris@gmail.com> wrote:
>
>
> On Sat, Jul 10, 2010 at 5:57 PM, Alexey Brazhe <brazhe@gmail.com> wrote:
>>
>> >Sure, M**0.5 is cho_factor(M). For other non-integers I am not sure what
>> >matrix exponentiation could possibly mean.
>>
>> >Are you sure you don't mean array exponentiation?
>>
>> Indeed, I needed to raise a matrix (not array) to power 1/2 (in fact,
>> -1/2).
>> More precisely, I need to compute W(W^TW)^(-1/2).
>> cho_factor fails with "matrix not positive definite", and I don't know how
>> to avoid that
>
> Well, the question remains as to the precise meaning of the square root
> (what is the application?), but my guess is that if you use eigh to
> decompose (W^TW) into u*d*u^T then form u*(d^{-1/2}*u^T you will get what
> you need. Maybe ;) Zero elements of d, if any, will be a problem.
Part of some code I used where I don't find the cleaned up version
omega = np.dot(dummyall, dummyall.T)
ev, evec = np.linalg.eigh(omega) #eig doesn't work
omegainvhalf = evec/np.sqrt(ev)
print np.max(np.abs(np.dot(omegainvhalf,omegainvhalf.T) - omegainv))
# now we can use omegainvhalf in GLS (instead of the cholesky)
no guarantees,
Josef
>
> You can also use the svd if the previous interpretation is correct, since if
> W = u*d*v the whole expression above reduces to u*d^(.5)*v.
>
> Chuck
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
More information about the SciPy-User
mailing list