在 Python 中求解具有线性不等式约束的最小二乘

Question

我正在尝试解决受 Python 中不等式约束的线性系统影响的最小二乘问题。我已经能够在 MatLab 中解决这个问题，但是对于我正在工作的项目，我们所有的代码库都应该是 Python，所以我正在寻找一种等效的方法来解决它，但一直无法解决。

问题的一些背景：

我的图像具有原始数字 (DN) 形式的像素值，我想提出一条回归线，用于模拟 DN 与图像中表面的真实反射率值之间的线性关系。

我有 6 个已知的反射率和对应的 DN，因此我建立了一个线性方程组：

import numpy as np
A = np.array([[1, 19039],[1, 47792], [1, 9672], [1, 32521], [1, 11409], [1, 58843]])
b = np.array([[0.05938044], [0.27213514], [0.00252875], [0.18535543], [0.01959069], [0.52605937]])
b = np.squeeze(b)

我将目标函数设置为 (Ax-b)*0.5 的 2-范数

def f(, x):
    # function to minimize
    y = np.dot(A, x) - b
    return (np.dot(y, y))*.5

然后我想添加我的约束。由于表面反射率值介于 0-1 之间，我想确保我的图像中的最小 DN 值在通过估计斜率和截距系数从 DN 转换为反射率时具有大于 0 的反射率值，并且最大 DN 值图像的反射率被映射到低于或等于 1。

根据我正在实现的paper，我可以将0 <= slope*DN + intercept <= 1的需求分成两部分：

slope*DN_max + intercept <= 1

-slope*DN_min - intercept <= 0

因此，我创建了两个矩阵，一个包含最小 DN 值和截距系数 (1)，另一个包含最大 DN 值和截距系数 (1)。我将它们制成矩阵，因为在实践中我将校准不止一个图像，因此我将有超过两列和两行（我将有 n x n*2 矩阵，其中 n = 图像数），但是对于这个简化的示例，我将只使用一张图片。

img_min = 0
img_max = 65536
C_min = np.array([[-1, img_min]])
C_max = np.array([[1, img_max]])

def cons_min(x):
    # constraint function for the minimum pixel values
    return np.array((C_min @ x))

def cons_max(x):
    # constraint function for the maximum pixel values
    return np.array((C_max @ x))-1

然后我尝试使用 optimize.minimize 来求解系数。

con1 = [{'type': 'ineq', 'fun': cons_min},
        {'type': 'ineq', 'fun': cons_max}
       ]
initial = [0, 0]
result = optimize.minimize(f, 
                           initial, 
                           method='SLSQP', 
                           constraints=con1
                           )

在 MatLab 中使用 lsqlin(A,B, C, c) 函数得到的结果是

intercept = 0.0000000043711483450798073327817072630808
slope = 0.0000069505505573854644717126521902273

但是使用我得到的 optimize.minimize 函数

[-2.80380803e-17,  1.52590219e-05]

列表的第一个值是截距，第二个是斜率。

我认为这可能是我设置的约束有问题，但我尝试过使用它们，但并没有以任何方式改善结果。

有没有更好的方法在 Python 中解决这个问题？或者我目前的方法有什么可以改进的地方吗？

谢谢。

Answer 1

感谢 cmets 中 sascha 的建议，我终于找到了解决方案：

import cvxpy as cp
import numpy as np

A = np.array([[1, 19039],[1, 47792], [1, 9672], [1, 32521], [1, 11409], [1, 58843]])
b = np.array([[0.05938044], [0.27213514], [0.00252875], [0.18535543], [0.01959069], [0.52605937]])
b = np.squeeze(b)
C_min = np.array([[-1, 0]])
C_max = np.array([[1, 65535]])

x = cp.Variable(A.shape[1])

objective = cp.Minimize(0.5 * cp.sum_squares(A@x-b))
constraints = [C_min@x <= 0, C_max@x <= 1]

prob = cp.Problem(objective, constraints)
result = prob.solve(solver=cp.ECOS)

intercept = x.value[0]
slope = x.value[1]

import matplotlib.pyplot as plt
%matplotlib inline

plt.figure()
plt.scatter(A[:, 1], b)
plt.plot(A[:, 1], np.multiply(A[:, 1], slope) + intercept)

根据我的限制，这给了我最适合的线

如果我检查并比较原始 MatLab 解决方案和 cvxpy 解决方案之间的残差，我发现在这个示例中 cvxpy 解决方案稍微好一点（尽管非常小）。

# MatLab estimated values for the slope and intercept
ML_inter = 0.0000000043711483450798073327817072630808
ML_slope = 0.0000069505505573854644717126521902273

# get the residuals for each data point
c_res = []
ml_res = []
for i in range(A.shape[0]):
    residual = (np.multiply(A[i, 1], x.value[1]) + x.value[0]) - b[i]
    c_res.append(residual)
    residual = (np.multiply(A[i, 1], ML_slope) + ML_inter) - b[i]
    ml_res.append(residual)

# calculate the sum of squares
ss_cvx = np.sum(np.array(c_res)**2)
ss_ml = np.sum(np.array(ml_res)**2)

print("Sum of squares for cvx:    ", ss_cvx)
print("Sum of squares for matlab: ", ss_ml)
print("Sum of squares is lower for CVX solution? ", ss_cvx < ss_ml)

# Sum of squares for cvx:     0.03203817995131549
# Sum of squares for matlab:  0.032038181467959566
# Sum of squares is lower for CVX solution?  True