Distinguish between linear and nonlinear model parameters

[stestoll]

For multi-parameter models, such as multi-Gauss distribution models, it would be great to distinguish linear and nonlinear parameters. Then we could use separable nonlinear least-squares for the fitting, just as we do for regularization. This could result in performance increase, and possibly also in less failures.

For an n-Gauss model, this would reduce the nonlinear parameters from 3n-1 to 2n.

[luisfabib]

This is a very good idea and we should do some benchmarks to see the potential performance and quality changes before any major implementation.

This also raises the question of whether we should have a dedication separable non-linear least-squares solver. A function with some simplified interface which allows simple calls from other functions and which avoids code duplication.

[luisfabib]

I have done some preliminary tests to assess whether there is potential in using a nested approach while differentiating linear from non-linear parameters.

With the implementation above, I could observe the following:

• There is no potential speedup from using a nested approach. Furthermore, in most cases the nested approach leads to slower execution times.

• The nested approach seems to be more stable between noise realizations. However, the standard approach seems to find the true solution more often

• The key decision to find the optimal solution is still to use the 'MultiStart' option to find the global minimum.

Further study is still needed to take an informed decision on the matter.

clear,clc,clf

t = linspace(0,5,600);
r = linspace(2,6,400);

% This one always fails
P = dd_gauss2(r,[4 0.5 0.4 3.5 0.2]);
% This fails less often
% P = dd_gauss2(r,[4 0.5 0.4 4.5 0.2]);

K = dipolarkernel(t,r);
V = K*P + whitegaussnoise(t,0.01);

% Standard non-linear optimization
%--------------------------------------
info = dd_gauss2;
par0 = info.Start;
lower = info.Lower;
upper = info.Upper;
tic
[parfit,Pfit] = fitparamodel(V,@dd_gauss2,r,K,par0,lower,upper,'multistart',1,'TolFun',1e-8);
toc
Vfit = K*Pfit;

% Nested optimization
%--------------------------------------
model = @(t,par)myfcn(t,par,r,K,V);
par0 = par0([1 2 4 5]);
lower = lower([1 2 4 5]);
upper = upper([1 2 4 5]);
tic
parfit2 = fitparamodel(V,model,t,par0,lower,upper,'multistart',1,'TolFun',1e-8);
toc
[Vfit2,Pfit2] = model(t,parfit2);

%Plots
subplot(211)
plot(t,V,'k.',t,Vfit,'r',t,Vfit2,'b--')
axis tight,grid on

subplot(212)
plot(r,P,'k',r,Pfit,'r',r,Pfit2,'b--')
axis tight,grid on

function [Vfit,Pfit] = myfcn(t,par,r,K,V)

%Basis functions
p1 = dd_gauss(r,par(1:2));
p2 = dd_gauss(r,par(3:4));

% Optimize the amplitudes of the basis functions
amps = lsqnonneg([K*p1 K*p2],V);

% Construct resulting distribution
Pfit = amps(1)*p1 + amps(2)*p2;
Vfit = K*Pfit;

end

[luisfabib]

Closed by the new overhaul of fitmultimodel at 579494d