Re: fixed points piecewise-linear fitting

[Ben Abbott]

On Mar 20, 2012, at 8:11 AM, Martin Helm wrote:

> 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-")

After more testing I concluded the problem isn't a linear one. Meaning that 
your linear solution isn't working as I'd hoped.

For my example the RMSE was about 70 for your solution and about 10 for my 
fsolve() approach.

I did find that your linear approach provided a reasonable initial guess, so I 
am using that. I'll clean things up and attach the script later.

Ben