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
采纳的回答
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.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Matrix Indexing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
