How to pass pre-generated values to an anonymous function to speed up a fit?

Question

In Matlab I am using the combination of fit and fittype function (from the curve fitting toolbox) to fit a signal using an anonymous function.

My issue is that, this anonymous function use the integral function of matlab, and this is extremely time consuming. I am therefore trying to get rid of this function inside the anonymous function.

An amazing improvement would be to pre-calculate these integral and pass it to the fit function, but I do not know how to do this.

I tried to define the vector that contain these integration values and pass it to the fit by defining it as 'dependent' or 'problem' variable in fittype but that did not do the trick. I also tried concatenating the coordinate and the result of this integral. When doing a two column vector it tried to fit a surface (which it is not the case). When concatenating in a one column vector it complain that the data set and the coordinate system should be the same size.

As an alternative I also tried using sum and trapz instead of integral but fittype told me that then the function did not fulfill the criteria of an anonymous function.

I am running short on ideas and on ways to ask search engines for new ideas and would welcome all new point of view on how to solve this.

---

For curious peoples here is the original anonymous function I am using:

function [ y ] = zeusEquation(varargin)
 psi0=varargin{1};% in radian 
 a1 = varargin{2};% parameter to be fitted
 a2 = varargin{3};%parameter to be fitted
 a3 = varargin{4};% average of background signal/ offset in gray level
 a4 = varargin{5};% slope on the signal background
 a5 = varargin{6};% decentering of the signal (in meter)
 a6 = varargin{7};% width of the signal in radian
 x = varargin{8};%Coordinate system in meter!!!!
 y = zeros(size(x));
 %fundamental bloc equation needed
 sinIntegral = @(t) 2.*sin(t)./(pi.*t);
 u = @(f, e, x) f.*abs(x-e);
 %Actual signal calculation
 for i = 1:max(size(x))
 if u(a6,a5,x(i))~=0
 y(i) = a3...
 +a4.*(x(i)-a5)...
 +(1-integral(sinIntegral,0,u(a6,a5,x(i))))...
 .*(a1.*sign(x(i)-a5).*sin(psi0)-...
 (1-a2).*integral(sinIntegral,0,u(a6,a5,x(i))).*(1-cos(psi0)));
 else
 y(i) = a3+a4.*(x(i)-a5);
 end
 end
end

Note: I cross asked this question on mathwork.

• ZF007: ..thumbs up for the manual move from SO ;-)
• Cris Luengo: An anonymous function is one without a name. Yours is called zeusEquation and hence not an anonymous function. An anonymous function in MATLAB is a lambda, and it can capture variables of the scope it is defined in. However, regular named functions defined inside another function can do so too, sort of combining the lambda and named function concepts. If I have time later I'll try to write an answer to describe the capture concept.

Answer 1

A named function is declared like this:

function out = function_name(in)

If it is written in a file of the same name, and the file is on the MATLAB path, then it is accessible from the base workspace and any function you call. If it is written as a local function (an additional function in a file defining a different function) then it is only accessible to the other functions in that file. If it is in a sub-directory called private, then it is only accessible to the functions defined in the parent directory.

An anonymous function is defined as like this:

function_handle = @(p) p

Here, function_handle is a variable that contains the anonymous function. In other languages this is called a lambda. This lambda can capture variables from its environment:

data = 1:10;
function_handle = @(p) p.*data

function_handle now contains a copy of data. Changing data or clearing it after creating the anonymous function won't affect the copy of data inside it.

The interesting thing for you here is that a named function can also capture data. This is true only for nested functions (which we didn't mention above). The difference between a nested function and a local function is one measly end statement, it's easy to miss it. This is a local function:

function out = main(in)
% code ...
function r = local(p)
% code ...

This is identical to:

function out = main(in)
% code ...
end
function r = local(p)
% code ...
end

The end statements ending a function are optional.

This is a nested function:

function out = main(in)
% code ...
 function r = nested(p)
 % code ...
 end
end

Here the end statements are mandatory. nested is inside main. nested also has access to variables declared in main. For example, nested can use in.

main can also create a function handle of nested and return it. This function handle then will point at the nested function, and capture its environment. Much like an anonymous function:

function out = main(in)
 data = integral(in,0,1);
 out = @r;
 function r = nested(p)
 r = p .* data;
 end
end

Now, calling main returns a handle to the function nested, which has captured data computed in main. You can use this to create a function handle that you can pass to fit. nested would be your zeusEquation, and main would simply package it with the precomputed integral.

Code review

There are a few things to note about your code:

• You compute sin(psi0) and 1-cos(psi0) repeatedly inside your loop. psi0 is a constant within this loop, and consequently you should compute these values once outside of the loop.

• You compute integral(sinIntegral,0,u(a6,a5,x(i))) twice for each x(i). This being an expensive function, you should compute it only once and re-use the result.

• for i = 1:max(size(x)) is long-hand for for i = 1:length(x). This only works correctly in your loop if x is a vector. In this case, you can use numel(x) instead of length(x). numel is faster.

These three points above lead to:

sinpsi0 = sin(psi0);
cospsi0_1 = 1-cos(psi0);
for i = 1:numel(x)
 if u(a6,a5,x(i))~=0
 int = integral(sinIntegral,0,u(a6,a5,x(i)));
 y(i) = a3 + a4.*(x(i)-a5) ...
 + (1-int) .* (a1.*sign(x(i)-a5).*sinpsi0 ...
 - (1-a2).*int.*cospsi0);
 else
 y(i) = a3 + a4.*(x(i)-a5);
 end
end

• integral must be called with scalar limits, but other components can be computed vectorised. It would be possible to compute only the integral within the loop, and the rest without a loop.

  This leads to:

  int = zeros(size(x));
  for i = find(u(a6,a5,x)~=0)
   int(i) = integral(sinIntegral,0,u(a6,a5,x(i)));
  end
  y = a3 + a4.*(x-a5) ...
   + (1-int) .* (a1.*sign(x-a5).*sin(psi0) ...
   - (1-a2).*int.*(1-cos(psi0)));
  

  This should be faster, but I haven't tested it. (I assume a1 and the like are all scalars?)

• You define your functionzeusEquation using varargin, then proceed to unpack the varargin cell array with the assumption that all your input arguments are there. This is correct in your situation of course, but you might as well do:

  function y = zeusEquation(psi0,a1,a2,a3,a4,a5,a6,x)
  

  In this way, if zeusEquation is called with fewer (or more!) arguments, you'll get a sensical error message, rather than an "index out of bounds" error.

• ALC: Dear Cris, Thank you for your answer, I tested your version after the fourth bullet: it definitively is faster: over 200 tries my former version takes ~6 hours, your solution takes~ 3 hours. This is a great improvement but is there a way to do even better? Regarding the first part of your answer: I feel you are trying to tell me something very important but I am afraid I have a hard time understanding it. Are you suggesting that instead of saving zeusEquation in a separate file I should write it as a nested function in the same script where I call fit?
• Cris Luengo: @ALC: Yes, writing it as a nested function of the function that calls fit is one way of doing it. You can then calculate the integral for all possible limits in the main function, and use that result in the zeusEquation. --- The other way of doing it is writing a function zeusEquationGenerator that computes the integral and returns a function handle to zeusEquation as its output argument. When you call fit you use the new function to create the function handle you pass into it. --- Either way you'd have what you are asking for: pre-computed integral values.
• ALC: Dear Chris, So far I have not been able to implement the nested function solution properly. As soon as I rely on the precomputed integral value the calculus can only be done for the same independant variable. but fittype seems to not like this. I have tried to trick the thing with an if loop and putting the original code but it seems more clever than this. I am still trying to understand what I am doing wrong. As for your second suggestion: I though doing so would recalculate every time the integral, What did I miss?
• Cris Luengo: @ALC: Sorry, I saw your comment but I forgot to answer. If your integral must be recomputed during the optimization procedure, then of course you can't pre-compute the integral. The above was written assuming the integral result is a constant. You must re-generate the function handle for every new fitting problem, so that the pre-computed integral can be updated with the new problem's constants. The part of the answer under "Code Review" does recompute the integral. Those snippets are written to illustrate other things.