How do I vectorize this code?

Question

I have written a recursive function, however, it takes a lot of time. Hence I vectorized it, but it does not yield the same result as the recursive function. This is my non-vectorized code:

function visited = procedure_explore( u, adj_mat, visited )
visited(u) = 1;
neighbours = find(adj_mat(u,:));
for ii = 1:length(neighbours)
    if (visited(neighbours(ii)) == 0)
        visited = procedure_explore( neighbours(ii), adj_mat, visited );
    end
end
end

This is my vectorized code:

function visited = procedure_explore_vec( u, adj_mat, visited )
visited(u) = 1;
neighbours = find(adj_mat(u,:));
len_neighbours=length(neighbours);
visited_neighbours_zero=visited(neighbours(1:len_neighbours)) == 0;
if(~isempty(visited_neighbours_zero))
    visited = procedure_explore_vec( neighbours(visited_neighbours_zero), adj_mat, visited );
end
end

This is the test code

function main
    adj_mat=[0 0 0 0;
             1 0 1 1;
             1 0 0 0;
             1 0 0 1];
    u=2;
    visited=zeros(size(adj_mat,1));
    tic
    visited = procedure_explore( u, adj_mat, visited )
    toc
    visited=zeros(size(adj_mat,1));
    tic
    visited = procedure_explore_vec( u, adj_mat, visited )
    toc
end

This is the algorithm I'm trying to implement:

If vectorization is impossible, a mex solution would also be good.

Update benchmark: This benchmark is based on MATLAB 2017a. It shows that the original code is faster than other methods

Speed up between original and logical methods is 0.39672
Speed up between original and nearest methods is 0.0042583

Full code

function main_recersive
    adj_mat=[0 0 0 0;
             1 0 1 1;
             1 0 0 0;
             1 0 0 1];
    u=2;
    visited=zeros(size(adj_mat,1));
    f_original=@()(procedure_explore( u, adj_mat, visited ));
    t_original=timeit(f_original);

    f_logical=@()(procedure_explore_logical( u, adj_mat ));
    t_logical=timeit(f_logical);

    f_nearest=@()(procedure_explore_nearest( u, adj_mat,visited ));
    t_nearest=timeit(f_nearest);

    disp(['Speed up between original and logical methods is ',num2str(t_original/t_logical)])
    disp(['Speed up between original and nearest methods is ',num2str(t_original/t_nearest)])    

end

function visited = procedure_explore( u, adj_mat, visited )
    visited(u) = 1;
    neighbours = find(adj_mat(u,:));
    for ii = 1:length(neighbours)
        if (visited(neighbours(ii)) == 0)
            visited = procedure_explore( neighbours(ii), adj_mat, visited );
        end
    end
end

function visited = procedure_explore_nearest( u, adj_mat, visited )
    % add u since your function also includes it.
    nodeIDs = [nearest(digraph(adj_mat),u,inf) ; u];
    % transform to output format of your function
    visited = zeros(size(adj_mat,1));
    visited(nodeIDs) = 1;

end 

function visited = procedure_explore_logical( u, adj_mat )
   visited = false(1, size(adj_mat, 1));
   visited(u) = true;
   new_visited = visited;
   while any(new_visited)
      visited = any([visited; new_visited], 1);
      new_visited = any(adj_mat(new_visited, :), 1);
      new_visited = and(new_visited, ~visited);
   end
end

Answer 1

Here's a fun little function that does a non-recursive breadth-first search on the graph.

function visited = procedure_explore_logical( u, adj_mat )
   visited = false(1, size(adj_mat, 1));
   visited(u) = true;
   new_visited = visited;

   while any(new_visited)
      visited = any([visited; new_visited], 1);
      new_visited = any(adj_mat(new_visited, :), 1);
      new_visited = and(new_visited, ~visited);
   end
end

In Octave, this runs about 50 times faster than your recursive version on a 100x100 adjacency matrix. You'll have to benchmark it on MATLAB to see what you get.

Answer 2

You can think of your adjacency matrix as a list of paths of length exactly one. You can generate paths of other lengths n by taking it to the n-th power up to the rank of your matrix. (adj_mat^0 is the identity matrix)

In a graph with n nodes, the longest path cannot be longer than n-1, therefore you can sum all the powers together for a reachability analysis:

adj_mat + adj_mat^2 + adj_mat^3
ans =
   0   0   0   0
   4   0   1   3
   1   0   0   0
   3   0   0   3

This is the number of (different) ways that you can use to go from one node to another. For simple reachablitity, check whether this value is greater than zero:

visited(v) = ans(v, :) > 0;

Depending on your definition, you may have to change columns and rows in the result (i.e. take ans(:, v)).

For performance, you can use the lower powers to make higher ones. For example something like A + A^2 + A^3 + A^4 + A^5 would be efficiently calculated:

A2 = A^2;
A3 = A2*A
A4 = A2^2;
A5 = A4*A;
allWalks= A + A2 + A3 + A4 + A5;

Note: If you want to include the original node as being reachable from itself, you should include the identity matrix in the sum.

This minimizes the number of matrix multiplications, also MATLAB will likely execute a matrix square faster than a regular multiplication.

In my experience, matrix multiplication is relatively fast in MATLAB and this will yield the result (reachability) vector for all the nodes in the graph at once. If you are only interested in a small subset of a large graph, this is probably not the best solution.

See also this answer: https://stackoverflow.com/a/7276595/1974021