Fastest Mixing Markov Chain

Introduction

This example is derived from the results in Boyd, Diaconis, and Xiao (2004), section 2. Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a connected graph with vertices $\mathcal{V}=\{1,\dots ,n\}$ and edges $\mathcal{E}\subseteq \mathcal{V}\times \mathcal{V}$. Assume that $(i,i)\in \mathcal{E}$ for all $i=1,\dots ,n$, and $(i,j)\in \mathcal{E}$ implies $(j,i)\in \mathcal{E}$. Under these conditions, a discrete-time Markov chain on $\mathcal{V}$ will have the uniform distribution as one of its equilibrium distributions. We are interested in finding the Markov chain, constructing the transition probability matrix $P\in {\mathbf{R}}_{+}^{n\times n}$, that minimizes its asymptotic convergence rate to the uniform distribution. This is an important problem in Markov chain Monte Carlo (MCMC) simulations, as it directly affects the sampling efficiency of an algorithm.

The asymptotic rate of convergence is determined by the second largest eigenvalue of $P$, which in our case is $\mu (P):={\lambda }_{max}(P-\frac{1}{n}\mathbf{1}{\mathbf{1}}^T)$ where ${\lambda }_{max}(A)$ denotes the maximum eigenvalue of $A$. As $\mu (P)$ decreases, the mixing rate increases and the Markov chain converges faster to equilibrium. Thus, our optimization problem is

$$\begin{array}{ll}\underset{P}{\text{minimize}}& {\lambda }_{max}(P-\frac{1}{n}\mathbf{1}{\mathbf{1}}^T)\\ \text{subject to}& P\ge 0,P\mathbf{1}=\mathbf{1},P=P^T\\ & P_{ij}=0,(i,j)\notin \mathcal{E}.\end{array}$$

The element $P_{ij}$ of our transition matrix is the probability of moving from state $i$ to state $j$. Our assumptions imply that $P$ is nonnegative, symmetric, and doubly stochastic. The last constraint ensures transitions do not occur between unconnected vertices.

The function ${\lambda }_{max}$ is convex, so this problem is solvable in CVXR. For instance, the code for the Markov chain in Figure 2 below (the triangle plus one edge) is

P <- Variable(c(n,n))
ones <- matrix(1, nrow = n, ncol = 1)

obj <- Minimize(lambda_max(P - 1/n))
constr1 <- list(P >= 0, P %*% ones == ones, P == t(P))
constr2 <- list(P[1,3] == 0, P[1,4] == 0)
prob <- Problem(obj, c(constr1, constr2))
result <- psolve(prob)

where we have set $n=4$. We could also have specified $P\mathbf{1}=\mathbf{1}$ with sum_entries(P,1) == 1, which uses the sum_entries atom to represent the row sums.

Example

In order to reproduce some of the examples from Boyd, Diaconis, and Xiao (2004), we create functions to build up the graph, solve the optimization problem and finally display the chain graphically.

## Boyd, Diaconis, and Xiao. SIAM Rev. 46 (2004) pgs. 667-689 at pg. 672
## Form the complementary graph
antiadjacency <- function(g) {
    n <- max(as.numeric(names(g)))   ## Assumes names are integers starting from 1
    setNames(lapply(as.character(1:n), function(v) {
        setdiff(1:n, g[[v]])
    }), 1:n)
}

## Fastest mixing Markov chain on graph g
FMMC <- function(g, verbose = FALSE) {
    a <- antiadjacency(g)
    n <- length(names(a))
    P <- Variable(c(n, n))
    o <- rep(1, n)
    objective <- Minimize(cvxr_norm(P - 1.0/n, "2"))
    constraints <- list(P %*% o == o, t(P) == P, P >= 0)
    ## Zero out transition probabilities for non-edges (excluding diagonal)
    non_edge_constraints <- do.call(c, lapply(names(a), function(i) {
        idx <- as.numeric(i)
        js <- setdiff(a[[i]], idx)
        lapply(js, function(j) P[idx, j] == 0)
    }))
    constraints <- c(constraints, non_edge_constraints)
    prob <- Problem(objective, constraints)
    opt_val <- psolve(prob)
    if(verbose)
        cat("Status: ", status(prob), ", Optimal Value = ", opt_val, "\n")
    list(status = status(prob), value = opt_val, P = value(P))
}

disp_result <- function(states, P, tolerance = 1e-3) {
    if(!("markovchain" %in% rownames(installed.packages()))) {
        rownames(P) <- states
        colnames(P) <- states
        print(P)
    } else {
        P[P < tolerance] <- 0
        P <- P/apply(P, 1, sum)   ## Normalize so rows sum to exactly 1
        mc <- new("markovchain", states = states, transitionMatrix = P)
        plot(mc)
    }
}

Results

Table 1 from Boyd, Diaconis, and Xiao (2004) is reproduced below.

We reproduce the results for various rows of the table.

g <- list("1" = 2, "2" = c(1,3), "3" = c(2,4), "4" = 3)
result1 <- FMMC(g, verbose = TRUE)

Status:  optimal , Optimal Value =  0.8451543 

disp_result(names(g), result1$P)

Row 1, line graph

g <- list("1" = 2, "2" = c(1,3,4), "3" = c(2,4), "4" = c(2,3))
result2 <- FMMC(g, verbose = TRUE)

Status:  optimal , Optimal Value =  0.7071068 

disp_result(names(g), result2$P)

Row 2, triangle plus one edge

g <- list("1" = c(2,4,5), "2" = c(1,3), "3" = c(2,4,5), "4" = c(1,3), "5" = c(1,3))
result3 <- FMMC(g, verbose = TRUE)

Status:  optimal , Optimal Value =  0.7559289 

disp_result(names(g), result3$P)

Bipartite 2 plus 3

g <- list("1" = c(2,3,5), "2" = c(1,4,5), "3" = c(1,4,5), "4" = c(2,3,5), "5" = c(1,2,3,4,5))
result4 <- FMMC(g, verbose = TRUE)

Status:  optimal , Optimal Value =  0.4588315 

disp_result(names(g), result4$P)

Square plus central point

Extensions

It is easy to extend this example to other Markov chains. To change the number of vertices, we would simply modify n, and to add or remove edges, we need only alter the constraints in constr2. For instance, the bipartite chain in Figure 3 is produced by setting $n=5$ and

constr2 <- list(P[1,3] == 0, P[2,4] == 0, P[2,5] == 0, P[4,5] == 0)

Session Info

R version 4.5.2 (2025-10-31)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.3

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Los_Angeles
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] markovchain_0.10.3 Matrix_1.7-4       CVXR_1.8.1        

loaded via a namespace (and not attached):
 [1] piqp_0.6.2            expm_1.0-0            jsonlite_2.0.0       
 [4] compiler_4.5.2        highs_1.12.0-3        Rcpp_1.1.1           
 [7] slam_0.1-55           parallel_4.5.2        cccp_0.3-3           
[10] yaml_2.3.12           fastmap_1.2.0         clarabel_0.11.2      
[13] lattice_0.22-9        igraph_2.2.2          knitr_1.51           
[16] htmlwidgets_1.6.4     backports_1.5.0       Rcplex_0.3-8         
[19] checkmate_2.3.4       gurobi_13.0-1         osqp_1.0.0           
[22] rlang_1.1.7           xfun_0.56             S7_0.2.1             
[25] RcppParallel_5.1.11-1 otel_0.2.0            cli_3.6.5            
[28] magrittr_2.0.4        Rglpk_0.6-5.1         digest_0.6.39        
[31] grid_4.5.2            gmp_0.7-5.1           lifecycle_1.0.5      
[34] ECOSolveR_0.6.1       scs_3.2.7             evaluate_1.0.5       
[37] codetools_0.2-20      Rmosek_11.1.1         stats4_4.5.2         
[40] rmarkdown_2.30        tools_4.5.2           pkgconfig_2.0.3      
[43] htmltools_0.5.9      

References

Boyd, S., P. Diaconis, and L. Xiao. 2004. “Fastest Mixing Markov Chain on a Graph.” SIAM Review 46 (4): 667–89.