Help vectorizing code to improve speed
显示 更早的评论
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 个评论
the cyclist
2023-8-30
That line is not just assembling G; it is also all the indexing and matrix creation that that line contains.
Can you give an example of the parameters used to call that function, so that we don't have to guess?
Michael Jones
2023-8-30
Michael Jones
2023-8-30
KSSV
2023-8-30
You should initilaze the variables which you are filling in the loop. What is the function wrapTo360?
Michael Jones
2023-8-30
采纳的回答
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
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
relerr = norm(dop-dop2)/norm(dop)
6 个评论
Michael Jones
2023-8-30
Thank you! This is what I was looking for. I tried it with my use case and it is faster for the sizes I'm working with.
Here's a screenshot from Matlab Profiler from the vectorized code:
Here's a screenshot from Matlab Profilers from the original with a for loop:
The Profiler ran a script which called the get_dop function 1260 times with [mxn] sizes for az/el of [136x1441] and the mask density was pretty high.
Bruno Luong
2023-8-30
编辑：Bruno Luong
2023-8-30
An alternative (could be slighly faster, and more TMW "approval")
[P, rc] = pagemldivide(pagemtimes(A.*M,'ctranspose',A,'none'), eye(4));
Michael Jones
2023-8-30
Bruno Luong
2023-8-30
What Mac is it? Does it have M1 or M2 processor?
Michael Jones
2023-8-30
Michael Jones
2023-8-30
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
