Uniform generation of random bipartite bi-regular graphs?

The problem of generating random $d$-regular graphs uniformly at random has been extensively studied. Some of the proposed algorithms are quite sophisticated. In general, the problem becomes more difficult as $d$ grows as a function of $n$. To the best of my knowledge, the problem is still open for $d = omegaleft(sqrt{n}right)$.

However, if one assumes that $d$ is a constant, then there is a simple and efficient algorithm called the configuration model:

https://en.wikipedia.org/wiki/Configuration_model

It can be easily adapted to the case of biregular graphs with $M$ vertices of degree $k$ on the left side, and $N$ vertices of degree $d$ on the right side. Of course, for this to be possible at all we need $Mk=Nd$.

The algorithm is as follows:

• Assign each vertex on the left side $k$ ``half-edges", and each vertex on the right $d$ half-edges.

• Choose a uniformly random perfect matching of size $Mk$ matching all of the half-edges on the left side to the half-edges on the right side.

• This tells you which vertices to connect to which, but there may be double edges. Repeat the above procedure until there are no double edges.

It is pretty easy to see that the resulting distribution is uniform by counting the number of perfect matching that can yield a given graph. It is harder to analyze the probability that there are no double edges, which is critical if we want to bound the expected runtime. However, it turns out that for constant $k$ and $d$, the probability that there are no double edges is a constant independent of $n$, and so we only need to repeat the algorithm $O(1)$ times.

Because it seems nobody will answer this question any time soon, I will give you what I tried for the same problem. I guess my answer lacks a formal proof but it successfully solved my case. The following python style pseudocode generate edges which contains tuples of vertices.

# Assuming there are M a_nodes and N b_nodes
cur_b_nodes = set(range(N))  # b nodes with the lowest degrees
next_b_nodes = set()  # b nodes with the second lowset degrees

edges = []
for a_node in range(M):
    for _ in range(k):
        if not cur_b_nodes:
            cur_b_nodes = next_b_nodes
            next_b_nodes = set()
        b_node = sample(cur_out_nodes, 1)[0]  # Random sampling
        cur_b_nodes.remove(out_node)
        next_b_nodes.add(out_node)

        edges.append((a_node, b_node))

I think this is an approach with the lowest time complexity. Any comments will be welcome.