In [16]:

    

import math
import numpy

S = numpy.matrix([
  [ 2.263E-01, 1.991E-01, 3.086E-01, 4.942E-01, 3.409E-01],
  [ 1.991E-01, 7.264E-01, 6.258E-01, 5.923E-01, 5.284E-01],
  [ 3.086E-01, 6.258E-01, 7.201E-01, 7.445E-01, 4.964E-01],
  [ 4.942E-01, 5.923E-01, 7.445E-01, 1.252E+00, 8.759E-01],
  [ 3.409E-01, 5.284E-01, 4.964E-01, 8.759E-01, 8.980E-01]
  ])

N = S.shape[0]
MAX_ITERATION = 20
THRESHOLD = 1e-10

Z = numpy.matrix(numpy.identity(N))
Y = S

lam = 1/numpy.linalg.norm(S, 'fro')
print("||S||_{F} = %e, lambda = %e" % (numpy.linalg.norm(S, 'fro'), lam))

for i in range(MAX_ITERATION):
    
    #print("Z_%d:" % (i+1))
    #print(Z)
    #print(numpy.linalg.norm(Z, 'fro'))

    #print("Y%d:" % (i+1))
    #print(Y)
    #print(numpy.linalg.norm(Y, 'fro'))
    #print(numpy.linalg.norm(Y[0:4,0:4], 'fro'))
    #print(numpy.linalg.norm(Y[0:4,4:8], 'fro'))
    #print(numpy.linalg.norm(Y[4:8,0:4], 'fro'))
    #print(numpy.linalg.norm(Y[4:8,4:8], 'fro'))

    X = lam*Y*Z
    #print("X_%d" % (i+1))
    #print(X)
    #print(numpy.linalg.norm(X, 'fro'))
    
    e = numpy.linalg.norm(X-numpy.identity(N), 'fro')
    print("error(%2d) = %e" % (i+1, e))
    if e < THRESHOLD:
        break
        
    T = 1/2*(3*numpy.identity(N)-X) # Second order.
    #T = 1/8*(15*numpy.identity(N)-10*X+3*X*X) # Third order.
    #T = 1/16*(35*numpy.identity(N)-35*X+21*X*X-5*X*X*X) # Fourth order.
    #print("T_%d" % (i+1))
    #print(T)
    #print(numpy.linalg.norm(T, 'fro'))
    #print(numpy.linalg.norm(T[0:4,0:4], 'fro'))
    #print(numpy.linalg.norm(T[0:4,4:8], 'fro'))

    Z = Z*T
    Y = T*Y
    
X = math.sqrt(lam)*Z
XT = math.sqrt(lam)*Y # Eq. (26) uses 1/sqrt(lam), which seems wrong.
print("X (S^{-1/2}):")
print(X)
#print("X^{-1} (S^{1/2}):")
#print(X.I)
print("XT (S^{1/2}:")
print(XT)


    

    




||S||_{F} = 3.105005e+00, lambda = 3.220607e-01
error( 1) = 1.880865e+00
error( 2) = 1.760351e+00
error( 3) = 1.571237e+00
error( 4) = 1.357474e+00
error( 5) = 1.173267e+00
error( 6) = 9.859161e-01
error( 7) = 7.892269e-01
error( 8) = 5.668276e-01
error( 9) = 2.861967e-01
error(10) = 6.729184e-02
error(11) = 3.472321e-03
error(12) = 9.053228e-06
error(13) = 6.147116e-11
X (S^{-1/2}):
[[ 8.95768031  2.75087178 -3.25182993 -1.71753856 -1.3528134 ]
 [ 2.75087178  3.49657629 -2.74524671 -0.09339462 -1.20494078]
 [-3.25182993 -2.74524671  4.43552359 -0.52724687  1.08002201]
 [-1.71753856 -0.09339462 -0.52724687  2.26993536 -0.73883184]
 [-1.3528134  -1.20494078  1.08002201 -0.73883184  2.2050112 ]]
XT (S^{1/2}:
[[ 0.26132527  0.01458755  0.19394894  0.29995205  0.17380703]
 [ 0.01458755  0.67762785  0.39255941  0.21432909  0.25878146]
 [ 0.19394894  0.39255941  0.61611805  0.35506152  0.15070095]
 [ 0.29995205  0.21432909  0.35506152  0.89815577  0.42818143]
 [ 0.17380703  0.25878146  0.15070095  0.42818143  0.77121538]]

In [ ]: