An efficient way to use a markov chain to detect how cyclical a graph is?

Question

I've been attempting to generate continuous adjacency matrices for undirected acyclic graphs (ideally trees) using a neural network, and I was thinking that I could solve this problem by creating a cost function with a differentiable, continuous score of how cyclical a graph is given its adjacency matrix. (If anyone has good ideas for how to solve this problem, I'm also interested in non markov chain-based solutions!) The solution I've come up works correctly (for small examples only), giving a score of 0 for trees and (0,1) scores for random dense adjacency matrices. A first order markov chain wouldn't work because on an undirected graph, it will take a step and then return the way it came, and ideally to detect cycles, you reach the node you started at without cheating :)

So far the only solution I've come up with is a second order markov chain representing a progressive random walk, where you can get the probability that a node can be returned to (without ever directly returning the direction you came from).

Here's my exploding implementation:

import pytorch

def expand_adj_matrix(adj_matrix):
    n = adj_matrix.shape[0]
    expanded_matrix = torch.zeros((n**2, n**2))

    for i in range(n):
        for j in range(n):
            if adj_matrix[i][j] > 0:  # Check if there is an edge
                for k in range(n):
                    if k != i:  # Avoiding immediate reversal
                        # Use the weight of the edge from the original adjacency matrix
                        expanded_matrix[n*i + j][n*j + k] = adj_matrix[j][k]
    return expanded_matrix

def cyclical_score_nobackforth(adj_matrix, max_power=9, eps=1e-6):
    expanded_matrix = expand_adj_matrix(adj_matrix)
    # Normalizing the expanded matrix
    norm_expanded_matrix = expanded_matrix / (expanded_matrix.sum(1, keepdim=True) + eps)

    score = 0
    for k in range(2, max_power + 1):
        # Raising the matrix to power k
        matrix_power = torch.matrix_power(norm_expanded_matrix, k)
        trace = matrix_power.trace()
        score += trace

    return score / max_power


#Make a random adjacency matrix and score it
N=20
y1=torch.rand((50,N))
adj_matrix=y1.T@y1
for i in range(N):
    adj_matrix[i,i]=0

norm_adj_matrix=torch.tensor(adj_matrix/adj_matrix.sum(0)).T.float()
cyclical_score_nobackforth(norm_adj_matrix,5)

But the obvious problem is that the second order markov chain explodes in size as it scales by the number of nodes squared squared and I need to generate graphs with hundreds to thousands of nodes! I'm wondering if it would be possible to shrink the size of the second order markov chain, since you really only need to know if the chain is making a move like A-B-A so perhaps things could be summed.

Any thoughts on this problem?