使用Scipy求解带约束线性方程组的异常问题求助

需要求解带约束的线性方程组 (AX = b)，其中：

• (A) 是8×8分块对角矩阵
• (X = [x1, y1, x2, y2, x3, y3, x4, y4]^T)，(b) 是8×1向量
• 约束条件：
  • (0.5(y1 + y2) = 0)
  • (0.5(x3 + x4) = 0)
  • (0.5(y3 + y4) = 0)

尝试过两种方法均未得到正确结果：

• 将约束整合为11×8矩阵，用伪逆求解，解不满足约束
• 改用Scipy的SLSQP优化器构造最小二乘问题，求解器显示成功且解满足约束，但代入 (AX) 后结果与 (b) 偏差极大，接近全零向量

原问题代码 #

import numpy as np
from scipy import optimize

A = np.array([
    [-261.60,  11.26,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [   4.07, -12.75,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,  -158.63,  -5.65,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,    -2.81, -12.14,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0, -265.99,  19.29,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,   12.59, -12.34,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0, -166.25, -12.63],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0,   -8.40, -11.14]
])

b = np.array([
     -6.95,
     16.35,
     -0.96,
     16.35,
     19.19,
    -15.85,
    -12.36,
    -15.63]).reshape(-1, 1)

def objective_function(x):
    return np.sum((np.dot(A, x) - b)**2)

def constraints(x):
    return np.array([
        0.5*(x[1] + x[3]),  # 0.5*(y1 + y2) = 0
        0.5*(x[4] + x[6]),  # 0.5*(x3 + x4) = 0
        0.5*(x[5] + x[7])   # 0.5*(y3 + y4) = 0
    ])

cons = {'type': 'eq', 'fun': constraints}
options = {'disp': True, 'verbose': 4}

res = optimize.minimize(objective_function, np.zeros(8), method='SLSQP', constraints=cons, options={'disp': True})

x_optimized = res.x

# Checking if solution was found
print(np.matmul(A, x_optimized))

原代码输出结果 #

[ 0.01723619  0.00041891  0.01781721 -0.00034223  0.0113301  -0.00848171
  0.01008367  0.00781168]

问题根源分析 #

核心问题是矩阵维度不匹配导致目标函数计算错误：

• (b) 被reshape为(8,1)列向量，但优化输入的 (x) 是一维数组(8,)，计算 (A@x) 得到一维数组(8,)，与列向量 (b) 做减法时，numpy会广播为(8,8)矩阵，平方求和后得到的是所有元素的平方和，完全偏离了原本要最小化的 (|Ax - b|_2^2) 目标。

解决方案 #

修正后的SLSQP优化代码 #

import numpy as np
from scipy import optimize

A = np.array([
    [-261.60,  11.26,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [   4.07, -12.75,      0.0,    0.0,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,  -158.63,  -5.65,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,    -2.81, -12.14,     0.0,    0.0,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0, -265.99,  19.29,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,   12.59, -12.34,     0.0,    0.0],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0, -166.25, -12.63],
    [    0.0,    0.0,      0.0,    0.0,     0.0,    0.0,   -8.40, -11.14]
])

# 保持b为一维数组，与Ax输出维度匹配
b = np.array([
     -6.95,
     16.35,
     -0.96,
     16.35,
     19.19,
    -15.85,
    -12.36,
    -15.63
])

def objective_function(x):
    # 计算真实的残差平方和
    return np.sum((A @ x - b)**2)

def constraints(x):
    # 去掉冗余的0.5，约束等价且计算更高效
    return np.array([
        x[1] + x[3],  # y1 + y2 = 0
        x[4] + x[6],  # x3 + x4 = 0
        x[5] + x[7]   # y3 + y4 = 0
    ])

cons = {'type': 'eq', 'fun': constraints}
# 增加迭代次数确保收敛充分
options = {'disp': True, 'maxiter': 1000}

res = optimize.minimize(objective_function, np.zeros(8), method='SLSQP', constraints=cons, options=options)

x_optimized = res.x

# 验证结果
print("优化得到的X：", x_optimized)
print("Ax的结果：", A @ x_optimized)
print("与b的误差：", A @ x_optimized - b)

更高效的带约束线性最小二乘方法 #

使用scipy.optimize.lsq_linear专门求解带约束的线性最小二乘问题，效率更高：

from scipy.optimize import lsq_linear

# 构造约束矩阵C和向量d，满足Cx = d
C = np.array([
    [0, 1, 0, 1, 0, 0, 0, 0],
    [0, 0, 0, 0, 1, 0, 1, 0],
    [0, 0, 0, 0, 0, 1, 0, 1]
])
d = np.array([0, 0, 0])

res_lsq = lsq_linear(A, b, bounds=(-np.inf, np.inf), constraints={'type': 'eq', 'fun': lambda x: C @ x - d})
print("lsq_linear得到的X：", res_lsq.x)
print("Ax的结果：", A @ res_lsq.x)

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