I think the whole procedure only makes sense if x1 is a subset of x2, but here we go:
scale = 1.5;
x1 = [0,4,6,10,15,20]*scale;
y1 = [18,17.5,13,12,8,10];
x2 = [0,10.5,28]*scale;
y2= [18.2,10.6,10.3];
%% y2 is the target function and y1 is the function to be shifted
f = y1;
g = interp1(x2,y2,x1);
ff = [0,ones(size(x1))];
A = [f.',-eye(numel(x1));-f.',-eye(numel(x1))];
b = [g.',-g.'];
sol = linprog(ff,A,b);
a = sol(1)
% Compare with least-squares solution
a2 = sum(f.*g)/sum(f.^2)
And the next step is
min: max_i ( abs( a*f_i - g_i ) )
? :-)
