Underestimation of f(x) by using a piecewise linear function

Question

I am trying to check if there is any Matlab/Python procedure to underestimate f(x) by using a piecewise linear function g(x). That is g(x) needs to be less or equal to, f(x). See the picture and code below. Could you please help to modify this code to find how to underestimate this function?

 x = 0.000000001:0.001:1;
 y = abs(f(x));

 %# Find section sizes, by using an inverse of the approximation of the derivative
 numOfSections = 5;
 totalRange = max(x(:))-min(x(:));

 %# The relevant nodes
 xNodes = x(1) + [ 0 cumsum(sectionSize)];
 yNodes = abs(f(xNodes));

 figure;plot(x,y);
 hold on;
 plot (xNodes,yNodes,'r');
 scatter (xNodes,yNodes,'r');
 legend('abs(f(x))','adaptive linear interpolation');

Why is that question you linked relevant? What do you mean by "underestimates" and "overestimates"? You want a function `g(x)` which always has a value `g(x) < f(x)`? Why not simply define `g(x) = f(x) - 1` (or any other constant)? You have just plotted two plots with different resolutions. In both cases, both your y coordinates are identical for a given x coordinate — Pranav Hosangadi, Jul 27 '22 at 01:03
This is more of a math question than a programming question. In general, a piece-wise linear convex hull of a convex function (as shown in your question) is not straightforward. For example, see [math.stackexchange.com](https://math.stackexchange.com/questions/3762278/how-to-find-convex-hull-of-functions0) — mhopeng, Jul 27 '22 at 02:24
One possible approach is: sample the _value_ and the _derivative_ of the orignal function at certain points, and use linear pieces that pass through those points with those slopes — Luis Mendo, Jul 27 '22 at 10:20
@mhopeng I do not think so because bioconjugation gives me a continuous non-linear function, and I need a piecewise linear function. — Juan, Jul 27 '22 at 15:00
@LuisMendo ok, but how about error tolerance, or which is the better way to distribute the sample points? — Juan, Jul 27 '22 at 16:57
If you want to guarantee a maximum error or optimize the choice of points, that's a different issue, which seems harder (and wasn't contained in your original question) — Luis Mendo, Jul 27 '22 at 19:25
You can use 1st order Taylor approximation (essentially what @Luis Mendo suggested). Bounding the error term with your desired maximum error will give you the sampling interval, which is not necessarily uniform. Should be easy to implement. This is not a programming question in fact, suggest to ask on math stack exchange. — Yuri Feldman, Jul 29 '22 at 15:55

score 3 · Answer 1 · answered Jul 29 '22 at 14:59

This approach is based on Luis Mendo's comment. The idea is the following:

Select a number of points from the original curve, your final piecewise linear curve will pass through these points
For each point calculate the equation of the tangent to the original curve. Because your graph is convex, the tangents of consecutive points in your sample will intersect below the curve
Calculate, for each set of consecutive tangents, the x-coordinate of the point of intersection. Use the equation of the tangent to calculate the corresponding y-coordinate

Now, after reordering the points, this gives you a piecewise linear approximation with the constraints you want.

h = 0.001;
x = 0.000000001:h:1;
y = abs(log2(x));

% Derivative of function on all the points
der = diff(y)/h;

NPts = 10; % Number of sample points

% Draw the index of the points by which the output will pass at random
% Still make sure you got first and last point
idx = randperm(length(x)-3,NPts-2);
idx = [1 idx+1 length(x)-1]; 
idx = sort(idx);
x_pckd = x(idx);
y_pckd = y(idx);
der_pckd = der(idx);

% Use obscure math to calculate the points of intersection
xder = der_pckd.*x_pckd;
x_2add = -(diff(y_pckd)-(diff(xder)))./diff(der_pckd);
y_2add = der_pckd(1:(end-1)).*(x_2add-(x_pckd(1:(end-1))))+y_pckd(1:(end-1));

% Calculate the error as the sum of the errors made at the middle points
Err_add = sum(abs(y_2add-interp1(x,y,x_2add))); 

% Get final x and y coordinates of interpolant
x_final = [reshape([x_pckd(1:end-1);x_2add],1,[]) x_pckd(end)];
y_final = [reshape([y_pckd(1:end-1);y_2add],1,[]) y_pckd(end)];

figure;
plot(x,y,'-k');
hold on
plot(x_final,y_final,'-or')

You can see in my code that the points are drawn at random. If you want to do some sort of optimization (e.g. what is the set of points that minimizes the error), you can just run this a high amount of time and keep track of the best contender. For example, 10000 random draws see the rise of this guy:

Nice! If you don’t know if the curve is concave or convex, you can look for zero crossings of the 2nd derivative. Those points are part of the piecewise linear curve. In between, you’d use your method or the one in the OP, depending on the sign of the 2nd derivative. — Cris Luengo, Jul 29 '22 at 16:06
Oh, and if you create a piecewise linear function that is above the input function, and one that is below, with nodes at the same x positions, then you can average the two to get a good, non-biased approximation to the function! — Cris Luengo, Jul 29 '22 at 16:10
Thanks! It is also surprisingly time efficient, in addition to needing few sample points to be precise. I didn't expect it to go as well — BillBokeey, Jul 29 '22 at 16:12

Underestimation of f(x) by using a piecewise linear function

1 Answers1