Re: fixed points piecewise-linear fitting

[Martin Helm]

Am 20.03.2012 02:18, schrieb Ben Abbott:
> On Mar 19, 2012, at 8:42 PM, Martin Helm wrote:
>
>> I want to add my poor man's solution for least squares fit of the
>> piecewise polynomials
>> returns the node values as first argument (makes only sense for the
>> piecewise linear case) and the function values at all nodes as second
>> argument.
>> Simplistic but fast.
>>
>> function [y, yy] = ppfit2(x, xf, yf, method)
>>  a = zeros (length (xf), length (x));
>>  for ii = 1:length (x)
>>    p = zeros (length (x), 1);
>>    p(ii) = 1;
>>    a(:, ii) = interp1 (x(:), p, xf(:), method);
>>  endfor
>>  y = a\yf(:);
>>  if (nargout==2)
>>    yy = interp1 (x(:), y, xf(:), method);
>>  endif
>> endfunction
>>
>> %!demo
>> %! clf
>> %! x = [0 1 2 3];
>> %! xf = linspace(0, 3, 1000);
>> %! yf = xf.^3 + rand(size(xf))*2.5;
>> %! plot (xf, yf, "y-")
>> %! hold on
>> %! [y, yy] = ppfit2(x,xf,yf, "linear");
>> %! plot (xf, yy, "r-")
>> %! [y, yy] = ppfit2(x,xf,yf, "cubic");
>> %! plot (xf, yy, "g-")
>> %! [y, yy] = ppfit2(x,xf,yf, "spline");
>> %! plot (xf, yy, "b-")
> Nice!
>
> I ran a comparison between your linear least squares solution and my fsolve 
> approach. With the default optimization values they are about the same speed.
>
> However, I prefer linear least squares, so I'm modifying my local function.
>
> Ben
>
If you want it in a way that the x nodes and the xf nodes do not need to
have the same end points add an "extrap" to the the two interp1 calls
(as you did in your own example).

I think the "Ansatz" to use finite elements for solving the 1d Laplace
equation (piecewise linear) or biharmonic equation (piecewise cubic) is
worth to check it.

Your own "ppolyfit" code has of course the benefit that you can in
principle make it using any possible distance function which is very
nice and I like it.
I tried myself something similar with fminunc, but it is very slow (p is
the norm 1,2, Inf) and it is no production code (no possibility to pass
options to fminunc or proper error handling, documentation ...).

function [y, yy] = ppfit(x, xf, yf, method, p)
  dist = @(y) norm (interp1 (x, y, xf, method, "extrap") - yf, p);
  y = fminunc (dist, interp1 (xf, yf, x));
  if (nargout==2)
    yy = interp1 (x, y, xf, method, "extrap");
  endif
endfunction

%!demo
%! clf
%! x = [0 1 2 3];
%! xf = linspace(0, 3, 100);
%! yf = xf.^3 + rand(size(xf))*2.5;
%! plot (xf, yf, "g-")
%! hold on
%! [y, yy] = ppfit(x,xf,yf, "linear", 1);
%! plot (xf, yy, "r-")
%! [y, yy] = ppfit(x,xf,yf, "cubic", Inf);
%! plot (xf, yy, "b-")
%! [y, yy] = ppfit(x,xf,yf, "spline", 2);
%! plot (xf, yy, "b-")