在 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)

我将 objective 函数设置为 (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

感谢 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

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

Solving Least Squares with Linear Inequality Constraints in Python

python

optimization

constraints

scipy

least-squares