Computing the local clustering coefficient for a directed graph

Question

I am trying to replicate Mathematica's LocalClusteringCoefficient function (as well as compute $C_i^{out}$) for a directed graph. Which was presumed that Mathematica calculates as the cyclic clustering coefficient, $C_i^{cyc}$ in this post: What exactly does LocalClusteringCoefficient compute for directed graphs? The equation comes from page 13 of this paper, also linked in Szabolcs's post:

$$C_i^{cyl}=$$ $$(A)_{ii}^3\over d_i^{in}d_i^{out}-d_i^{\leftarrow\rightarrow}$$ Where:
$$d_i^{in} = (A^T)_i 1$$ $$d_i^{out} = (A)_i 1$$ $$d_i^{\leftarrow\rightarrow}=(A)_{ii}^2$$

A = Adjacency Matrix, $(A)_i$ = $i$-th row of A, 1 = N-dimensional column vector $(1,1,....1,)^T$

For a graph g:

g = Graph[{1, 2, 3, 4}, {1 -> 2, 2 -> 3, 3 -> 1, 1 -> 4}, VertexLabels ->"Name"]

Here is what I've tried:

s = Range[VertexCount[g]];
one = Table[1, Max[VertexCount[g]]];

din = Transpose[AdjacencyMatrix[g]][[#]] & /@ s;
dout = AdjacencyMatrix[g][[#]] & /@ s;
dlink = b2.b2.b2

numer = dout.dout.dout
denom = (din.one).(dout.one) - dlink;
numer/denom
%.one

With these results:

(* {1/3, 2/3, 1/3, 0} *)

However, LocalClusteringCoefficient[g] yields:

(*{1/2, 1, 1, 0}*)

Answer 1

Your code is incorrect. Here is a direct implementation:

This is the adjacency matrix:

am = AdjacencyMatrix[g];

$(A^3)_{ii}$ in the formula refers to the $i$th element of the diagonal of its cube:

numerator = Normal@Diagonal@MatrixPower[am, 3]

You can get the in- and out-degrees of the vertices using VertexInDegree and VertexOutDegree.

We also need the degrees of reciprocal connections. A simple way to get these is to first construct an adjacency matrix of reciprocal connections, then sum up its rows:

Total[am Transpose[am]]

Then we are ready to compute the denominator of the formula:

denominator = VertexOutDegree[g] VertexInDegree[g] - Total[am Transpose[am]]

Finally just take the ratio:

numerator/denominator
(* {1/2, 1, 1, Indeterminate} *)

The indeterminate comes from 0/0. One way to avoid it is to use

div[0,0] = 0;
div[x_,y_] := Divide[x,y]

MapThread[div, {numerator, denominator}]

You can do it in many other ways of course.

---

Your code:

It would help if next time you commented the code, and indicated what each line is meant to do, like I did above with my code.

s = Range[VertexCount[g]];

one = Table[1, Max[VertexCount[g]]];

I do not understand Max[VertexCount[g]]. VertexCount[g] is a number. What is the point of Max?

din = Transpose[AdjacencyMatrix[g]][[#]] & /@ s

Why are you defining a function and mapping it? The result is simply `Transpose@Adja

dout = AdjacencyMatrix[g][[#]] & /@ s;

Same comment as above. Also, it is confusing to me to see the adjacency matrix and its transpose be called dout and din. These names suggest degree vectors.

dlink = b2.b2.b2

b2 is not defined, so I can't comment on this.

numer = dout.dout.dout

The result is a matrix. The formula asks for its diagonal. You are not taking the diagonal.

denom = (din.one).(dout.one) - dlink;

din.one is a weird way to obtain in-degrees. If you need to do summation, use Total on the adjacency matrix.

(din.one).(dout.one) is the dot product of the degree vectors, and the result is a number. The formula has an element-wise product and the result should be a vector, not a number.

(In case this is what confused you: there is no indication of Einstein summation convention being used in this paper.)