我正在寻找解决 python 中的最小二乘问题，以便我最小化 0.5*||np.dot(A,x) - b||^2 受限于约束 np.dot(C,x) = d，边界为 0<x<np.inf，其中
A : shape=[m, n]
C : shape=[m, n]
b : shape=[m]
d : shape=[m]
都是已知的矩阵。
不幸的是，scipy.optimize.lsq_linear() 函数似乎只适用于上限/下限约束：
minimize 0.5 * ||A x - b||**2
subject to lb <= x <= ub
不是等式约束。此外，我想为解决方案添加界限，使其仅是正数。有没有使用scipy 或NumPy 内置函数的简单或干净的方法？
【问题讨论】：
标签： python numpy optimization scipy least-squares
如果你坚持scipy.optimize.lsq_linear，你可以在目标函数中添加一个惩罚项，即
minimize 0.5 * ||A x - b||**2 + beta*||Cx-d||**2
subject to lb <= x <= ub
并选择足够大的标量beta，以便将您的约束优化问题转换为具有简单框约束的问题。但是，对于约束优化问题，使用scipy.optimize.minimize 更方便：
import numpy as np
from scipy.optimize import minimize
# ... Define your np.ndarrays A, C, b, d here
# constraint
cons = [{'type': 'eq', 'fun': lambda x: C @ x - d}]
# variable bounds
bnds = [(0, None) for _ in range(A.shape[1])]
# initial point
x0 = np.ones(A.shape[1])
# res.x contains your solution
res = minimize(lambda x: np.linalg.norm(A@x-b)**2, x0=x0, bounds=bnds, constraints=cons)
我使用了欧几里得范数。用m = 2000 和n = 1000 为玩具示例计时，在我的机器上会产生 3 分钟和 10 秒。
请注意，目标梯度和约束雅可比都由有限差分近似，这需要对目标和约束函数进行多次评估。因此，对于大问题，代码可能太慢了。在这种情况下，值得提供精确的梯度，雅可比和黑森。此外，由于
# f(x) = ||Ax-b||^2 = x'A'Ax - 2b'Ax + b'b
# grad f(x) = 2A'Ax - 2A'b
# hess f(x) = 2A'A
通过预先计算A'A、b'A 和b'b，我们可以显着减少评估函数所需的时间。另外，我们只需要计算一次A'Ax：
from scipy.optimize import NonlinearConstraint
# Precalculate the matrix-vector products
AtA = A.T @ A
btA = b.T @ A
btb = b.T @ b
Atb = btA.T
# return the objective value and the objective gradient
def obj_and_grad(x):
AtAx = AtA @ x
obj = x.T @ AtAx - 2*btA @ x + btb
grad = 2*AtAx - 2*Atb
return obj, grad
#constraint d <= C@x <= d (including the jacobian and hessian)
con = NonlinearConstraint(lambda x: C @ x, d, d, jac=lambda x: C, hess=lambda x, v: v[0] * np.zeros(x.size))
# res.x contains your solution
res = minimize(obj_and_grad, jac=True, hess=lambda x: AtA2, x0=x0, bounds=bnds, constraints=(con,), method="trust-constr")
这里，选项jac=True 告诉minimize 我们的函数obj_and_grad 返回一个包含目标的元组
函数和梯度。我们进一步选择了 'trust-constr' 方法，因为它是唯一支持 hessian 的可用方法。为上述同一个玩具示例再次计时，在我的机器上产生 7 秒。
【讨论】：