Help vectorizing code to improve speed

5 次查看（过去 30 天）
Michael Jones
Michael Jones 2023-8-30
评论： Michael Jones 2023-8-30
I have a function (below) that will be called many times. The profiler reveals that the slowest line in the function (by a lot) is:
G = [x(mask,i),y(mask,i),z(mask,i),ones(nnz(mask),1)];
I was surprised since this is simply assembling the G matrix. The matrix operations (rcond, and "\") are much faster. I tried bulding the G matrix outside of the for loop, but it is three dimensional, and I had to do a squeeze() on the matrix slices inside the loop, which was even slower.
Can anyone help vectorize this code more than it is already or help build the G matrix more efficiently?
Any other efficiency tips are appreciated.
Thank you!
function dop = get_dop(az, el, sat_mask, minel)
%Calulate Pdop using a vector a az/el values from visible satellites to a single ground point.
%Input:
% az: vector of azimuths for each satellite (row) at each time (column)
% el: vecotr of elevations for each satellite (row) at each time (column)
% sat_mask: vector of rows of az/el to include
% min_el: minimum elevation angle to include satellite row in calculation
%Output:
% pdop: vector of pdop values corresponding to az/el values for the given ground point
%
sat_mask = boolean(sat_mask);
az = squeeze(az(:,sat_mask,:));
el = squeeze(el(:,sat_mask,:));
el_mask = el >= minel;
az = deg2rad(az);
el = deg2rad(el);
[x,y,z] = sph2cart(az,el, ones(size(az)));
dop = nan(1,size(el,2));
for i = 1:size(el,2)
%Elevation mask to get when the sats are actually in view
mask = el_mask(:,i);
if nnz(mask) >= 3 %Ensure that at least 3 satellites are available before proceeding
G = [x(mask,i),y(mask,i),z(mask,i),ones(nnz(mask),1)]; %matrix of all the x,y,z coordinates for in view satellites
if rcond(G'*G) > 1e-8 %Rcond check to ensure matrix is invertible. If not, skip and dop remains NaN.
P = (G'*G)\eye(4);
dop(i) = sqrt(P(1,1)+P(2,2)+P(3,3));
end
end
end
end
  5 个评论
显示 3更早的评论
KSSV
KSSV 2023-8-30
You should initilaze the variables which you are filling in the loop. What is the function wrapTo360?
Michael Jones
Michael Jones 2023-8-30
WrapTo360 is from the mapping toolbox (and maybe some others) https://www.mathworks.com/help/map/ref/wrapto360.html
It wraps angles back into the interval [0, 360].

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bruno Luong
Bruno Luong 2023-8-30
编辑：Bruno Luong 2023-8-30
Ran in:
Here is a vectorized version of the for-loop.
It is not faster according to this kind of test size and mask density, I'm little bit surprise by that.
m = 50;
n = 2000;
x = rand(m,n);
y = rand(m,n);
z = rand(m,n);
el_mask = rand(m,n) > 0.2;
tic
dop = nan(1,n);
for i = 1:size(el_mask,2)
%Elevation mask to get when the sats are actually in view
mask = el_mask(:,i);
if nnz(mask) >= 3 %Ensure that at least 3 satellites are available before proceeding
G = [x(mask,i),y(mask,i),z(mask,i),ones(nnz(mask),1)]; %matrix of all the x,y,z coordinates for in view satellites
if rcond(G'*G) > 1e-8 %Rcond check to ensure matrix is invertible. If not, skip and dop remains NaN.
P = (G'*G)\eye(4);
dop(i) = sqrt(P(1,1)+P(2,2)+P(3,3));
end
end
end
toc
Elapsed time is 0.024660 seconds.
tic
[m,n] = size(el_mask);
sz3 = [m,1,n];
X = reshape(x,sz3);
Y = reshape(y,sz3);
Z = reshape(z,sz3);
I = ones(sz3);
A = [X,Y,Z,I];
M = reshape(el_mask,sz3);
[P,rc] = pageinv(pagemtimes(A.*M,'ctranspose',A,'none'));
dop2 = sqrt(P(1,1,:)+P(2,2,:)+P(3,3,:));
dop2(sum(M,1)<3 | rc<=1e-8) = NaN;
dop2 = reshape(dop2,[1,n]);
toc
Elapsed time is 0.015918 seconds.
relerr = norm(dop-dop2)/norm(dop)
relerr = 1.7188e-16
  6 个评论
显示 4更早的评论
Michael Jones
Michael Jones 2023-8-30
The Mac is an M1 Max and the PC is Intel i9-12900k. Both are running Matlab r2023a.
Michael Jones
Michael Jones 2023-8-30
Its unfortunate I can't generate a Mex file from the the vectorized version as it looks like pagemldivide and pageinv aren't supported for Mex generation. However, the vecotrized version is faster than the Mex compiled version of the original code, so it's still a win.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by