How to vectorize matrix definition without using nested functions

I have a matrix H with dimensions N x M. I also have a few (1 x N) vectors denoted here as n1 and n2, and a few (1 x M) vectors denoted here as m1 and m2.

I want to define H such that H_(n,m) is a function of n1(n), n2(n), m1(m), and m2(m) without a for loop. Below is my attempt with arrayfun:

function [H] = solveH(N, M)
  % Initialize n1, n2, m1, and m2
  n1 = <some code here>;
  n2 = <some code here>;
  m1 = <some code here>;
  m2 = <some code here>;
  
  % NMat(i) and MMat(i) contain all subscript pairs for H
  [MMat, NMat] = meshgrid(m, n);
  
  H = arrayfun(@calcH, NMat, MMat);
  
  % inputs n and m are scalars
  function outputH = calcH(n, m)
    outputH = <some function of n1(n), n2(n), m1(m), and m2(m)>;
  end
end

This works fine, except that it is actually slower (!) than using a for loop because n1, n2, m1, and m2 have a scope in multiple functions and are slow like global variables.

I also considered turning n1, n2, m1, and m2 into matrices and doing element-by-element operations, but in reality I have something like five (1 x N) vectors and five (1 x M) vectors, and using repmat that many times wastes memory and makes it even slower.

What approach should I use so that my program executes the fastest?

EDIT: Here is most of the function for those who want to see exactly what operations I'm using.

% Solves the bold matrices H and E analytically
% ka      wavenumber in air    (1 x M vector)
% ks      wavenumber in bar    (1 x M vector)
% a      air gap  of HCG      (scalar)
% s      bar width of HCG    (scalar)
% nBar      refractive index    (scalar)
% period    period of HCG      (scalar)
% lambda    wavelength of light    (scalar)
% M      orders - region II    (scalar)
% N      orders - regions I, III    (scalar)
function [H, E, beta] = solveHE(ka, ks, a, s, nBar, period, lambda, N, M, polarization)
  
  n = 1:1:N;  
  m = 1:1:M;
  
  % Precalculate constants
  gamma = conj( 2*pi * sqrt(1/lambda^2 - (n-1).^2/period^2) ); % (1 x N)
  beta = conj( sqrt((2*pi*nBar/lambda)^2 - ks.^2) );  % (1 x M vector)
  p = (2*pi/period) * (n-1);    % (1 x N vector)
  aa = a/2;
  ss = s/2;
  dif = period - aa;
  nInv = (1 / nBar^2);
  
  % NMat(i) and MMat(i) contain all subscript pairs for H and E
  [MMat, NMat] = meshgrid(m, n);
  
  % Efficient vectorization?
  [H, E] = arrayfun(@calcHE, NMat, MMat);
  
  % Calculate H_(n,m) and E_(n,m)
  function [resH, resE] = calcHE(n, m)
    d = ~(n-1);  % delta = 1 for n = 0 (index 1), 0 otherwise
    d = 1 / period * (2-d);    % precalculate
    
    hAir = d*cos(ks(m)*ss) * 2 * eval(ka(m),aa,p(n),aa);
    hBar = d*cos(ka(m)*aa) * (eval(ks(m),ss,p(n),dif) -                      eval(ks(m),-ss,p(n),aa));
    
    resH = hAir + hBar;
    
    if (strcmp(polarization, 'TM'))
      resE = beta(m)/gamma(n) * (hAir + nInv*hBar);
    else
      resE = gamma(n)/beta(m) * resH;
    end
  end
  
  % All parameters are scalars
  % 
  % sin(KU - PA)     sin(KU + PA)
  % ------------  +  ------------ 
  %    2(K-P)           2(K+P)
  function z = eval(K, U, P, A)
    z = sin(K*U-P*A)/(2*(K-P)) + sin(K*U+P*A)/(2*(K+P));
  end
end

Answer 1

[Matt J]

using repmat that many times wastes memory and makes it even slower.

You might also be able to reduce memory requirements by replacing REPMAT calls with BSXFUN, but we would have to see the details of the operations being done to know how applicable that would be.

In any case, you shouldn't need to be calling repmat once you have NMat and MMat. Simply doing n1(NMat) will turn n1 into an array and similarly for the others. That will save overhead, but not memory.

[Matt J]

You could also try

 H=bsxfun(@calcH,(1:N).', 1:M);

but that might not be the most optimal use of bsxfun.

[Vince]

Good point on indexing to avoid repmat. Now it runs so much faster--thanks!