Detailed Explanation of Python Implementation of Jacobi Iterative Algorithm

import numpy as np
import time

1.1 Jacobi iterative algorithm

def Jacobi_tensor_V2(A,b,Delta,m,n,M):
start=time.perf_counter()# Start timing 
find=0# Used to mark whether it converges within the specified number of steps 
X=np.ones(n)# Iteration starting point 
x=np.ones(n)# Used to store the intermediate results of iterations 
d=np.ones(n)# Used for storage Ax**(m-2) Diagonal part of 
m1=m-1
m2=2-m
for i in range(M):
print('X',X)
a=np.copy(A)
# Get Ax**(m-2)
for j in range(m-2):
a=np.dot(a,X)
# Get d  And  (2-m)Dx**(m-2)+(L'+U')x**(m-2)
for j in range(n):
d[j]=a[j,j]
a[j,j]=m2*a[j,j]
# Iterative update 
for j in range(n):
x[j]=(b[j]-np.dot(a[j],X))/(m1*d[j])
# Judge whether the accuracy requirements are met 
if np.max(np.fabs(X-x))<Delta:
find=1
break 
X=np.copy(x)
end=time.perf_counter()# End timing 
print(' Time: ',end-start)
print(' Iteration ',i)
return X,find,i,end-start

1.2 Generating functions for tensor A and vector b:

def Creat_A(m,n):# Generating tensor A
size=np.full(m, n)
X=np.ones(n)
while 1:
# Randomly generate tensors of a given shape A
A=np.random.randint(-49,50,size=size)
# Judge Dx**(m-2) Whether it is not singular, if it is, it meets the requirements and jumps out of the loop 
D=np.copy(A)
for i1 in range(n):
for i2 in range(n):
if i1!=i2:
D[i1,i2]=0
for i in range(m-2):
D=np.dot(D,X)
det=np.linalg.det(D)
if det!=0:
break
# Will A Diagonal plane tensor expansion of 10 Times, so that the diagonal surface is dominant 
for i1 in range(n):
for i2 in range(n):
if i1==i2:
A[i1,i2]=A[i1,i2]*10
print('A:')
print(A)
return A
# By A And the given X According to Ax**(m-1)=b Generating vector b
def Creat_b(A,X,m):
a=np.copy(A)
for i in range(m-1):
a=np.dot(a,X)
print('b:')
print(a)
return a

1.3 Generating function of symmetric tensor S:

def Creat_S(m,n):# Generate symmetric tensor B
size=np.full(m, n)
S=np.zeros(size)
print('S',S)
for i in range(4):
# Generate n Is a vector a
a=np.random.random(n)*np.random.randint(-5,6)
b=np.copy(a)
# Right a Go on m-1 Second outer product, get rank 1 Symmetric tensor b
for j in range(m-1):
b=outer(b,a)
# Will be different b Low rank symmetric tensor is obtained by superposition S
S=S+b
print('S:')
print(S)
return S
def outer(a,b):
c=[]
for i in b:
c.append(i*a)
return np.array(c)
return a

1.4 Experiment 1

def test_1():
Delta=0.01# Precision 
m=3#A Order of 
n=3#A Dimension of 
M=200# Maximum number of iteration steps 
X_real=np.array( [2,3,4])
A=Creat_A(m,n) 
b=Creat_b(A,X_real,m)
Jacobi_tensor_V2(A,b,Delta,m,n)