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');

Answer 1

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

1. Select a number of points from the original curve, your final piecewise linear curve will pass through these points
2. 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
3. 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: