def jacobi_iteration(A,b,init_x,iters=10):
dim = init_x.shape[0]
x = init_x
result = np.zeros((iters+1,dim))
result[0,:] = x.reshape(-1)
diag_A = np.diag(A)
D_inv = np.diag(1.0/diag_A)
LU = np.diag(diag_A) - A
B = D_inv @ LU
D_inv_b = D_inv @ b
for i in range(iters):
x = B @ x + D_inv_b
result[i+1,:] = x.reshape(-1)
return result
%/accordion%
Gauss-Seidel
(然而并没什么卵用,不必Jacobi快多少,还不能并行计算了。。。)
%accordion% Numpy实现 %accordion%
def gauss_seidel_iteration(A,b,init_x,iters=10):
dim = init_x.shape[0]
x = init_x
result = np.zeros((iters+1,dim))
result[0,:] = x.reshape(-1)
diag_A = np.diag(A)
D_inv = np.diag(1.0/diag_A)
LU = np.diag(diag_A) - A
B = D_inv @ LU
D_inv_b = D_inv @ b
for i in range(iters):
for d in range(dim):
x[d] = B[d,:] @ x + D_inv_b[d]
result[i+1,:] = x.reshape(-1)
return result
def gradient_desc(A,b,init_x,iters=10):
dim = init_x.shape[0]
x = init_x
result = np.zeros((iters+1,dim))
result[0,:] = x.reshape(-1)
for i in range(iters):
ax = A @ x
grad = ax - b
flat_x = x.reshape(-1)
lr = np.dot(flat_x,flat_x)/np.dot(ax.reshape(-1),flat_x)
x -= grad * lr
result[i+1,:] = x.reshape(-1)
return result