equilibrium.jl

using LinearAlgebra: qr

function equilibrium(Q::AbstractMatrix)
    @assert size(Q, 1) == size(Q, 2)
    R = qr(Q').R
    n = size(Q, 1)
    v = zeros(eltype(Q), n)
    v[end] = 1
    for i in n-1:-1:1
        s = zero(eltype(Q))
        for j in i+1:n
            s += R[i,j] * v[j]
        end
        v[i] = -s / R[i,i]
    end
    return v ./ sum(v)
end