[SciPy-User] raising a matrix to float power
Alexey Brazhe
brazhe@gmail....
Sun Jul 11 04:13:09 CDT 2010
Yes, my problem seems to be solved with this:
#--------------------------------
def winvhalf(X):
e, V = linalg.eigh(X)
return dot(V, dot(inv(diag(e+0j)**0.5),V.T))
## I need W = W(W^T W)^{-1/2}
x = dot(B.T,B)
x = winvhalf(x)
B = dot(B, real(x))
#----------------------------------
On Sun, Jul 11, 2010 at 4:24 AM, <josef.pktd@gmail.com> wrote:
> 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
> >
> >
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20100711/31607dd0/attachment.html
More information about the SciPy-User
mailing list