Directed graph similarity calculation -- optimizing Mathematica code

I am looking to optimize the following computation. First, some explanations:

Explanation

I have a directed graph whose adjacency matrix may only contain 1s and 0s. The similarity of nodes i and j is defined as follows:

• count how many nodes have an edge pointing to both i and j, and count how many nodes do both i and j point to

• divide the sum of these two counts by the total number of nodes either i or j is connected to

I need to find the maximum similarity of each node in graph (its similarity with any other node except itself), and compute the average of these numbers.

The graphs I'm working with are not necessarily connected, but they don't have components of size 1 (i.e. nodes that are not connected to any other nodes).

The code

meanMaxSim[g_Graph] :=
 Module[{am, pList, similarityMatrix},
  am = AdjacencyMatrix[g];
  pList = Transpose@Join[am, Transpose[am]];
  similarityMatrix = Outer[Total[#1 #2]/Total[Sign[#1 + #2]] &, pList, pList, 1];
  Mean[Max /@ dropDiagonal[similarityMatrix]]
]

(Sign is used to keep 0s and convert anything larger to 1s.)

This function just removes the elements on the diagonal of an n by n matrix, ending up with an n by n - 1 matrix:

dropDiagonal[matrix_] := 
 MapThread[Drop[#1, {#2}] &, {matrix, Range@Length[matrix]}]

Graph[] is only available in Mathematica 8, but this is irrelevant here as most of the time is spent in Outer[] anyway. If you have an earlier version, just assume that you have the adjacency matrix as a SparseArray as starting point.

I am looking for modifications that speed this up significantly (more than 3x), but any other comment on the code is most welcome. Hopefully a speedup is possible without a significant increase in code complexity.

Notes

Doing the computations with machine-numbers is okay, but I didn't bother doing that since it doesn't provide a significant speedup.

A note about self-loops: they are allowed in the graph, and yes, they're double-counted in this

Minor improvements, but consider using:

pList = ArrayFlatten[{{am\[Transpose], am}}];

and

dropDiagonal = MapThread[Delete, {#, Range@Length@#}] &;

I also tried dropDiagonal = MapIndexed[Drop, #] & but this is slower, at least in v7.

Replace:

Total[#1 #2]/Total[Sign[#1 + #2]] &

With:

#.#2 / Tr@Unitize[#1 + #2] &

I believe this is an improvement to your compiled function:

Compile[{{p, _Real, 2}}, 
  With[{b = p\[Transpose]},
    #.b / Total@Sign[# + b] & /@ p
  ]
]

StackExchange Code Review Q#2385, answer score: 5