Vectorize 2d chi-square grid search

Question

Best-fit linear parameters A and B (y=Ax+b) correspond to the minimum of the chi-square function over these parameters. I want to do a brute force grid search for the global minimum (guaranteed because 2-parameter linear chi-square is a paraboloid) and have achieved it with 3 nested loops (below) but want to avoid loops (i.e., vectorize using array properties).

Chi-square (weighted least squares) is defined as (pseudocode):

Chi-square(k,j) = sum (y[i]-(A[k]*x[i]+B[j]))/yerr[i])^2.

Below is Matlab code that fills a 100 x 100 grid with chi-square values over the 10,000 combinations of A and B parameter values (100 values each). There are three data arrays: x, y and yerr.

Thanks for any help towards a loopless version of a 2-parameter linear chi-square grid!

Keith

 % create parameter grid
  a = linspace(85,110,100);
  b = linspace(10,35,100);
  [A,B] = meshgrid(a,b);

  % calculate chi-square over parameter grid
  chi2(100,100) = zeros;

  for k = 1:100;
      for j = 1:100;
          for i = 1:length(y)
          chi2a = ((y(i)-a(k)*x(i)-b(j))/yerr(i)).^2;
          chi2(k,j) = chi2(k,j)+chi2a;
          end
      end
  end

Divakar Divakar · Accepted Answer · 2017-05-24T15:29:29

We could bsxfun it -

x3d = reshape(x,1,1,numel(x));
y3d = reshape(y,1,1,numel(y));
yerr3d = reshape(yerr,1,1,numel(yerr));
p0 = bsxfun(@minus, bsxfun(@minus,y3d,bsxfun(@times,a(:),x3d)), b);
p1 = bsxfun(@rdivide, p0, yerr3d);
out = sum(p1.^2,3);

With MATLAB's implicit expansion, computing p0 and p1 would simplify to -

p0 = ((y3d - a(:).*x3d) - b);
p1 = p0 ./yerr3d;

Timings -

% Setup
N = 2000;
x = rand(N,1);
y = rand(N,1);
yerr = rand(N,1);

a = linspace(85,110,100);
b = linspace(10,35,100);

we get -

----------- With loopy method -------------------------
Elapsed time is 1.056787 seconds.
----------- With BSXFUN method -------------------------
Elapsed time is 0.109601 seconds.

Vectorize 2d chi-square grid search

1 Answers