如何用scipy.optimize.minimize求解带行约束与非负性的最小二乘问题

没问题，我来一步步拆解这个问题，帮你用scipy.optimize.minimize实现带约束的最小二乘求解。首先得明确：minimize是针对一维向量优化的，所以我们需要把矩阵T扁平化，再把约束条件转化为向量形式的约束。

问题转化 #

你的原模型是：

P = M*T + error，其中P(n×m), M(n×k), T(k×m)

目标是最小化Frobenius范数平方（等价于最小二乘误差）：
$$\min_{T} |P - M T|_F^2$$

约束条件有两个：

1. T的所有元素非负：$T_{i,j} \geq 0$ 对所有i,j
2. M*T的每一行和为100：$\sum_{j=1}^m (M T)_{i,j} = 100$ 对所有i=1..n

实现步骤 #

1. 导入依赖 #

首先需要导入必要的库：

import numpy as np
from scipy.optimize import minimize

2. 定义目标函数 #

我们把矩阵T扁平化成长度为k*m的一维向量T_vec，在目标函数里再把它reshape回k×m的矩阵。目标函数计算Frobenius范数的平方：

def objective(T_vec, M, P, k, m):
    # 把扁平化的向量转回矩阵T
    T = T_vec.reshape(k, m)
    # 计算误差矩阵的Frobenius范数平方
    error = P - M @ T
    return np.sum(error ** 2)  # 等价于np.linalg.norm(error, 'fro')**2

3. 定义约束条件 #

等式约束：M*T的每一行和为100

我们需要构造一个约束函数，返回M*T的行和 - 100，让这个结果等于0：

def constraint_row_sum(T_vec, M, P, k, m, n):
    T = T_vec.reshape(k, m)
    # 计算M*T的每一行和：M@T 是n×m矩阵，乘以全1的m维向量得到n维行和向量
    row_sums = M @ T @ np.ones(m)
    # 返回行和与100的差值，约束要求这个差值为0
    return row_sums - 100 * np.ones(n)

然后把这个约束包装成minimize能识别的字典形式：

constraints = {
    'type': 'eq',
    'fun': constraint_row_sum,
    'args': (M, P, k, m, n)
}

边界约束：T的所有元素非负

每个元素的下界是0，上界没有限制，用bounds参数：

# 生成k*m个(0, None)的边界，对应T的每个元素
bounds = [(0, None) for _ in range(k*m)]

4. 初始化变量 #

为了让优化收敛更快，建议选择一个合理的初始值。比如可以先计算无约束的最小二乘解，再把负元素设为0：

# 无约束最小二乘解：T_unconstrained = M^+ P，其中M^+是M的伪逆
T_unconstrained = np.linalg.lstsq(M, P, rcond=None)[0]
# 把负元素设为0，作为初始值
T_initial = np.maximum(T_unconstrained, 0)
# 扁平化初始值
T_initial_vec = T_initial.flatten()

或者简单用全1矩阵作为初始值：

T_initial_vec = np.ones(k*m)

5. 调用minimize求解 #

现在把所有参数传入minimize函数，选择合适的优化器。因为有等式约束和边界约束，推荐用SLSQP算法（这是minimize中支持同时处理等式、不等式和边界约束的算法之一）：

# 获取矩阵维度
n, m = P.shape
k = M.shape[1]

# 调用优化
result = minimize(
    fun=objective,
    x0=T_initial_vec,
    args=(M, P, k, m),
    method='SLSQP',
    bounds=bounds,
    constraints=constraints,
    options={'maxiter': 1000, 'disp': True}  # disp=True显示优化过程信息
)

6. 处理结果 #

优化完成后，从结果中提取最优的T矩阵：

if result.success:
    # 把扁平化的最优向量转回矩阵
    T_opt = result.x.reshape(k, m)
    print("最优T矩阵：")
    print(T_opt)
    
    # 验证约束是否满足
    row_sums_opt = M @ T_opt @ np.ones(m)
    print("\nM*T的行和：")
    print(row_sums_opt)
    
    # 计算最优误差
    final_error = np.linalg.norm(P - M @ T_opt, 'fro')**2
    print(f"\n最小二乘误差：{final_error:.4f}")
else:
    print("优化失败，原因：", result.message)

关键注意点 #

• 选择合适的优化器：SLSQP是处理这类带边界和等式约束问题的常用选择，其他比如trust-constr也可以，但SLSQP更轻量。
• 初始值很重要：用无约束解的非负版本作为初始值，能大大加快收敛速度，避免陷入局部最优。
• 约束函数的维度：等式约束函数返回的向量长度必须等于约束的数量（这里是n个，对应n行的行和约束）。

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