K-th order neighbors in graph

I would like to do this, but in Matlab with a graph (it is not a directed graph): https://stackoverflow.com/questions/18393842/k-th-order-neighbors-in-graph-python-networkx

I have a graph in which I want to efficiently find a list of all K-th order neighbors of a node (j). K-th order neighbors are defined as all nodes which can be reached from the node in question in exactly K hops.

For example,

s = [1 1 1 1 2 2 2 2 2 8 8 12 6];
t = [3 5 4 2 6 10 7 9 8 11 12 13 7];
G = graph(s,t);
plot(G)

If j=1, I want to get :

for k=1

[2 3 4 5]

for k=2

[6 7 8 9 10]

for k=3

[11 12]

k=4

[13]

k=5

empthy

Answer 1

[Matt J]

function  result=distk(G,i,k)
%  G -graph
%  i - start node
%  k - number of steps

   A=adjacency(G);

    v=A(i,:);
    prior=false(size(v)); prior(i)=true;
    for n=1:k-1
       prior=prior|v;
       v=v*A;
    end
    
    result  = find(v&~prior);
    
end

[Bruno Luong]

Yes I wonder that same thing.

The k-neighbor of the node #j is simply

A = adjacency(G);
Ak = A^k;
i = find(Ak(:,j))

[Rub Ron]

It does not work.

s = [1 1 1 1 2 2 2 2 2 8 8 12];
t = [3 5 4 2 6 10 7 9 8 11 12 13];
G = graph(s,t);
plot(G)

k=3;

A = adjacency(G);
Ak = A^k;
ix = find(Ak(:,1))

output:

ix =

     2
     3
     4
     5
    11
    12

for k=3, ix should be [11, 12]

[Matt J]

Since the graph is non-directed, there is nothing invalid about the solution. A k=3-step path from 1 to 2 , for example, is 1-2-8-2

[Rub Ron]

Sorry, it is not what I was looking for. I eddited my post to make it clearer.

[Matt J]

If the rule is that k must also be the shortest path length to the resulting node set, then

k=3; %steps
i=1; %start node

v=A(i,:);
accum=logical(v);
for n=1:k-1
    
   v=v*A;

   if n<k-1
    accum=accum|v;
   end
    
end

result  = find(v&~accum)

result =

    11    12

[Rub Ron]

Thanks, but for k=2

result =

     1     6     7     8     9    10

and for k=1

result =

  1×0 empty double row vector

[Matt J]

for k=1:5

    i=1; %start node
    v=A(i,:);
    accum=false(size(v)); accum(i)=true;
    for n=1:k-1


       if n<k-1
        accum=accum|v;
       end

       v=v*A;

    end
    result  = find(v&~accum)

end

result =

     2     3     4     5


result =

     6     7     8     9    10


result =

    11    12


result =

    13


result =

  1×0 empty double row vector

[Rub Ron]

Thanks Matt, indeed it works for the given example, however, if I test the solution for other graph:

s = [1 1 1 1 2 2 2 2 2 8 8 12 14 14 1 14];
t = [3 5 4 2 6 10 7 9 8 11 12 13 7 6 8 15];
G = graph(s,t);

I get this:(for k=1:5)

result =
     2     3     4     5     8

result =
     2     6     7     8     9    10    11    12

result =
     6     7     9    10    11    12    13    14

result =
    13    14    15

result =
    15

But, I would expect this:

result =
     2     3     4     5     8

result =
     6     7     9    10    11    12

result =
     13    14

result =
     15

result =

  1×0 empty double row vector

[Matt J]

OK, how about this,

A=adjacency(G);

i=1;

for k=1:5

    i=1; %start node
    v=A(i,:);
    prior=false(size(v)); prior(i)=true;
    for n=1:k-1

       prior=prior|v;
       v=v*A;

    end
    result  = find(v&~prior)

end

[Rub Ron]

Now this works!

Bruno also provide this way.

d = distances(G,start);
% neighhbor for any specific k, e.g. k=2; neightbors i are
% i = find(d==k)

It works too. The only restriction is that if the graph egdes have weighted values, the output might not be as desired.

can you updated your answer with this comment to accept it?

[Matt J]

I've updated the answer. A drawback of Bruno's approach that I see is that its use of distances() will compute all distances, even those greater than the particular k of interest.

[Bruno Luong]

Those two comments on my answer is easily overcome.

If a fix small distance is requested, use the while-loop if BF method by a for-loop. The distance vector increases by maximum 1 each time.

Copy the weighted graph to an unweighted one and you are done.

UnweightedG = graph(spones(adjacency(WeightedG))); 

[Rub Ron]

Alternatively, one could remove the weight before using distances approach

G.Edges.Weight=[];

[Bruno Luong]

"A drawback of Bruno's approach that I see is that its use of distances() will compute all distances, even those greater than the particular k of interest."

I just did some benchmark, and Matlab graph distances seems to be pretty fast and beat all even for small k and large n. Let alone the case where k is large.

If that is matter at all.

[Bruno Luong]

Here is my benchmark code, feel free to play with it.

function BenchGNeighbor
start = 1; % starting node

% Graph corresponds to the 2D square grid of size 1000 x 1000
[X,Y] = ndgrid(1:1000);
n = numel(X);
I = sub2ind(size(X),X(1:end-1,:),Y(1:end-1,:));
J = I+1;
A = sparse(I,J,1,n,n);
I = sub2ind(size(X),X(:,1:end-1),Y(:,1:end-1));
J = I+size(X,1);
A = A + sparse(I,J,1,n,n);
A = A + A';

G = graph(A);

k = 100; %size(X,1);

% clean memory
clear X Y I J

k = min(k,size(A,1));

%% Matlab graph distances method
tic
d = distances(G, start);
rMatlab = find(d==k);
tMatlab = toc;

%% Bruno's algorithm
tic
n = size(A,1);
% initialize distance matrix
notdone = true(n,1);
i = start;
for q=1:k
    notdone(i) = false;
    [i,~] = find(A(:,i));
    if isempty(i)
        break
    end
    i  = unique(i);
    i = i(notdone(i));
end
rBruno = reshape(i,1,[]);
tBruno = toc;

%% Matt method
tic
i = start; % start node
v = A(i,:);
prior = false(size(v));
prior(i) = true;
for j=1:k-1
    prior = prior|v;
    v = v*A;
end
rMatt  = find(v & ~prior);
tMatt = toc;

% Check matching results
if ~isequal(rBruno, rMatlab) || ~isequal(rMatt, rMatlab)
    fprintf('Bug detected!!!\n');
end

fprintf('Timings for 2D-grid graph with number of nodes = %d, k=%d\n', n, k);
fprintf('\ttMatlab = %f [second]\n', tMatlab);
fprintf('\ttBruno  = %f [second]\n', tBruno);
fprintf('\ttMatt   = %f [second]\n', tMatt);

Answer 2

[Bruno Luong]

It sounds like you look at graph-distance and NOT what you described "K-th order neighbors are defined as all nodes which can be reached from the node in question in exactly K hops." The later problem is solved by my other answer.

If it is is the first case (graph distance) one can do by shortest path algorithms such as Bellman-Ford (BF) algorithm.

Graph example (yours):

s = [1 1 1 1 2 2 2 2 2 8 8 12 14 14 1 14];
t = [3 5 4 2 6 10 7 9 8 11 12 13 7 6 8 15];
G = graph(s,t);
plot(G)

BF algo (note: the bellow is not optimal implemented, since using full matrix with Inf).

start = 1; % starting node
B = adjacency(G);
B = full(B);
B(B==0) = Inf;
n = size(B,1);
d = inf(n,1);
d(start) = 0;
while any(isinf(d)) % careful graph must be connected, otherwise loop forever
    d = min(d,min(d+B,[],1).');
end

Or simply using MATLAB stock function

d = distances(G,start);
% neighhbor for any specific k, e.g. k=2; neightbors i are
% i = find(d==k)

Display result

n = size(G.Nodes,1);
nodes = setdiff(1:n,start)';
d = d(:);
dcell = accumarray(d(nodes), nodes, [], @(x) {x.'});
K = (1:length(dcell))';
for k=1:length(dcell)
    i = dcell{k};
    fprintf('k = %d: i = %s\n', k, mat2str(i));
end

Result

k = 1: i = [2 3 4 5 8]
k = 2: i = [6 7 9 10 11 12]
k = 3: i = [13 14]
k = 4: i = 15

[Bruno Luong]

In the above I code Bellman-Ford algorithm like this single-line while loop for illustration purpose:

% ...
% Initialization of d and B
while any(isinf(d)) % careful graph must be connected, otherwise loop forever
    d = min(d,min(d+B,[],1).');
end

A better implementation of Bellman-Ford algorithm that exploits the sparsity of ajadcent matrix to compute the shortest distance is

start = 1; % starting node
A = adjacency(G); % eventually with non-unit and non-negative weight
n = size(A,1);
% initialize distance matrix
d = inf(n,1);
du = 0;
i = start;
while true
    d(i) = du;
    [i,j,dij] = find(A(:,i));
    if isempty(i)
        break
    end
    [du, p] = sort(du(j)+dij);
    [i, p]  = sort(i(p)); % here we requires stable sorting, which is the case with MATLAB
    b = [true; diff(i,1,1)>0];
    i = i(b);
    du = du(p(b));
    b = du < d(i);
    i = i(b);
    du = du(b);
end
% Exit the while-loop d contains distance vector from start-node to all nodes

Answer 3

[Bruno Luong]

Assuming A is the transposed of the adjacent matrix (n x n).

j in (1:n) is the source node number.

This method should be faster than computing A^k proposed by Matt

x = zeros(n,1);
x(j) = 1;
for i=1:k
    x = A*x;
end
i = find(x)';

fprintf('%d-neighbor of %d is %s\n', k, j, mat2str(i))

If sparse form of the adjadcent matrix is readily available, it could be preferable to use it.

[Rub Ron]

Same here:

s = [1 1 1 1 2 2 2 2 2 8 8 12];
t = [3 5 4 2 6 10 7 9 8 11 12 13];
G = graph(s,t);
plot(G)

k=3;
A = adjacency(G);
n = size(A,1);
x = zeros(n,1);
x(1) = 1;
for i=1:k
    x = A*x;
end
i = find(x)'

produces:

i =

     2     3     4     5    11    12

it should be:[ 11 12]

[Bruno Luong]

Seems like you have a directed graph

s = [1 1 1 1 2 2 2 2 2 8 8 12];
t = [3 5 4 2 6 10 7 9 8 11 12 13];
G = digraph(s,t);
plot(G)

k=3;
A = adjacency(G)';
n = size(A,1);
x = zeros(n,1);
x(1) = 1;
for i=1:k
    x = A*x;
end
i = find(x)'

[Rub Ron]

No, unfortunately the graph is not directed. Do you think for directed graph the solution is close to this or might be totall different?

[Bruno Luong]

Have you tried my code with digraph? It gives [11 12] as you claim as the result that is correct.

If the graph is not directed then

[2 3 4 5 11 12] is correct answer. For example

1 => 2 => 1 => 2 

is clear three steps (edge) to reach target 2 from source 1. So 2 must be 3-neighbor of 1.

Take your pick.

[Rub Ron]

Sorry, it is not what I was looking for. I eddited my post to make it clearer.

[Bruno Luong]

When you allow to go from s to t, and forbid to go in the opposite direction from t to s, it's called digraph

s = [1 1 1 1 2 2 2 2 2 8 8 12];
t = [3 5 4 2 6 10 7 9 8 11 12 13];
G = digraph(s,t);
plot(G)
A = adjacency(G).';
n = size(A,1);
x = zeros(n,1);
x(1) = 1;
for k=1:5
    x = A*x;
    i = find(x)';
    fprintf('k = %d: i = %s\n', k, mat2str(i));
end

Output match your expectation

k = 1: i = [2 3 4 5]
k = 2: i = [6 7 8 9 10]
k = 3: i = [11 12]
k = 4: i = 13
k = 5: i = zeros(1,0)

[Rub Ron]

You might be focusing only in the example rather than in a general case. In my case the graphs are not directed, I cannot make them "nicely" directed. For intance, by changing the example graph to this one, the outputs are not the desired ones

s = [3 5 4 2 6 10 7 9 8 11 12 13];
t = [1 1 1 1 2 2 2 2 2 8 8 12];
G = digraph(s,t);

Answer 4

[Bruno Luong]

Given A an (n x n) symmetric ajadcent matrix (undirected graph), and k the order (integer scalar), I propose this algorithm to find index of all k-neigbor of node #start

n = size(A,1);
k = min(k,n);
notdone = true(n,1);
i = start;
for q=1:k
    notdone(i) = false;
    [i,~] = find(A(:,i));
    if isempty(i)
        break
    end
    i = unique(i);
    i = i(notdone(i));
end
neighbouridx = i; % column vector