使用Gradient Descent最小化欧几里得范数求解线性方程组出错求助

我尝试用Python的梯度下降法最小化$|Ax-b|^2$来求解以下线性方程组：

• $x - 2y + 3z = -1$
• $3x + 2y - 5z = 3$
• $2x - 5y + 2z = 0$

方程组的正确解为$x=0.2608$，$y=-0.0869$，$z=-0.4782$，但我通过梯度下降得到的结果不正确：

• $x=0.6397305590418562$
• $y=0.24242515764463315$
• $z=-0.4736842052386$

我的代码如下：

A = np.array([[1, -2, 3], [3, 2, -5], [2, -5, 2]])
X = np.array([[x, y, z]])
b = np.array([[-1], [3], [0]])

f = (A*X - b).transpose() * (A*X - b)
display(f)
print('---')

# Define the function and derivatives
g = (x + 1)**2 + (1 - 2*y)**2 + (3*z + 1)**2 + (3*x - 3)**2 + (2*y - 3)**2 + (-5*z - 3)**2 + 4*x**2 + 25*y**2 + 4*z**2
g = sp.expand(g)
print(g)
print('---')
g_x = sp.diff(g, x)
g_y = sp.diff(g, y)
g_z = sp.diff(g, z)
print(g_x)
print(g_y)
print(g_z)
print('---')

# Define the NumPy function
def g_func(x, y, z):
    return 14*x**2 - 16*x + 33*y**2 - 16*y + 38*z**2 + 36*z + 30

def grad_g(x, y, z):
    return 28*x - 16, 66*y - 16, 76*z + 36

# Running the Gradient Descent Method
cur_pos = (10, 10, 10)

def gradient_descent(f, grad, x0, alpha = 0.01, max_iter = 1000, tol = 1e-6):
    k = 0
    while k < max_iter and abs(grad(x0[0], x0[1], x0[2])[0]) > tol and abs(grad(x0[0], x0[1], x0[2])[1]) > tol and abs(grad(x0[0], x0[1], x0[2])[2]) > tol:
        x_der, y_der, z_der = grad(x0[0], x0[1], x0[2])
        x0 = x0[0] - alpha * x_der, x0[1] - alpha * y_der, x0[2] - alpha * z_der
        k += 1
    return x0[0], x0[1], x0[2], f(x0[0], x0[1], x0[2]), grad(x0[0], x0[1], x0[2])

gradient_descent(g_func, grad_g, cur_pos)

错误排查与修正 #

核心错误：目标函数手动展开完全错误 #

$|Ax-b|^2$的正确展开是$(Ax - b)^T(Ax - b)$，展开后会包含变量的交叉项（如$xy$、$xz$、$yz$），但你错误地将每个残差项拆成了单个变量的平方和（比如把第一个残差$x-2y+3z+1$拆成$(x+1)^2+(1-2y)2+(3z+1)^2$），丢失了所有交叉项，导致目标函数和梯度完全偏离正确形式。

次要错误 #

1. 矩阵运算错误：原代码中A*X是元素-wise乘法，矩阵乘法应使用@符号（或np.dot）；
2. 学习率设置不合理：原代码的alpha=0.01对于正确的目标函数来说过大，容易导致迭代震荡，无法收敛到最优解；
3. 符号变量未定义：代码中直接使用x,y,z但未通过sympy定义符号变量，属于语法遗漏。

修正后的代码 #

import numpy as np
import sympy as sp

# 定义符号变量
x, y, z = sp.symbols('x y z')

A = np.array([[1, -2, 3], [3, 2, -5], [2, -5, 2]])
X = sp.Matrix([x, y, z])
b = sp.Matrix([-1, 3, 0])

# 正确计算目标函数 ∥Ax - b∥²
residual = A @ X - b
g = residual.T @ residual
g = sp.expand(g[0])  # 转置后为1x1矩阵，提取数值元素
print("正确展开的目标函数：")
print(g)

# 计算梯度分量
g_x = sp.diff(g, x)
g_y = sp.diff(g, y)
g_z = sp.diff(g, z)
print("\n梯度分量：")
print(f"g_x: {g_x}")
print(f"g_y: {g_y}")
print(f"g_z: {g_z}")

# 转换为NumPy可调用的函数
def g_func(x_val, y_val, z_val):
    return float(g.subs({x: x_val, y: y_val, z: z_val}))

def grad_g(x_val, y_val, z_val):
    dx = float(g_x.subs({x: x_val, y: y_val, z: z_val}))
    dy = float(g_y.subs({x: x_val, y: y_val, z: z_val}))
    dz = float(g_z.subs({x: x_val, y: y_val, z: z_val}))
    return dx, dy, dz

# 梯度下降迭代函数
def gradient_descent(f, grad, x0, alpha=0.001, max_iter=10000, tol=1e-6):
    x_cur, y_cur, z_cur = x0
    k = 0
    while k < max_iter:
        dx, dy, dz = grad(x_cur, y_cur, z_cur)
        # 判断梯度是否足够小，满足收敛条件
        if abs(dx) < tol and abs(dy) < tol and abs(dz) < tol:
            break
        # 更新参数
        x_cur -= alpha * dx
        y_cur -= alpha * dy
        z_cur -= alpha * dz
        k += 1
    return x_cur, y_cur, z_cur, f(x_cur, y_cur, z_cur), grad(x_cur, y_cur, z_cur)

# 初始点与迭代
cur_pos = (10, 10, 10)
result = gradient_descent(g_func, grad_g, cur_pos)
print("\n梯度下降结果：")
print(f"x: {result[0]:.4f}, y: {result[1]:.4f}, z: {result[2]:.4f}")

运行修正后的代码，会得到与正确解一致的结果：$x≈0.2609$，$y≈-0.0870$，$z≈-0.4783$。

内容的提问来源于stack exchange，提问作者Orhan94