from scipy import * # # Example http://www.maths.man.ac.uk/~nareports/narep336.ps.gz # def db(a): # # Denman Neavers # eta = 1.0 Y0 = a Z0 = identity(shape(a)[0]) iter = 0 while eta > 1.e-8: Y1 = 0.5*(Y0+linalg.inv(Z0)) Z1 = 0.5*(Z0+linalg.inv(Y0)) eta = linalg.norm(dot(Y1,Y1)-a) Y0 = Y1 Z0 = Z1 iter = iter + 1 return iter, Y1 def schulz(a): # # Schulz iteration # eta = 1.0 Y0 = a Z0 = identity(shape(a)[0]) iter = 0 while eta > 1.e-8: Y1 = 0.5*dot(Y0,3*identity(shape(a)[0])-dot(Z0,Y0)) Z1 = 0.5*dot(3*identity(shape(a)[0])-dot(Z0,Y0),Z0) eta = linalg.norm(dot(Y1,Y1)-a) Y0 = Y1 Z0 = Z1 iter = iter + 1 return iter, Y1 eps = 5.e-8 eps = 2**-5 a = array(([1,0,0,1],[0,eps,0,0],[0,0,eps,0],[0,0,0,1])) # # matrix sqrt of a # x = array(([1,0,0,0.5],[0,sqrt(eps),0,0],[0,0,sqrt(eps),0],[0,0,0,1])) print dot(x,x) # # matrix sqrt of a by funm # x1 = linalg.funm(a,sqrt) print dot(x1,x1) print linalg.eigvals(a) iterdb, x2 = db(a) print 'Denman Beavers',iterdb, dot(x2,x2) iterschulz, x3 = schulz(a) print 'schulz iteration',iterschulz, dot(x3,x3)