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

编程问答 2022-08-15

问题描述

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

问题的一些背景：

我有像素值为原始数字（DN）形式的图像，并且我想拿出一条回归线来模拟DN与图像中表面的真实反射率值之间的线性关系。 / p>

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

import numpy as np
A = np.array([[1,19039],[1,47792],9672],32521],11409],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×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

Result with MatLab

但是我使用了optimize.minimize函数

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

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

Result with Python

我认为可能是我设置的约束存在问题，但是我尝试解决这些约束，但并没有以任何方式改善结果。

在Python中是否有更好的方法来解决此问题？还是我目前的方法可以改进？

谢谢。

解决方法

感谢sascha在评论中的建议，我设法找到了解决方案：

import cvxpy as cp
import numpy as np

A = np.array([[1,19039],[1,47792],9672],32521],11409],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[:,np.multiply(A[:,slope) + intercept)

哪个会根据我的约束给我最合适的线

enter image description here

如果我检查并比较原始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,x.value[1]) + x.value[0]) - b[i]
    c_res.append(residual)
    residual = (np.multiply(A[i,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
constraintsleast-squaresoptimizationoptimizationpythonscipyscipy