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.

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.