How to implement a constrained linear fit in Python?

Question

I'm trying to fit a linear model to a set of data, with the constraint that all the residuals (model - data) are positive - in other words, the model should be the "best overestimate". Without this constraint, linear models can be easily found with numpy's polyfit as shown below.

import numpy as np
import matplotlib.pyplot as plt

x = [-4.12179107e-01, -1.40664082e-01, -5.52301563e-06,  1.82898473e-01]
y = [-4.14846251, -3.31607886, -3.57827245, -5.09914559]

plt.scatter(x,y)
coeff = np.polyfit(x,y,1)
plt.plot(x,np.polyval(coeff,x),c='r',label='numpy-polyval')
plt.plot(x,np.polyval([-2,-3.6],x),c='g',label='desired-fit') #a rough guess of the desired result
plt.legend()

example1

Is there an efficient way to implement a linear fit with this type of constraint?

user_15187728 · Answer

This is a quadratic programming problem. There are several libraries (CVXOPT, quadprog etc.) that can be used to solve it. Here is an example using quadprog:

import numpy as np
import matplotlib.pyplot as plt
import quadprog

x = [-4.12179107e-01, -1.40664082e-01, -5.52301563e-06, 1.82898473e-01]
y = [-4.14846251, -3.31607886, -3.57827245, -5.09914559]

A = np.c_[x, np.ones(len(x))]
y = np.array(y)
G = A.T @ A
a = A.T @ y
C = A.T
b = y
coeffs = quadprog.solve_qp(G, a, C, b)[0]

plt.scatter(x, y)
plt.plot(x, np.polyval(coeffs, x), c='r')
plt.show()

This gives:

See e.g. this post for more information. In describes, in particular, how to set up a linear regression problem as a quadratic programming problem.

As a side note, the optimal line will always pass through one data point, but it need not pass through two such points. For example, take x = [-1., 0., 1.] and y = [1., 2., 1.].

user_1883316 · Answer

Yes, the best fit is a line through the top two points. I do an argsort to find the top Ys, compute the slope and y-intercept, and off we go:

import numpy as np
import matplotlib.pyplot as plt

x = [-4.12179107e-01, -1.40664082e-01, -5.52301563e-06,  1.82898473e-01]
y = [-4.14846251, -3.31607886, -3.57827245, -5.09914559]

plt.scatter(x,y)
coeff = np.polyfit(x,y,1)
model1 = np.polyval(coeff,x)
model1 += (y-model1).max()
print(model1)
print(sum((y-model1)**2))

z = np.argsort(y)
pt0 = (x[z[-1]],y[z[-1]])
pt1 = (x[z[-2]],y[z[-2]])
m = (pt1[1]-pt0[1])/(pt1[0]-pt0[0])
b = pt0[1]-m*pt0[0]

model2 = np.polyval([m,b],x)
print(model2)
print(sum((y-model2)**2))

plt.plot(x,model1,c='r',label='numpy-polyval')
plt.plot(x,model2,c='g',label='generated')
plt.legend()
plt.show()

Output:

How to implement a constrained linear fit in Python?

2 Answers