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

[the cyclist]

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]

Sure: here's a script that calls it with some dummy data. The sizes of the inputs are similar to my use case.

%Call get_dop
clear; clc;
for i =1:2000
    az(:,i) = wrapTo360(linspace(0,360,50)+i*5);
    el(:,i) = wrapTo360(linspace(0,90,50)+i*5);
end
sat_mask  = ones(size(az,1),1);
minel = 10;

dop = get_pdop(az,el,sat_mask,minel);

function dop = get_pdop(az, el, sat_mask, minel)
    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)
        mask = el_mask(:,i);
        if nnz(mask) >= 3 
            G = [x(mask,i),y(mask,i),z(mask,i),ones(nnz(mask),1)]; 
            if rcond(G'*G) > 1e-8
                P = (G'*G)\eye(4);
                dop(i) = sqrt(P(1,1)+P(2,2)+P(3,3));
            end
        end
    end
end

[Michael Jones]

I would imagine slicing x, y, and z is time consuming here, and I'm hoping there's a way to avoid that or do it better.

[KSSV]

You should initilaze the variables which you are filling in the loop. What is the function wrapTo360?

[Michael Jones]

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].

Answer 1

[Bruno Luong]

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

[Michael Jones]

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]

An alternative (could be slighly faster, and more TMW "approval")

[P, rc] = pagemldivide(pagemtimes(A.*M,'ctranspose',A,'none'), eye(4));

[Michael Jones]

Thanks, I'll try that as well.

Interestingly, the vectorized version is much faster (4sec) than the for loop version (19sec) on my Mac. On my PC it is only somewhat faster - 9sec for vectorized version vs 13.5sec for for loop.

[Bruno Luong]

What Mac is it? Does it have M1 or M2 processor?

[Michael Jones]

The Mac is an M1 Max and the PC is Intel i9-12900k. Both are running Matlab r2023a.

[Michael Jones]

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.