[SciPy-user] sqrtm

[Nils Wagner]

Hi all, 
 
Just another sqrtm implementation based on Pade iteration... 
 
def sqrtm2(A,p): 
# 
# Pade Iteration 
# 
  m,n = shape(A) 
  i=arange(1,p+1) 
  xi = 0.5*(1.+cos(0.5*(2*i-1)*pi/p)) 
  ai2 = 1.0/xi-1. 
  Y0 = A 
  Z0 = identity(n) 
  maxiter = 5 
  for k in arange(0,maxiter): 
 
    sum1= zeros((n,n),Complex) 
    sum2= zeros((n,n),Complex) 
    for j in arange(0,p): 
      sum1 = sum1 + linalg.inv(dot(Z0,Y0)+ai2[j]*identity(n))/xi[j] 
      sum2 = sum2 + linalg.inv(dot(Y0,Z0)+ai2[j]*identity(n))/xi[j] 
    Y1 = dot(Y0,sum1)/p 
    Z1 = dot(Z0,sum1)/p 
    Y0 = Y1 
    Z0 = Z1 
  return Y1,Z1 
 
# 
# test matrix 
# 
A=rand(5,5) 
# 
# Modify the spectrum 
# 
A = dot(A,A)/25.0+identity(5) 
# 
Y2,Z2 = sqrtm2(A,4) 
print linalg.norm(dot(Y2,Y2)-A) 
print linalg.norm(dot(Z2,Z2)-linalg.inv(A)) 
 
Any suggestions ? 
 
Nils