Speed optimization for girth-counting algorithm

Question

I have the following girth-counting algorithm on a large 0-1 matrix and I need to optimize it for speed:

function CountG8()
M = 512
N = 1024
H = rand(0:1,M,N)
girth_8 = 0

for i=1:M-3
    for j=1:N-1
        if H[i,j]==1
            for j1=j+1:N
                if H[i,j1]==1
                    for i1=i+1:M
                        if H[i1,j1]==1
                            for j2=1:N
                                if H[i1,j2]==1&&j2!=j1
                                    for i2=i+1:M
                                        if H[i2,j2]==1&&i2!=i1
                                            for j3=1:N
                                                if H[i2,j3]==1&&j3!=j2
                                                    for i3=i+1:M
                                                        if H[i3,j3]==1&&i3!=i2
                                                            if H[i3,j]==1
                                                                girth_8 += 1
                                                            end
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end
end

The algorithm runs without errors, but it is too slow.

The code is supposed to count every cycle of length 8 in a bipartite graph.

A bipartite graph is a graph of two groups of nodes. Left group, Check Nodes (CNs) and right group, Variable Nodes (VNs). Every CN from left group is allowed to be connected to a VN from right group. This is known as Directed Graph (DG). A DG has only cycles of even length 4,6,8,10,.. The Girth of the graph is the length of shortest cycle in it.

The algorithm walks through the graph in a BFS fashion searching for all possible nodes combination forming a cycle of the specified length (here 8). Matrix H in the code represents the graph. Each row is a CN and each column is a VN in the graph.

Here is a Bipartite Graph with its corresponding matrix as well as 2 cycles of length 4 and 6 shown:

EDIT: A comment on matrix power approach given in the answer.

In my case, the adjacency matrix is (M+N)x(M+N) = 1536 x 1536. Watch out for the complexity of (1536x1536)^8.

Answer 1

What you should be using is a neighboring matrix.

The key is graph theory.

You represent the graph with a matrix.

For example:
\$ \begin{bmatrix} 0 & 1 & 1 & 0 & 2 \\ 1 & 0 & 2 & 0 & 0 \\ 1 & 2 & 0 & 1 & 1 \\ 0 & 0 & 1 & 0 & 1 \\ 2 & 0 & 1 & 1 & 0 \\ \end{bmatrix} \$
Here is a graph and its corresponding neighbouring matrix. If you look at the first row you see that A has no one-step connection to itself. But it has a one-step connection to every other node. Node B has a one-step connection to A and 2 different one-step connections to node C. As you can see it all sums up if you keep going like this.

Now then, here's the mathemagic. If you multiply the matrix with itself the result will be a matrix which tells you all the connections of two-stop length. Continue multiplying with the original matrix until you've done it 8 times and you will have a matrix which tells you how many connections of eight-step length you have.

How is matrix multiplication done, you might ask?
Wikipedia
Another site (with examples)
Basically, you multiply rows with columns.(It would take another whole answer to explain it)

When you've gotten the neighbouring matrix of the eight power, you just sum up all the elements of the diagonal. That will tell you the number of cycles of length 8. (because the diagonal holds the information of a n-step way back to itself, where n is the power of which the matrix has been multiplied)

This method works for all graphs. All you need to do is make sure the graph is bipartite to begin with.