0

I'm running into an optimization issue. The problem has a dual objective. I have 2 clouds of points, the green ones $(x_i, y_i)$ and the red ones $(x_j, y_j)$. I'm looking for the equations of 2 lines, one being above the green points and the other being below the red points, with one constraint: they have to be parallel. As the objective function, I just summed the distances of each line to its respective cloud of points, the distance being just a difference (the constraints will ensure that they meet the required criteria). Here is the formulation of the problem:

$$\min_{a,b_u,b_l} {\sum_{i=1}^m {\left( \left( ax_i+b_u \right) - y_i \right)} + \sum_{j=1}^n {\left( y_j - \left( ax_j+b_l \right) \right)}} \ \ \ \text{s.t.}\ \left\{\matrix {ax_i+b_u\ge y_i \ \ \forall \ i \in [1,m] \\y_j\ge ax_j+b_l \ \ \forall \ j \in [1,n]}\right.$$

To simplify the notations, the indices $i$ and $j$ refer to the clouds (i.e. $x_i\ne x_j$ and $y_i\ne y_j$).

The program can be put in standard form: $$\min_{a,b_u,b_l} {a\left(\sum_{i=1}^m x_i - \sum_{j=1}^n x_j\right) + b_u m - b_l n} \ \ \ \text{s.t.}\ \left\{\matrix {-ax_i-b_u\le -y_i \ \ \forall \ i \in [1,m] \\ax_j+b_l\le y_j \ \ \forall \ j \in [1,n]}\right.$$

I implemented it with scipy's linprog:

import numpy as np
from scipy.optimize import linprog


def fit_both_lines(x_u, y_u, x_l, y_l):
m = len(x_u)
n = len(x_l)
c = np.array([np.sum(x_u) - np.sum(x_l), m, -n], dtype=float)


A_ub = np.array([(-x_u[i], -1, 0) for i in range(m)] + [(x_l[j], 0, 1) for j in range(n)])
b_ub = np.array([-y_u[i] for i in range(m)] + [y_l[j] for j in range(n)])

res = linprog(c, A_ub, b_ub, method="highs-ds")

slope, intercept_u, intercept_l = res.x

return (intercept_u, slope), (intercept_l, slope)

Here is a snippet to test it:

from numpy.random import randint, random
import matplotlib.pyplot as plt


N = 10
x_u = randint(0, 100, N)
y_u = random(N) * 100


x_l = randint(0, 100, N)
y_l = random(N) * 100


upper_line, lower_line = fit_both_lines(x_u, y_u, x_l, y_l)


plt.scatter(x_u, y_u, c="green")
plt.scatter(x_l, y_l, c="red")
plt.plot(x_u, upper_line[0] + upper_line[1] * x_u, c="green")
plt.plot(x_l, lower_line[0] + lower_line[1] * x_l, c="red")

Though, most of the time, the slope $a$ gets estimated to 0 which makes me think there is a mistake somewhere.

Here is a sample plot:

enter image description here

I also tried the absolute values as a distance function and linearized them, but I got the very same result... What did I do wrong?

Mark Morrisson
  • 123
  • I would have first tested on something I know the result such as red points $(0,0)$ and $(1,1)$ and green points $(0,-1)$ and $(1,0)$. Also it seems you use $m=n$ in your test. Try with $m\neq n$. I am not sure to understand what your + in the construction of A_ub and b_ub does in that case. – Marc Dinh Oct 17 '23 at 08:03
  • Well, I tested of course on fewer points and it fails. I've used more points on the picture for illustration purposes. In general, m and n will be different, but it won't pose any problem, it will just change the number of constraints (hence the dimensions of A_ub and b_ub). I've redone the calculations many times and I can't find where the problem is... – Mark Morrisson Oct 17 '23 at 12:25
  • If you are sure that A_ub and b_ub is correct, then it may comes from the default options of linprog. Do you know what are the default options for the bounds with linprog? I suggest you enforce no bounds for them to be sure that these default options do not intefere. – Marc Dinh Oct 17 '23 at 16:04
  • By default, the variables are unbounded in linprog. I'm not sure the optimization program is correct. If I fit only one line at a time (either red or green), it works fine, but if I try to fit both at the same time, it fails. Either there's a computation of coding mistake somewhere, or the problem is wrongly modelled. – Mark Morrisson Oct 17 '23 at 16:18
  • Ok, then can you provide the value of A_ub and b_ub you get with your program with the data I rpovided above? – Marc Dinh Oct 17 '23 at 16:30
  • Ah also, the formulation is correct to me. – Marc Dinh Oct 17 '23 at 17:49
  • With the set of points you suggested, the solution returned by the solver is two identical lines with intercept 0 and slope 0. But that is incorrect: I expect the red line to go through the 2 red points, and the green line to go through the 2 green points. Interesting case... – Mark Morrisson Oct 18 '23 at 15:49
  • Then you can check step by step. Your program is short enough to do that. With those points, write down the mathematical formulation you should get and compare to what you get from your program. – Marc Dinh Oct 19 '23 at 08:02
  • Yeah, I did that yesterday, and I found that (0, 0, 0) returned by linprog was not the optimal minimum for the objcetive function (0 vs -2 for the 2 parallel lines), so I suspect a problem with linprog. I tried with the "simplex" method, and the result was different... Very strange. – Mark Morrisson Oct 19 '23 at 08:22
  • 1
    Ahhhh, by checking scipy's reported issues, it seems that by default the decision variables are non-negative. Though, this is not very clear in the documentation. By specifying bounds=(None, None), it works as expected! – Mark Morrisson Oct 19 '23 at 12:05

0 Answers0