# Finding least upper bound line above a set of points using NMinimize

I have this set of points:

sTest={{265., -12.4317}, {269., -12.4839}, {270., -12.6383}, {274., \
-12.6947}, {278., -12.6996}, {279., -12.7755}, {300., -12.8652}, \
{309., -12.9276}, {313., -13.0571}, {314., -13.0811}, {315., \
-13.0985}, {316., -13.1501}, {322., -13.1908}, {323., -13.2129}, \
{324., -13.2311}, {325., -13.2757}, {326., -13.3849}, {332., \
-13.7497}, {335., -13.762}, {355., -13.7981}, {359., -13.9145}, \
{364., -14.0233}, {368., -14.0583}, {369., -14.1957}, {371., \
-14.2555}, {377., -14.4896}, {378., -14.5427}, {390., -14.5961}, \
{399., -14.6112}, {412., -14.7939}, {414., -14.805}, {415., \
-14.9441}, {423., -15.1411}, {424., -15.2414}, {445., -15.256}, \
{454., -15.3825}, {458., -15.4697}, {459., -15.6437}, {460., \
-15.6559}, {461., -15.7129}, {467., -15.8281}, {468., -15.9638}, \
{470., -15.987}, {489., -16.295}, {500., -16.3044}, {501., -16.3501}, \
{504., -16.3555}, {513., -16.4882}, {514., -16.7162}, {544., \
-16.7981}, {547., -16.9239}, {548., -16.9534}}


And I initially draw a line from the first to last point:

And notice there are some points above and below the line. I which to find the line such that all points are below the line and the sum of the vertical distance from each point to the line is a minimum. I can manually, find an approximate minimum using Manipulate:

which is the line $$-7.52-0.017x$$ giving a value of the sum of $$5.77$$. I would think I could use NMinimize to find the minimum line even if I have to set a limit on the values of $$a$$ and $$b$$ via:

sumF[a_, b_, buf_] :=
Sum[(a + b buf[[i, 1]]) - buf[[i, 2]], {i, 1, Length@buf}]

checkBuf[a_, b_, buf_] := Module[{},
VectorQ[Table[(a + b buf[[i, 1]]) - buf[[i, 2]],
{i, 1, Length@buf}], Positive]
]

NMinimize[{sumF[a, b, sTest], -8 < a < -6 && -1 < b < 0 &&
checkBuf[a, b, sTest] == True}, {a, b}]


where checkBuf returns True if all points are below the line. However, this gives the error:

NMinimize::nsol: There are no points that satisfy the constraints {False}.

I was wondering if this is not a problem NMinimize can solve or if it is, if I'm not coding it correctly and if so could someone help me do so?

• The appears to be a strong periodicity in the residuals: ListPlot[Transpose[{sTest[[All, 1]], LinearModelFit[sTest, x, x]["FitResiduals"]}], Joined -> True]. Is there a rationale for ignoring that?
– JimB
Commented Oct 25, 2022 at 16:42
• @JimB: I see what you mean. This data plots the log of the accuracy as a function of number of terms for a fractional power expansion $\sum_{n=1}^N c_i (z^{1/9})^n$. I assume the periodicity stems from the exponential terms. I'm interested in generating a simple linear expression $A(N)=a+bN$ which gives the minimum accuracy for $N$ terms (hense having all the points below the line.
– josh
Commented Oct 25, 2022 at 17:35

We may define a line: c0+ c1 x with parameters c0 and c1. To determine the parameters we request that the following sum is minimized:

Total[((c0  + c1 #[[1]]) - #[[2]]) & /@ sTest];


under the conditions that the line is above the points:

 And @@ (((c0  + c1 #[[1]]) >= #[[2]]) & /@ sTest)


With this we can determine the parameters c0 and c1:

line = c0 + c1 x /.
Minimize[{Total[((c0  + c1 #[[1]]) - #[[2]]) & /@ sTest],
And @@ (((c0  + c1 #[[1]]) >= #[[2]]) & /@ sTest)}, {c0, c1}][[2]]
Plot[line, {x, sTest[[1, 1]], sTest[[-1, 1]]}, Epilog -> Point[sTest]]


• Outstanding Daniel! Let me check it a bit on some data sets and get back to you later.
– josh
Commented Oct 25, 2022 at 22:03
• Works well for my application over varied data sets. Also have some fitted to quadratics nicely! Thanks!
– josh
Commented Oct 29, 2022 at 11:59
• Glad to be of some help. Commented Oct 29, 2022 at 14:48

It seems plausible that a line that connecting two points on the convex hull is the solution.

(* Find points in convex hull *)
chr = ConvexHullRegion[sTest][[1]];

(* Find the two points on the convex hull border that satisfy the constraints and
minimize the total error *)
minError = ∞;
besta = 0;
bestb = 0;
(* Loop through all possible pairs of points on the convex hull *)
Do[Do[
(* Slope *)
bb = (chr[[i, 2]] - chr[[j, 2]])/(chr[[i, 1]] - chr[[j, 1]]);
(* Intercept *)
aa = chr[[i, 2]] - bb chr[[i, 1]];

If[Min[aa + bb chr[[All, 1]] - chr[[All, 2]]] == 0,
error =
aa Length[sTest] + bb Total[sTest[[All, 1]]] -
Total[sTest[[All, 2]]];
If[minError == ∞, minError = error; besta = aa;
bestb = bb];
If[minError > error, minError = error; besta = aa; bestb = bb]],

{j, 1, i - 1}], {i, 2, Length[chr]}]

(* Show results *)
{besta, bestb}
(* {-7.8383, -0.0164702} *)

ListPlot[{sTest, {{Min[sTest[[All, 1]]],
besta + bestb Min[sTest[[All, 1]]]}, {Max[sTest[[All, 1]]],
besta + bestb Max[sTest[[All, 1]]]}}}, PlotStyle -> {Blue, Red},
Joined -> {False, True}]


Another alternative is to use NonlinearModelFit:

sTest = Rationalize[sTest, 0];
nlm = NonlinearModelFit[sTest, {a + b x, Min[a + b sTest[[All, 1]] - sTest[[All, 2]]] > 0},
{{a, -8}, {b, -1/100}}, x, WorkingPrecision -> 20]
nlm["BestFitParameters"]
Show[ListPlot[sTest], Plot[nlm[x], {x, Min[sTest[[All, 1]]], Max[sTest[[All, 1]]]},
PlotStyle -> Red]]


• Thanks Jim! I suspected a convex hull might be applicable. Will check it with some larger data sets and get back to you.
– josh
Commented Oct 25, 2022 at 22:04

A slightly generalized supplement of @DanielHuber's nice answer, which describes the line in vector form {p1,p2}+t {e1,e2} (includes possible case of infinte slope)

e = {e1, e2};
p = {p1, p2};


This approach minimizes the projection perpendicular to the unknown direction e and is therefore independent of the underlying xy-coordinate system

\[CapitalDelta] = Map[Cross[e] . (# - p) &, sTest];
mini = NMinimize[{\[CapitalDelta] . \[CapitalDelta],Map[# < 0 &, \[CapitalDelta]], e . e == 1}, {e1, e2, p1, p2}]

Show[{ListPlot[sTest],Graphics[InfiniteLine[{p, p + e}] /. mini[[2]]]   }]


A variation of Daniel Huber's answer is to use matrix-form LinearOptimization directly. Speed doesn't matter on the example case much, but this is 2.7 times faster on my system, and with ten million points, this grows to a tenfold improvement:

line = {1, x} . LinearOptimization[
{Length[sTest], Total[sTest[[All, 1]]]},
Transpose[{{1, #[[1]]}, -#[[2]]} & /@ sTest]];
Plot[line, {x, sTest[[1, 1]], sTest[[-1, 1]]}, Epilog -> Point[sTest]]