Chapter 8 Experiment 8: Mixing Times of Markov Chains and Matrix Structure

What Is This Experiment About?

Imagine taking a random walk on a map: at every step, where you go next is decided by a “transition rule.” After enough steps, the frequencies with which you appear at the various locations settle down—this is the stationary distribution.

The question is: how many steps does it take to “settle down”?

The answer is hidden in the eigenvalues of a matrix. Through Markov chains with five different structures, this experiment helps you build three core intuitions:

The shape of the matrix determines the speed of convergence: different transition rules correspond to different arrangements of the eigenvalues in the complex plane, and this “fingerprint” tells you directly how many steps are needed to settle down
The larger the spectral gap, the faster the convergence: $\text{gap} = 1 - |\lambda_1|$ measures “how far the second-largest eigenvalue is from 1”; the farther it is, the faster the chain settles down
What happens when random perturbations are added: if every step carries a little noise, how does the convergence behavior change?

How to Use This Experiment¶

This experiment consists of two interactive modules:

Module 1 (§1): a unified dashboard in which you can switch among the five chains; seven panels update together:

Top left: the positions of the eigenvalues in the complex plane (the “fingerprint plot”)
Top middle: the evolution of the heat map of the matrix power $\mathbf{P}^k$ (the “convergence animation”)
Bottom left: the decay curve of the TV distance (the “mixing-speed curve”)
Bottom row, second panel: the stationary distribution $\boldsymbol{\pi}$

Module 2 (§7): the perturbation experiment, comparing two models: a single perturbation vs injecting randomness at every step

Each of §2–§6 has a “concrete matrix” cell that uses a small matrix to show you clearly the structure of each chain.

Structure of the Experiment¶

Section	Topic	In one sentence
§0	Theoretical framework	Mathematical definitions of the stationary distribution, the mixing time, and the TV distance
§1	Unified dashboard	Switch among the five chains: eigenvalues + heat maps + TV curve
§2	Birth-death chain	Moves only one step at a time, like climbing stairs
§3	Cycle chain	Walks around a ring; even-length rings have a periodicity problem
§4	Google Matrix	Web page ranking: irregular links + teleportation
§5	Doubly stochastic matrix	Row and column sums all equal 1, so the stationary distribution must be uniform
§6	Hypercube	A random walk that flips bits; the higher the dimension, the slower
§7	Perturbation experiment	How the convergence behavior changes when noise is added

§0 Common Theoretical Framework¶

Before turning to the concrete examples, we first make a few key concepts clear. The material in this section is the “common language” of all the experiments that follow.

§0.1 What Is a Transition Matrix?¶

Suppose a system has $N$ possible states (for example, the weather may be “sunny,” “rainy,” or “cloudy”). The entry $P_{ij}$ of the transition matrix $\mathbf{P}$ is “the probability of jumping to state $j$ at the next step, given that the system is in state $i$ .”

Since a system that starts from state $i$ must go to some state, all these probabilities must add up to 1:

\sum_{j=0}^{N-1} P_{ij} = 1, \quad P_{ij} \geq 0

(1)

This is called a transition matrix (stochastic matrix). Note that it is each row whose sum is 1, because the first subscript $i$ of $P_{ij}$ is the “point of departure.”

Evolution of the state distribution: if the current state distribution is the row vector $\boldsymbol{\mu}^{(k)}$ , then after one step

\boldsymbol{\mu}^{(k+1)} = \boldsymbol{\mu}^{(k)} \mathbf{P}

(2)

and after $k$ steps

\boldsymbol{\mu}^{(k)} = \boldsymbol{\mu}^{(0)} \mathbf{P}^k

(3)

§0.2 What Is a Stationary Distribution?¶

After enough steps, the distribution no longer changes—this is the stationary distribution $\boldsymbol{\pi}$ :

\boldsymbol{\pi}^\top \mathbf{P} = \boldsymbol{\pi}^\top

(4)

This equation says that $\boldsymbol{\pi}$ multiplied by $\mathbf{P}$ is still $\boldsymbol{\pi}$ ; that is, $\boldsymbol{\pi}$ is an eigenvector of $\mathbf{P}^\top$ for the eigenvalue $\lambda_0 = 1$ .

In the language of linear systems:

(\mathbf{P}^\top - \mathbf{I})\boldsymbol{\pi} = \mathbf{0}

(5)

That is, we find the null space of $\mathbf{P}^\top - \mathbf{I}$ ; adding the normalization condition $\sum_i \pi_i = 1$ then determines $\boldsymbol{\pi}$ uniquely (for an irreducible chain).

§0.3 Eigenvalues and the Speed of Convergence¶

For an irreducible aperiodic chain, the eigenvalues of $\mathbf{P}$ are ordered as follows:

\underbrace{1}_{ \lambda_0 } > |\lambda_1| \geq |\lambda_2| \geq \cdots \geq |\lambda_{N-1}|

(6)

Diagonalizing $\mathbf{P}$ : $\mathbf{P}^k = \mathbf{1}\boldsymbol{\pi}^\top + \sum_{i=1}^{N-1} \lambda_i^k \mathbf{v}_i \mathbf{u}_i^\top$

When $k$ is large, $\lambda_i^k \to 0$ (because $|\lambda_i| < 1$ ), so every row of $\mathbf{P}^k$ tends to $\boldsymbol{\pi}^\top$ .

The speed of convergence is determined by $|\lambda_1|$ : the closer $|\lambda_1|$ is to 1, the more slowly $\lambda_1^k$ decays, and the more steps are needed to converge.

§0.4 The Spectral Gap and the Mixing Time¶

The spectral gap is the distance between $\lambda_0$ and $\lambda_1$ :

\text{gap} = 1 - |\lambda_1| \in (0, 1]

(7)

The larger the spectral gap → the farther the second-largest eigenvalue $|\lambda_1|$ is from 1 → the faster the convergence.

The total variation distance (TV distance) quantifies “how far the current distribution is from the stationary distribution”:

\|\boldsymbol{\mu}^{(k)} - \boldsymbol{\pi}\|_{\mathrm{TV}} = \frac{1}{2}\sum_{i=0}^{N-1}|\mu^{(k)}_i - \pi_i|

(8)

Intuitively: if you regard $\boldsymbol{\mu}^{(k)}$ and $\boldsymbol{\pi}$ as two probability distributions, the TV distance is “the total discrepancy between the two distributions”; its range is $[0, 1]$ , and 0 means they are identical.

The mixing time is the smallest number of steps needed for the TV distance to drop to the threshold $\varepsilon = 0.25$ :

t_{\mathrm{mix}} = \min\bigl\{k : \|\boldsymbol{\mu}^{(k)} - \boldsymbol{\pi}\|_{\mathrm{TV}} \leq 0.25\bigr\} \approx \frac{\log(1/2\varepsilon)}{\text{gap}} = \frac{\log 2}{\text{gap}}

(9)

This formula tells you that if the spectral gap shrinks by half, the mixing time roughly doubles.

# Colab: enable the custom widget manager so that FigureWidget and the sliders are interactive; skipped elsewhere
try:
    from google.colab import output
    output.enable_custom_widget_manager()
except ImportError:
    pass

# ============================================================
# Imports and global settings
# ============================================================
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import ipywidgets as widgets
from IPython.display import display

np.set_printoptions(precision=4, suppress=True, linewidth=100)

# --- Color scheme (Accent Mix) ---
C_BG    = "#F8F8F8"
C_GRID  = "#D6D6D6"
C_AXIS  = "#000000"
C_V1    = "#57068C"   # violet
C_V2    = "#006385"   # deep blue
C_T1    = "#2AD2C9"   # teal
C_T2    = "#8900E1"   # Ultra Violet
C_WARN  = "#FF5D47"   # orange
C_AUX   = "#AB82C5"   # light violet

# Colors for the chains (one per chain)
CHAIN_COLORS = {
    'birth':    C_V1,
    'cycle':    C_V2,
    'google':   C_WARN,
    'dbstoch':  C_T1,
    'hypercube':C_T2,
}

def print_header(title):
    print("=" * 60)
    print(f"  {title}")
    print("=" * 60)

def print_step(step, desc):
    print(f"\n▶ Step {step}: {desc}")
    print("-" * 40)

def compare_print(label, actual, expected_desc):
    print(f"\n[{label}]")
    print(f"  Computed value: {actual}")
    print(f"  Expected value: {expected_desc}")

print_header("Markov chain experiment | initialization complete")
print("  Color scheme: book palette")
print("  Interactive framework: plotly FigureWidget + ipywidgets")

§0.5 Mathematical Utility Functions¶

The following six functions are the computational foundation of the whole experiment; every later cell uses them. Next to each function is a description of the mathematical quantity it computes.

# ============================================================
# Shared mathematical tools
# ============================================================

# --- Step 1: stationary distribution (power iteration) ---
def stationary_dist(P, tol=1e-13, max_iter=5000):
    """Find the stationary distribution π of the transition matrix P by power iteration"""
    n = P.shape[0]
    v = np.ones(n) / n
    for _ in range(max_iter):
        nv = v @ P
        if np.max(np.abs(nv - v)) < tol:
            return nv
        v = nv
    return v

# --- Step 2: second-largest eigenvalue |λ₁| (numpy directly; most stable for N<=30) ---
def lambda1_abs(P):
    """
    Find the absolute value |λ₁| of the second-largest eigenvalue.
    For matrices of teaching size (N<=30), calling numpy.linalg.eigvals directly is the most reliable choice.
    This avoids the failure of a home-made power iteration on periodic chains (|λ₁|=1, e.g., λ=−1).
    """
    eigs = np.abs(np.linalg.eigvals(P))
    eigs_sorted = np.sort(eigs)[::-1]
    return float(eigs_sorted[1]) if len(eigs_sorted) > 1 else 0.0

# --- Step 3: total variation distance ---
def tv_dist(mu, pi):
    return 0.5 * np.sum(np.abs(mu - pi))

# --- Step 4: mixing-time estimate (handles periodic chains with |λ₁|≥1) ---
def mixing_time(lam, eps=0.25):
    """
    t_mix = ⌈log(1/2ε) / log(1/|λ₁|)⌉
    If |λ₁| >= 1 (a periodic chain: the cycle chain/hypercube with r=0),
    return None, meaning no convergence under this definition (a lazy probability is needed to break the periodicity).
    """
    if lam >= 1.0 - 1e-6:
        return None
    if lam < 1e-9:
        return 1
    return int(np.ceil(np.log(1 / (2 * eps)) / np.log(1 / lam)))

# --- Step 5: matrix power (repeated squaring) ---
def mat_pow(P, k):
    if k == 0:
        return np.eye(P.shape[0])
    if k == 1:
        return P.copy()
    half = mat_pow(P, k // 2)
    R = half @ half
    return R @ P if k % 2 else R

# --- Step 6: decay sequence of the TV distance (built-in cap so periodic chains do not blow up) ---
def tv_decay_series(P, pi, max_k):
    """
    Starting from δ₀ (all probability in state 0), compute the sequence of TV distances over max_k steps.
    max_k is truncated at 500 so that periodic chains do not make the computation too slow.
    """
    max_k = min(int(max_k), 500)
    mu = np.zeros(P.shape[0]); mu[0] = 1.0
    tvs = []
    for _ in range(max_k + 1):
        tvs.append(float(tv_dist(mu, pi)))
        mu = mu @ P
    return np.array(tvs)

# ============================================================
print_header("Utility functions | loaded")
print("  stationary_dist : stationary π by power iteration")
print("  lambda1_abs(P)      : |λ₁| via numpy eigvals (pi not needed)")
print("  tv_dist         : total variation distance ‖μ − π‖_TV")
print("  mixing_time     : t_mix; returns None for periodic chains")
print("  mat_pow         : P^k by repeated squaring")
print("  tv_decay_series : TV decay sequence (capped at 500 steps)")

§0.6 Constructing the Five Transition Matrices¶

Below we define the transition matrices of five well-known Markov chains. Behind each function lies a different physical or mathematical motivation:

Function	Model	An everyday analogy
`make_birth_death`	Birth-death chain	Random rises and falls in the length of a queue
`make_cycle`	Cycle chain	Walking randomly around a ring
`make_google`	Google Matrix	Browsing randomly among web pages
`make_doubly_stochastic`	Doubly stochastic matrix	Balanced flows in a logistics network
`make_hypercube`	Hypercube walk	Randomly flipping one bit of a binary password

# ============================================================
# Constructing the transition matrices of the five Markov chains
# ============================================================

def make_birth_death(S=8, p=0.5):
    """
    Birth-death chain
    State space: {0, 1, …, S}; matrix size (S+1)×(S+1)
    P[i][i+1] = p (to the right), P[i][i-1] = q = 1-p (to the left)
    Boundary: P[0][0]=q, P[0][1]=p; P[S][S-1]=q, P[S][S]=p
    """
    n = S + 1
    q = 1 - p
    P = np.zeros((n, n))
    P[0, 0] = q; P[0, 1] = p
    P[n-1, n-2] = q; P[n-1, n-1] = p
    for i in range(1, n-1):
        P[i, i-1] = q
        P[i, i+1] = p
    return P

def make_cycle(N=8, r=0.0):
    """
    Cycle chain
    Symmetric random walk on the ring C_N
    P[i][(i±1) mod N] = (1-r)/2, P[i][i] = r (lazy probability)
    Eigenvalues: λ_k = r + (1-r)cos(2πk/N), k = 0,1,…,N-1
    """
    P = np.zeros((N, N))
    prob = (1 - r) / 2
    for i in range(N):
        P[i, i] = r
        P[i, (i + 1) % N] = prob
        P[i, (i - 1) % N] = prob
    return P

def make_google(N=8, alpha=0.85, seed=42):
    """
    Google Matrix (PageRank)
    G = α·P_link + (1-α)·(1/N)·11^T
    P_link: Barabási–Albert preferential-attachment directed graph → row-stochastic matrix
      - each new node chooses its link targets with probability proportional to the in-degrees of the existing nodes
      - dangling nodes (out-degree = 0) jump uniformly to all nodes instead
    Rank-one correction (1-α)·(1/N)·11^T:
      - makes every entry of the matrix > 0 (irreducible)
      - guarantees a unique stationary distribution (the PageRank vector)
      - |λ_k| ≤ α for k ≥ 1 (lower bound on the spectral gap = 1-α)
    """
    rng = np.random.default_rng(seed)
    m = max(1, min(3, N // 4))   # each new node links to m existing nodes
    adj = np.zeros((N, N), dtype=float)
    # initially: m+1 nodes all linked to one another
    for ii in range(m + 1):
        for jj in range(m + 1):
            if ii != jj:
                adj[ii, jj] = 1.0
    in_deg = adj.sum(axis=0).copy()
    # BA preferential attachment: add the new nodes one at a time
    for new in range(m + 1, N):
        probs = in_deg[:new].copy()
        if probs.sum() == 0:
            probs = np.ones(new)
        probs = probs / probs.sum()
        targets = rng.choice(new, size=min(m, new),
                              replace=False, p=probs)
        for t in targets:
            adj[new, t] = 1.0
            in_deg[t] += 1
    # normalize each row (dangling nodes get a uniform row)
    for ii in range(N):
        s = adj[ii].sum()
        adj[ii] = adj[ii] / s if s > 0 else np.ones(N) / N
    P_link = adj
    # Google correction
    return alpha * P_link + (1 - alpha) / N * np.ones((N, N))
def make_doubly_stochastic(N=6, m=0.6):
    """
    Doubly stochastic matrix
    All row sums and column sums equal 1; the stationary distribution is always uniform, π = (1/N)1
    Built from the mixing strength m: P[i][i] ≈ m, P[i][(i+1)%N] ≈ m/N
    The remaining entries are filled uniformly so that all row and column sums equal 1
    """
    P = np.full((N, N), (1 - m) / N)
    for i in range(N):
        P[i, i] += m * (1 - 1/N)
        P[i, (i + 1) % N] += m / N
    return P

def make_hypercube(n=3, r=0.0):
    """
    Hypercube walk
    State space: {0,1}ⁿ; size 2ⁿ × 2ⁿ
    At each step one bit, chosen uniformly, is flipped; lazy probability r
    Eigenvalues: λ_k = r + (1-r)(1 - 2k/n), multiplicity C(n,k)
    Mixing time ∼ (n/2)ln n (cutoff phenomenon)
    """
    sz = 1 << n
    P = np.zeros((sz, sz))
    for i in range(sz):
        P[i, i] = r
        for b in range(n):
            P[i, i ^ (1 << b)] += (1 - r) / n
    return P

# ============================================================
print_header("Matrix construction functions | loaded")
chains_info = [
    ("Birth-death", "make_birth_death(S, p)",    "tridiagonal, (S+1)×(S+1)"),
    ("Cycle",     "make_cycle(N, r)",           "circulant matrix, N×N"),
    ("Google",    "make_google(N, alpha)",       "rank-one correction, N×N"),
    ("Doubly stoch.", "make_doubly_stochastic(N,m)","row and column sums = 1, N×N"),
    ("Hypercube",   "make_hypercube(n, r)",        "Kronecker product, 2ⁿ×2ⁿ"),
]
for name, func, note in chains_info:
    print(f"  {name:<10}  {func:<32}  # {note}")

§1 Unified Comparison Dashboard¶

This dashboard lets you compare the five chains with the same “measuring stick.” Switch chains or adjust the parameters, and the seven plots update together.

How to Read the Seven Panels¶

Top row (matrix structure)

Eigenvalues (complex plane): each point is an eigenvalue. The orange star on the far right is always $\lambda_0 = 1$ (the dominant eigenvalue). The closer the other points are to the origin = the faster the mixing; points lying on the unit circle = a periodicity problem.
Heat maps of $\mathbf{P}^k$ (four panels): the shade of each cell of the matrix indicates the size of the transition probability (white = 0, dark violet = high). From left to right, $k = 1, t_{\text{mix}}/3, t_{\text{mix}}, 3t_{\text{mix}}$ . When the colors of all the rows become alike, the chain has converged.

Bottom row (convergence analysis)

TV distance decay: the horizontal axis is the number of steps $k$ , and the vertical axis is the distance from the stationary distribution. The orange dot marks the mixing time $t_{\text{mix}}$ (where the TV distance drops to 0.25). The green dashed line is the threshold $\varepsilon = 0.25$ .
Stationary distribution $\boldsymbol{\pi}$ : the frequency with which each state is eventually visited. A uniform distribution = all states are on an equal footing; a skewed distribution = some states are more “important.”

The Metrics Bar¶

The title bar of the dashboard shows three numbers:

$|\lambda_1|$ : the absolute value of the second-largest eigenvalue; the smaller, the better
Spectral gap $= 1 - |\lambda_1|$ : the larger it is, the faster the convergence
$t_{\text{mix}}$ : the number of steps needed for the TV distance to drop to 0.25; the smaller, the better

# ============================================================
# §1 Unified comparison dashboard
# ============================================================

# --- Step 1: heat-map helper: turn a matrix into a heatmap trace ---
def matrix_heatmap_trace(M, title=""):
    n = M.shape[0]
    # custom color scale: white (0) → violet (1)
    colorscale = [
        [0.0, "#FFFFFF"],
        [0.5, "#AB82C5"],
        [1.0, "#57068C"]
    ]
    return go.Heatmap(
        z=M,
        colorscale=colorscale,
        zmin=0, zmax=1,
        showscale=False,
        xgap=0.5, ygap=0.5,
    )

# --- Step 2: build the dashboard ---
def build_dashboard(P, chain_name, color):
    pi    = stationary_dist(P)
    lam   = lambda1_abs(P)
    gap   = 1 - lam
    tmix  = mixing_time(lam)
    N     = P.shape[0]

    # steps for the heat maps: k=1, t_mix/3, t_mix, 3*t_mix
    _tmix_safe = tmix if tmix is not None else 30
    k_steps = [
        1,
        max(2, _tmix_safe // 3),
        _tmix_safe,
        min(_tmix_safe * 3, 400)
    ]

    # eigenvalues (numpy)
    eigs = np.linalg.eigvals(P)

    # TV decay
    _tm = tmix if tmix is not None else 50
    max_k = min(_tm * 5 + 30, 300)
    tv_series = tv_decay_series(P, pi, max_k)

    # --- Build the subplot layout: 2 rows × 5 columns ---
    # top row: eigenvalue plane | heat maps ×4
    # bottom row: TV decay | stationary distribution
    fig = make_subplots(
        rows=2, cols=5,
        column_widths=[0.22, 0.195, 0.195, 0.195, 0.195],
        row_heights=[0.5, 0.5],
        subplot_titles=(
            "Eigenvalues (complex plane)",
            f"P<sup>{k_steps[0]}</sup>",
            f"P<sup>{k_steps[1]}</sup>",
            f"P<sup>{k_steps[2]}</sup> (t_mix)",
            f"P<sup>{k_steps[3]}</sup>",
            "TV distance decay", "Stationary distribution π", "", "", ""
        ),
        specs=[
            [{"type": "scatter"}, {"type": "heatmap"}, {"type": "heatmap"},
             {"type": "heatmap"}, {"type": "heatmap"}],
            [{"type": "scatter"}, {"type": "bar"}, {"type": "scatter"},
             {"type": "scatter"}, {"type": "scatter"}],
        ]
    )

    # --- Eigenvalue plane (row=1, col=1) ---
    theta_circle = np.linspace(0, 2*np.pi, 200)
    fig.add_trace(go.Scatter(
        x=np.cos(theta_circle), y=np.sin(theta_circle),
        mode='lines',
        line=dict(color=C_GRID, width=1, dash='dot'),
        showlegend=False, hoverinfo='skip'
    ), row=1, col=1)
    # non-dominant eigenvalues
    mask = np.abs(eigs - 1.0) > 0.01
    fig.add_trace(go.Scatter(
        x=eigs[mask].real, y=eigs[mask].imag,
        mode='markers',
        marker=dict(color=color, size=8, opacity=0.8,
                    line=dict(color='white', width=1)),
        name='λ_k (k≥1)', showlegend=False,
        hovertemplate='λ = %{x:.4f} + %{y:.4f}i<extra></extra>'
    ), row=1, col=1)
    # λ_0 = 1
    fig.add_trace(go.Scatter(
        x=[1], y=[0], mode='markers',
        marker=dict(color=C_WARN, size=12, symbol='star',
                    line=dict(color='white', width=1)),
        name='λ₀ = 1', showlegend=False,
        hovertemplate='λ₀ = 1 (stationary)<extra></extra>'
    ), row=1, col=1)

    # --- Heat maps (row=1, col=2,3,4,5) ---
    for idx, k in enumerate(k_steps):
        Pk = mat_pow(P, k)
        fig.add_trace(matrix_heatmap_trace(Pk), row=1, col=idx+2)

    # --- TV decay (row=2, col=1) ---
    ks = np.arange(len(tv_series))
    fig.add_trace(go.Scatter(
        x=ks, y=tv_series,
        mode='lines', line=dict(color=color, width=2),
        name='TV distance', showlegend=False,
        hovertemplate='k=%{x}, TV=%{y:.4f}<extra></extra>'
    ), row=2, col=1)
    fig.add_trace(go.Scatter(
        x=[0, max_k], y=[0.25, 0.25],
        mode='lines', line=dict(color=C_T1, width=1, dash='dash'),
        showlegend=False, hoverinfo='skip'
    ), row=2, col=1)
    # mark t_mix
    fig.add_trace(go.Scatter(
        x=[tmix], y=[0.25],
        mode='markers+text',
        marker=dict(color=C_WARN, size=10),
        text=[f'k={tmix}'], textposition='top right',
        showlegend=False
    ), row=2, col=1)

    # --- Stationary distribution (row=2, col=2) ---
    fig.add_trace(go.Bar(
        x=list(range(N)), y=pi,
        marker_color=color, opacity=0.8,
        showlegend=False,
        hovertemplate='state %{x}: π=%{y:.4f}<extra></extra>'
    ), row=2, col=2)

    # --- Layout settings ---
    fig.update_layout(
        title=dict(
            text=f"<b>{chain_name}</b>&nbsp;&nbsp;|λ₁| = {lam:.4f}&nbsp;&nbsp;spectral gap = {gap:.4f}&nbsp;&nbsp;t_mix = {'∞ (periodic)' if tmix is None else tmix}",
            font=dict(size=14, color=C_AXIS)
        ),
        plot_bgcolor=C_BG, paper_bgcolor=C_BG,
        height=520, width=900,
        margin=dict(l=40, r=20, t=80, b=40),
    )
    # eigenvalue plane: equal aspect ratio
    fig.update_xaxes(scaleanchor='y', scaleratio=1,
                     range=[-1.3, 1.3], row=1, col=1,
                     gridcolor=C_GRID, zerolinecolor=C_GRID)
    fig.update_yaxes(range=[-1.3, 1.3], row=1, col=1,
                     gridcolor=C_GRID, zerolinecolor=C_GRID)
    # TV plot
    fig.update_xaxes(title_text='Step k', row=2, col=1,
                     gridcolor=C_GRID)
    fig.update_yaxes(title_text='TV distance', range=[0, 0.55],
                     row=2, col=1, gridcolor=C_GRID)
    # stationary-distribution plot
    fig.update_xaxes(title_text='State', row=2, col=2,
                     gridcolor=C_GRID)
    fig.update_yaxes(title_text='π_i', row=2, col=2,
                     gridcolor=C_GRID)
    # heat maps: flip the y-axis (row 0 at the top)
    for col in range(2, 6):
        fig.update_yaxes(autorange='reversed', row=1, col=col)

    return fig, lam, gap, tmix

print_header("Unified comparison dashboard | functions defined")
print("  build_dashboard(P, chain_name, color)")
print("  Output: plotly FigureWidget + metric values")

# ============================================================
# §1 Interactive controls: dashboard for switching among the five chains
# ============================================================

# --- Step 1: define the parameter configuration of each chain ---
CHAIN_CONFIGS = {
    'Birth-Death': {
        'color': CHAIN_COLORS['birth'],
        'params': [
            widgets.IntSlider(value=8, min=4, max=14, step=1,
                              description='States S:', style={'description_width':'100px'}),
            widgets.FloatSlider(value=0.5, min=0.1, max=0.9, step=0.05,
                                description='P(right) p:', style={'description_width':'100px'},
                                readout_format='.2f'),
        ],
        'builder': lambda ps: make_birth_death(int(ps[0]), ps[1])
    },
    'Cycle': {
        'color': CHAIN_COLORS['cycle'],
        'params': [
            widgets.IntSlider(value=8, min=4, max=20, step=1,
                              description='Ring length N:', style={'description_width':'100px'}),
            widgets.FloatSlider(value=0.0, min=0.0, max=0.5, step=0.05,
                                description='Lazy prob. r:', style={'description_width':'100px'},
                                readout_format='.2f'),
        ],
        'builder': lambda ps: make_cycle(int(ps[0]), ps[1])
    },
    'Google Matrix': {
        'color': CHAIN_COLORS['google'],
        'params': [
            widgets.IntSlider(value=7, min=4, max=12, step=1,
                              description='Nodes N:', style={'description_width':'100px'}),
            widgets.FloatSlider(value=0.85, min=0.5, max=0.99, step=0.01,
                                description='Damping α:', style={'description_width':'100px'},
                                readout_format='.2f'),
        ],
        'builder': lambda ps: make_google(int(ps[0]), ps[1])
    },
    'Doubly Stochastic': {
        'color': CHAIN_COLORS['dbstoch'],
        'params': [
            widgets.IntSlider(value=6, min=3, max=12, step=1,
                              description='States N:', style={'description_width':'100px'}),
            widgets.FloatSlider(value=0.6, min=0.05, max=1.0, step=0.05,
                                description='Mixing m:', style={'description_width':'100px'},
                                readout_format='.2f'),
        ],
        'builder': lambda ps: make_doubly_stochastic(int(ps[0]), ps[1])
    },
    'Hypercube': {
        'color': CHAIN_COLORS['hypercube'],
        'params': [
            widgets.IntSlider(value=3, min=2, max=4, step=1,
                              description='Dimension n:', style={'description_width':'100px'}),
            widgets.FloatSlider(value=0.0, min=0.0, max=0.5, step=0.05,
                                description='Lazy prob. r:', style={'description_width':'100px'},
                                readout_format='.2f'),
        ],
        'builder': lambda ps: make_hypercube(int(ps[0]), ps[1])
    },
}

# --- Step 2: build the UI ---
chain_selector = widgets.ToggleButtons(
    options=list(CHAIN_CONFIGS.keys()),
    value='Birth-Death',
    description='Markov chain:',
    style={'description_width': '80px', 'button_width': '110px'},
)

output_fig   = widgets.Output()
param_box    = widgets.VBox([])
metrics_html = widgets.HTML(value='')

def update_params(_=None):
    cfg = CHAIN_CONFIGS[chain_selector.value]
    param_box.children = cfg['params']
    for sl in cfg['params']:
        sl.observe(update_figure, names='value')
    update_figure()

def update_figure(_=None):
    cfg  = CHAIN_CONFIGS[chain_selector.value]
    ps   = [sl.value for sl in cfg['params']]
    P    = cfg['builder'](ps)
    fig, lam, gap, tmix = build_dashboard(P, chain_selector.value, cfg['color'])
    metrics_html.value = (
        f"<div style='font-family:monospace; font-size:13px; "
        f"background:#F0EAF6; padding:8px 16px; border-radius:6px; margin:4px 0;'>"
        f"Matrix size: <b>{P.shape[0]}×{P.shape[0]}</b>&emsp;"
        f"|λ₁| = <b>{lam:.4f}</b>&emsp;"
        f"spectral gap = <b>{gap:.4f}</b>&emsp;"
        f"t_mix = <b>{'∞ (periodic; increase the lazy probability r)' if tmix is None else str(tmix)+' steps'}</b>"
        f"</div>"
    )
    with output_fig:
        output_fig.clear_output(wait=True)
        fig.show()

chain_selector.observe(update_params, names='value')
update_params()

display(widgets.VBox([
    chain_selector,
    param_box,
    metrics_html,
    output_fig,
]))

§2 The Birth-Death Chain¶

Physical Intuition: Hospital Beds¶

Imagine a hospital with $S+1$ beds numbered $0, 1, \ldots, S$ . Each day one of two things happens at random: a new patient is admitted (a “birth,” a step to the right) with probability $p$ , or a current patient is discharged (a “death,” a step to the left) with probability $q = 1-p$ . At the boundary (bed 0 or bed S) the probabilities are the same; only the direction is restricted.

This is the birth-death chain—at each step it can only move between adjacent states.

Matrix Structure¶

\mathbf{P} = \begin{pmatrix} q & p & 0 & \cdots & 0 \\ q & 0 & p & \cdots & 0 \\ 0 & q & 0 & \cdots & 0 \\ \vdots & & & \ddots & \vdots \\ 0 & 0 & 0 & q & p \end{pmatrix}

(10)

Tridiagonal matrix: the nonzero entries lie only on the main diagonal and on the one diagonal on each side of it
White cells (zeros) = states that cannot be reached in one step

Detailed Balance: Why Are All the Eigenvalues Real?¶

The birth-death chain satisfies the detailed balance condition:

\pi_i \cdot P_{i,i+1} = \pi_{i+1} \cdot P_{i+1,i}

(11)

Intuitively: the “flow” from state $i$ to $i+1$ equals the “flow” from $i+1$ back to $i$ . It is like balanced traffic on a two-way road—what goes in equals what comes out.

From detailed balance one can prove that $\mathbf{P}$ can be “symmetrized,” so its eigenvalues are all real, lying in $[-1, 1]$ .

Stationary distribution: $\pi_i = \pi_0 (p/q)^i$ , a geometric distribution (skewed to one side when $p \neq 0.5$ ).

A Concrete Example (First, What the Matrix Looks Like)¶

# ============================================================
print_header("§2.0 Birth-death chain | a concrete matrix (S=4, p=0.4)")
# ============================================================

# --- Step 1: construct and print the 5×5 matrix ---
print_step(1, "Construct make_birth_death(S=4, p=0.4)")
P_ex = make_birth_death(S=4, p=0.4)   # 5×5 matrix
print("\n  p = 0.4 (to the right), q = 0.6 (to the left)\n")
print("  P =\n")
print("          s0    s1    s2    s3    s4")
for i, row in enumerate(P_ex):
    row_str = '  '.join(f'{v:.2f}' for v in row)
    print(f"  s{i}      {row_str}")

# --- Step 2: check entry by entry against the formulas ---
print_step(2, "Entry-by-entry check against the formulas P[i,i-1]=q, P[i,i+1]=p")
checks = [
    ("P[0,0] = q (bounce at boundary)",  P_ex[0,0], 0.6),
    ("P[0,1] = p (advance at boundary)",  P_ex[0,1], 0.4),
    ("P[1,0] = q (to the left)",      P_ex[1,0], 0.6),
    ("P[1,2] = p (to the right)",      P_ex[1,2], 0.4),
    ("P[4,3] = q (left at boundary)",  P_ex[4,3], 0.6),
    ("P[4,4] = p (bounce at boundary)",  P_ex[4,4], 0.4),
    ("P[2,0] = 0 (no direct jump)",P_ex[2,0], 0.0),
]
all_ok = True
for label, val, expected in checks:
    ok = abs(val - expected) < 1e-10
    all_ok = all_ok and ok
    print(f"  {label:<34}  computed={val:.2f}  expected={expected:.2f}  {'✓' if ok else '✗'}")
print(f"\n  All checks: {'✓ passed' if all_ok else '✗ failed'}")

# --- Step 3: heat map (visualizing the structure) ---
print_step(3, "Heat map: tridiagonal structure (nonzero entries near the diagonal)")
colorscale_ex = [[0.0,"#FFFFFF"],[0.5,"#AB82C5"],[1.0,"#57068C"]]
fig_ex = go.Figure(go.Heatmap(
    z=P_ex,
    colorscale=colorscale_ex, zmin=0, zmax=1,
    text=[[f'{v:.2f}' for v in row] for row in P_ex],
    texttemplate='%{text}',
    textfont=dict(size=13),
    xgap=2, ygap=2,
    colorbar=dict(title='Probability', thickness=12)
))
fig_ex.update_layout(
    title='Transition matrix P of the birth-death chain (S=4, p=0.4): tridiagonal structure',
    xaxis=dict(title='Destination state j', tickvals=list(range(5)),
               ticktext=[f'state {j}' for j in range(5)],
               gridcolor=C_GRID),
    yaxis=dict(title='Starting state i', tickvals=list(range(5)),
               ticktext=[f'state {i}' for i in range(5)],
               autorange='reversed', gridcolor=C_GRID),
    plot_bgcolor=C_BG, paper_bgcolor=C_BG,
    height=380, width=460,
    margin=dict(l=80, r=60, t=60, b=60)
)
fig_ex.show()
print("  ▸ White cells = probability 0 (not reachable); dark violet = probability 0.6; light violet = probability 0.4")
print("  ▸ Nonzero entries lie only on the main diagonal and the diagonals just above and below it: tridiagonal structure")

# ============================================================
print_header("§2 Birth-death chain | checking detailed balance and real eigenvalues")
# ============================================================

# --- Step 1: construct the matrix ---
print_step(1, "Construct the birth-death chain P (S=6, p=0.6)")
P_bd = make_birth_death(S=6, p=0.6)
pi_bd = stationary_dist(P_bd)
print("P =\n", P_bd)

# --- Step 2: check the stationary distribution ---
print_step(2, "Stationary distribution π (geometric distribution, ratio p/q = 0.6/0.4 = 1.5)")
compare_print(
    "π[i] / π[i-1]",
    np.round(pi_bd[1:] / pi_bd[:-1], 4),
    "should all equal p/q = 1.5"
)

# --- Step 3: check detailed balance ---
print_step(3, "Detailed balance: π_i · P_{i,i+1} = π_{i+1} · P_{i+1,i}")
lhs = pi_bd[:-1] * P_bd[:-1, 1:]  # π_i * p
rhs = pi_bd[1:] * P_bd[1:, :-1]   # π_{i+1} * q (only the diagonal entries are used)
compare_print(
    "π_i·p - π_{i+1}·q",
    np.round(np.diag(lhs) - np.diag(rhs), 10),
    "should all be 0 (detailed balance holds)"
)

# --- Step 4: eigenvalues (all real) ---
print_step(4, "Eigenvalues: should all be real, lying in [-1, 1]")
eigs_bd = np.sort(np.linalg.eigvals(P_bd).real)[::-1]
compare_print(
    "Eigenvalues (largest to smallest)",
    np.round(eigs_bd, 4),
    "λ_0=1, all other |λ_k| < 1, all real (a reversible chain can be symmetrized)"
)
lam_bd = lambda1_abs(P_bd)
compare_print(
    "Mixing-time estimate",
    mixing_time(lam_bd),
    f"t_mix ≈ log(2) / gap = {np.log(2)/(1-lam_bd):.1f}"
)

§3 The Cycle Chain¶

Physical Intuition: A Circular Track¶

$N$ runners stand on a circular track; at each step each one independently moves forward one position with probability $1/2$ and back one position with probability $1/2$ (the lazy probability $r$ is the probability of staying in place).

The ring structure brings an important property: the “right neighbor” of the last state $N-1$ is state 0.

Why Does the DFT Diagonalize It?¶

The transition matrix is a circulant matrix: each row is the previous row cyclically shifted one position to the right. Circulant matrices have a perfect property—the column vectors of the DFT (discrete Fourier transform) matrix are exactly their eigenvectors, and the eigenvalues can be computed exactly:

\lambda_k = r + (1-r)\cos\!\left(\frac{2\pi k}{N}\right), \quad k = 0, 1, \ldots, N-1

(12)

This means that the eigenvalues are all real (the transition matrix is symmetric) and lie in the interval $[2r-1, 1]$ of the real axis; only $\lambda_0 = 1$ (and, when $r = 0$ and $N$ is even, $\lambda_{N/2} = -1$ ) lies on the unit circle.

The Periodicity Problem of Even-Length Rings¶

When $N$ is even and $r=0$ , $\lambda_{N/2} = \cos(\pi) = -1$ .

What does this mean? $(-1)^k$ oscillates back and forth between odd and even steps and never decays! As a result, the TV distance does not decrease monotonically; instead it jumps up every other step.

The remedy: add a lazy probability $r > 0$ , so that $|\lambda_{N/2}| < 1$ ; the periodicity disappears and the TV distance becomes monotonically decreasing.

Connection with §9.2.2

For $N=2$ , the DFT matrix $\frac{1}{\sqrt{2}}\begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}$ is the Hadamard gate of quantum mechanics (its columns form the Hadamard basis; for general $N$ , the DFT matrix corresponds to the quantum Fourier transform). In the special case $N=2$ , the transition matrix of the cycle chain is $\begin{pmatrix} r & 1-r \\ 1-r & r \end{pmatrix}$ (stepping forward and stepping backward both lead to the other state; the code make_cycle writes these two directions into the same entry, so for $N=2$ its off-diagonal entries are only $(1-r)/2$ , and the function applies only for $N \geq 3$ ), which is exactly an analog of the quantum Hadamard gate ( $r=0.5$ ).

A Concrete Example (First, What the Matrix Looks Like)¶

# ============================================================
print_header("§3.0 Cycle chain | a concrete matrix (N=6, r=0)")
# ============================================================

# --- Step 1: construct and print the 6×6 matrix ---
print_step(1, "Construct make_cycle(N=6, r=0)")
P_cy = make_cycle(N=6, r=0)
print("\n  At each step move left or right with probability 1/2 (lazy probability r=0)\n")
print("          s0    s1    s2    s3    s4    s5")
for i, row in enumerate(P_cy):
    row_str = '  '.join(f'{v:.2f}' for v in row)
    print(f"  s{i}      {row_str}")

# --- Step 2: check the circulant structure ---
print_step(2, "Circulant structure: each row is a cyclic shift of the previous row")
checks = [
    ("P[0,1] = P[0,5] = 1/2 (left-right symmetric)", P_cy[0,1], P_cy[0,5], 0.5),
    ("P[5,0] = 1/2 (the ring closes up)",           P_cy[5,0], P_cy[5,0], 0.5),
    ("P[3,2] = P[3,4] = 1/2",             P_cy[3,2], P_cy[3,4], 0.5),
    ("P[0,0] = 0 (never stays in place)",           P_cy[0,0], P_cy[0,0], 0.0),
]
for label, v1, v2, expected in checks:
    ok = abs(v1-expected)<1e-10 and abs(v2-expected)<1e-10
    print(f"  {label:<46}  ={'✓' if ok else '✗'}")

# --- Step 3: comparison with the lazy probability r=0.2 ---
print_step(3, "Comparison: adding the lazy probability r=0.2")
P_cy_lazy = make_cycle(N=6, r=0.2)
print("\n  P (r=0.2) =\n")
print("          s0    s1    s2    s3    s4    s5")
for i, row in enumerate(P_cy_lazy):
    row_str = '  '.join(f'{v:.2f}' for v in row)
    print(f"  s{i}      {row_str}")
print("  ▸ The diagonal changes from 0 to 0.20 (probability of staying in place)")
print("  ▸ The off-diagonal entries shrink from 0.50 to 0.40")

# --- Step 4: heat-map comparison (r=0 vs r=0.2) ---
print_step(4, "Heat-map comparison: r=0 (left) vs r=0.2 (right)")
colorscale_ex = [[0.0,"#FFFFFF"],[0.5,"#AB82C5"],[1.0,"#57068C"]]
fig_cy = make_subplots(rows=1, cols=2,
    subplot_titles=('r=0 (not lazy, diagonal is 0)', 'r=0.2 (lazy, diagonal nonzero)'))
for col_idx, (P_show, title) in enumerate([(P_cy, 'r=0'), (P_cy_lazy, 'r=0.2')], 1):
    fig_cy.add_trace(go.Heatmap(
        z=P_show,
        colorscale=colorscale_ex, zmin=0, zmax=0.6,
        text=[[f'{v:.2f}' for v in row] for row in P_show],
        texttemplate='%{text}', textfont=dict(size=11),
        xgap=2, ygap=2, showscale=(col_idx==2)
    ), row=1, col=col_idx)
    fig_cy.update_yaxes(autorange='reversed', row=1, col=col_idx)
fig_cy.update_layout(
    title='Transition matrix P of the cycle chain (N=6)',
    plot_bgcolor=C_BG, paper_bgcolor=C_BG,
    height=320, width=680,
    margin=dict(l=40, r=40, t=60, b=40)
)
fig_cy.show()

# ============================================================
print_header("§3 Cycle chain | DFT eigenvalues and periodicity")
# ============================================================

N_c = 8

# --- Step 1: check the exact DFT formula for the eigenvalues ---
print_step(1, f"Cycle chain C_{N_c}, r=0: checking the DFT eigenvalue formula")
P_cyc = make_cycle(N=N_c, r=0.0)
eigs_cyc_num = np.sort(np.linalg.eigvals(P_cyc).real)[::-1]
k_vals = np.arange(N_c)
eigs_cyc_dft = np.cos(2 * np.pi * k_vals / N_c)  # analytic formula for r=0
eigs_cyc_dft_sorted = np.sort(eigs_cyc_dft)[::-1]
compare_print(
    "Numerical eigenvalues",
    np.round(eigs_cyc_num, 4),
    f"DFT formula cos(2πk/{N_c}): {np.round(eigs_cyc_dft_sorted, 4)}"
)

# --- Step 2: the problem with even-length rings ---
print_step(2, f"C_{N_c} (even): λ_{{N/2}} = {eigs_cyc_num[-1]:.4f}, which makes the TV distance oscillate")
pi_cyc = stationary_dist(P_cyc)
mu = np.zeros(N_c); mu[0] = 1.0
tv_vals = []
for _ in range(30):
    tv_vals.append(tv_dist(mu, pi_cyc))
    mu = mu @ P_cyc
print("  TV distance, first 10 steps:", np.round(tv_vals[:10], 4))
print("  TV distance between odd and even steps =>", "oscillating" if tv_vals[1] > tv_vals[0] else "monotone")

# --- Step 3: add the lazy probability r=0.1 ---
print_step(3, "Add the lazy probability r=0.1 to break the even-length periodicity")
P_cyc_lazy = make_cycle(N=N_c, r=0.1)
eigs_lazy = np.sort(np.linalg.eigvals(P_cyc_lazy).real)[::-1]
compare_print(
    f"With laziness added, λ_{{N/2}}",
    f"{eigs_lazy[-1]:.4f}",
    f"r + (1-r)·cos(π) = 0.1 + 0.9·(-1) = -0.8, |λ| < 1"
)
lam_lazy = lambda1_abs(P_cyc_lazy)
compare_print(
    "|λ₁| (lazy chain)",
    f"{lam_lazy:.4f}",
    f"r + (1-r)·cos(2π/{N_c}) = 0.1 + 0.9·{np.cos(2*np.pi/N_c):.4f} = {0.1 + 0.9*np.cos(2*np.pi/N_c):.4f}"
)

§4 Google Matrix (PageRank)¶

Motivation: How Do You Judge How Important a Web Page Is?¶

In the 1990s, search engines faced a problem: of the billions of pages on the web, which ones are “truly important”? Simply counting how many times a page is linked to (its in-degree) is not enough— because spam sites can link to one another to inflate their counts.

The idea of Larry Page and Sergey Brin (the founders of Google) was: let a random surfer judge for you.

Simulate a person who clicks links at random as they move among web pages. If they end up spending most of their time on some page, that page is “important”—this is PageRank.

The behavior of this random surfer is exactly the stationary distribution of a Markov chain!

The Underlying Structure: An Irregular Directed Graph¶

The real web is an irregular directed graph: some pages (such as Wikipedia) are linked to by thousands of pages (high in-degree), while others are linked to only once or twice. This kind of distribution, in which “a few nodes are extremely important,” is called a power-law distribution, and it is a pervasive feature of real networks.

We generate such irregular graphs with the Barabási–Albert (BA) preferential attachment model: each newly added node chooses its link targets with probability proportional to the in-degrees of the existing nodes— “the rich get richer,” and nodes with high in-degree are more likely to receive new links.

Two Technical Problems and Google’s Solutions¶

Problem 1: dangling nodes

Some web pages have no outgoing links at all (such as pure image pages); a surfer who arrives there is “stuck” and cannot continue.

Solution: make the outgoing links of these nodes uniform—equivalent to “jumping at random to any page.”

Problem 2: disconnected components

The web may consist of several subgraphs that are not connected to one another, so that the stationary distribution is not unique (or does not exist at all).

Solution: add a teleportation probability $1-\alpha$ —with probability $1-\alpha$ the surfer jumps at random to any page, regardless of where they currently are. This is the Google correction:

\mathbf{G} = \alpha \mathbf{P}_{\text{link}} + (1-\alpha) \cdot \frac{1}{N} \mathbf{1}\mathbf{1}^\top

(13)

The Mathematical Effect: A Rank-One Correction¶

$(1-\alpha)\cdot\frac{1}{N}\mathbf{1}\mathbf{1}^\top$ is a rank-one matrix: all its entries are the same, equal to $(1-\alpha)/N$ .

After this matrix is added, all the entries of $\mathbf{G}$ are strictly greater than 0, which guarantees irreducibility (any two web pages can reach each other with positive probability), and hence a unique stationary distribution.

More importantly, this correction “compresses” all the non-dominant eigenvalues to $|\lambda_k| \leq \alpha$ :

\text{spectral gap} = 1 - |\lambda_1| \geq 1 - \alpha

(14)

$\alpha = 0.85$ (Google’s actual setting) means that the spectral gap is at least 0.15, so convergence takes roughly $\lceil \log 2 / 0.15 \rceil \approx 5$ steps.

The Trade-off between Speed and Quality¶

The larger $\alpha$ is (close to 1): the more “faithful” PageRank is (the more it respects the link structure), but the slower the convergence
The smaller $\alpha$ is (close to 0): convergence is fast, but PageRank tends to the uniform distribution and loses its meaning

A Concrete Example (First, What the Matrix Looks Like)¶

# ============================================================
print_header("§4.0 Google Matrix | concrete examples with a BA graph and an SBM")
# ============================================================
# This cell shows two kinds of irregular directed graphs as P_link,
# illustrating how PageRank turns "link structure" into "node importance"

import numpy as np

# --- Step 1: Barabási–Albert (BA) preferential-attachment graph ---
print_step(1, "BA graph: preferential attachment produces a power-law degree distribution")

def make_ba_stochastic(N=8, m=2, seed=42):
    """
    Barabási–Albert directed graph → stochastic matrix
    N : final number of nodes
    m : each new node links to m existing nodes (preferential attachment)
    Returns: an N×N row-stochastic matrix (each row sums to 1)
    """
    rng = np.random.default_rng(seed)
    # initialization: m+1 nodes, every pair linked to each other
    adj = np.zeros((N, N), dtype=float)
    for i in range(m + 1):
        for j in range(m + 1):
            if i != j:
                adj[i, j] = 1.0
    in_degree = adj.sum(axis=0).copy()   # in-degree
    # add the new nodes one at a time
    for new in range(m + 1, N):
        # preferential attachment: choose targets with probability proportional to in-degree
        probs = in_degree[:new].copy()
        if probs.sum() == 0:
            probs = np.ones(new)
        probs = probs / probs.sum()
        targets = rng.choice(new, size=m, replace=False, p=probs)
        for t in targets:
            adj[new, t] = 1.0      # new → t
            in_degree[t] += 1
    # handle nodes with out-degree 0 (dangling nodes get uniform outgoing jumps)
    for i in range(N):
        if adj[i].sum() == 0:
            adj[i] = np.ones(N) / N
        else:
            adj[i] /= adj[i].sum()
    return adj

N_ba = 8
P_ba_link = make_ba_stochastic(N=N_ba, m=2, seed=7)
print(f"\n  BA graph, N={N_ba} nodes, m=2 (each new node links to 2 existing nodes)")
print("  P_link (row-stochastic matrix; nonzero means there is a link):\n")
print("       " + '  '.join(f'node{j}' for j in range(N_ba)))
for i, row in enumerate(P_ba_link):
    row_str = '  '.join(f'{v:.3f}' if v > 0.001 else '  —  ' for v in row)
    print(f"  node{i}  {row_str}")

# in-degree statistics
in_deg = (P_ba_link > 0.001).sum(axis=0)
print(f"\n  Number of links into each node (in-degree):")
for j in range(N_ba):
    bar = '█' * int(in_deg[j])
    print(f"    node{j}: {in_deg[j]:2d}  {bar}")
print("  ▸ The in-degree distribution is uneven—this is the embryonic form of a power-law distribution")

# --- Step 2: add Google's α correction ---
print_step(2, "Apply the Google correction: G = α·P_link + (1-α)·(1/N)·11^T")
alpha_ba = 0.85
G_ba = alpha_ba * P_ba_link + (1 - alpha_ba) / N_ba * np.ones((N_ba, N_ba))
pi_ba = stationary_dist(G_ba)
idx_sorted = np.argsort(pi_ba)[::-1]   # PageRank from highest to lowest

print(f"\n  α = {alpha_ba}, PageRank (stationary distribution π) ranking:")
print(f"  {'Rank':>4}  {'Node':>6}  {'PageRank':>10}  {'In-deg':>6}  {'Bar'}")
print("  " + "-" * 52)
for rank, j in enumerate(idx_sorted):
    bar = '▓' * int(pi_ba[j] * 80)
    print(f"  {rank+1:>4}  node{j:<3}  {pi_ba[j]:>10.4f}  {in_deg[j]:>6}  {bar}")
compare_print(
    "Sum of the PageRanks",
    f"{pi_ba.sum():.6f}",
    "should be 1.0 (normalized probability distribution)"
)

# --- Step 3: Stochastic Block Model (SBM) ---
print_step(3, "SBM: a graph with community structure (2 communities, dense within and sparse between)")

def make_sbm_stochastic(sizes, p_in=0.7, p_out=0.05, seed=42):
    """
    Stochastic Block Model → row-stochastic matrix
    sizes  : list of the numbers of nodes in the communities, e.g. [4, 4]
    p_in   : link probability within a community
    p_out  : link probability between communities
    """
    rng = np.random.default_rng(seed)
    N = sum(sizes)
    adj = np.zeros((N, N), dtype=float)
    # build the community labels
    labels = []
    for k, s in enumerate(sizes):
        labels.extend([k] * s)
    labels = np.array(labels)
    # generate the directed edges at random
    for i in range(N):
        for j in range(N):
            if i == j:
                continue
            p = p_in if labels[i] == labels[j] else p_out
            if rng.random() < p:
                adj[i, j] = 1.0
    # fill in the dangling nodes and normalize
    for i in range(N):
        if adj[i].sum() == 0:
            adj[i] = np.ones(N) / N
        else:
            adj[i] /= adj[i].sum()
    return adj, labels

P_sbm_link, sbm_labels = make_sbm_stochastic(
    sizes=[4, 4], p_in=0.75, p_out=0.06, seed=13
)
N_sbm = P_sbm_link.shape[0]
print(f"\n  SBM, N={N_sbm} nodes, 2 communities of 4 members each")
print("  p_in=0.75 (within communities), p_out=0.06 (between communities)")
print("\n  P_link (community A = nodes 0-3, community B = nodes 4-7):\n")
print("       " + '  '.join(
    f'{"A" if j<4 else "B"}{j%4}' for j in range(N_sbm)
))
for i, row in enumerate(P_sbm_link):
    grp = 'A' if i < 4 else 'B'
    row_str = '  '.join(f'{v:.2f}' if v > 0.001 else '  — ' for v in row)
    print(f"  {grp}{i%4}     {row_str}")

# --- Step 4: PageRank of the SBM ---
print_step(4, "PageRank of the SBM: importance reflects community size and link density")
alpha_sbm = 0.85
G_sbm = alpha_sbm * P_sbm_link + (1 - alpha_sbm) / N_sbm * np.ones((N_sbm, N_sbm))
pi_sbm = stationary_dist(G_sbm)
pi_A = pi_sbm[:4].sum()
pi_B = pi_sbm[4:].sum()
print(f"\n  Total PageRank of community A (nodes 0-3): {pi_A:.4f}")
print(f"  Total PageRank of community B (nodes 4-7): {pi_B:.4f}")
print(f"  ▸ The difference between the communities reflects the difference in their link densities")

# --- Step 5: heat-map comparison (BA vs SBM, P_link vs G) ---
print_step(5, "Heat-map comparison: P_link (underlying structure) vs G (after the Google correction)")
colorscale_ex = [[0.0,"#FFFFFF"],[0.4,"#AB82C5"],[1.0,"#57068C"]]

fig_go2 = make_subplots(
    rows=2, cols=2,
    subplot_titles=(
        'BA graph: P_link (underlying, power-law in-degree)',
        'BA graph: G = α·P_link + (1-α)·U (Google correction)',
        'SBM: P_link (underlying, community structure)',
        'SBM: G (after the Google correction; communities still visible)',
    ),
    vertical_spacing=0.12,
    horizontal_spacing=0.08,
)

def add_heatmap(fig, M, row, col, show_text=True, show_scale=False, zmax=1.0):
    n = M.shape[0]
    text = [[f'{v:.2f}' if v > 0.005 else '' for v in r] for r in M]
    fig.add_trace(go.Heatmap(
        z=M,
        colorscale=colorscale_ex, zmin=0, zmax=zmax,
        text=text if show_text else None,
        texttemplate='%{text}' if show_text else None,
        textfont=dict(size=9),
        xgap=1, ygap=1, showscale=show_scale,
    ), row=row, col=col)
    fig.update_yaxes(autorange='reversed', row=row, col=col)

add_heatmap(fig_go2, P_ba_link, 1, 1, show_text=True,  show_scale=False)
add_heatmap(fig_go2, G_ba,      1, 2, show_text=True,  show_scale=True)
add_heatmap(fig_go2, P_sbm_link,2, 1, show_text=True,  show_scale=False)
add_heatmap(fig_go2, G_sbm,     2, 2, show_text=True,  show_scale=True)

# add community dividing lines to the SBM heat maps
for col in [1, 2]:
    fig_go2.add_shape(
        type='line', x0=3.5, x1=3.5, y0=-0.5, y1=7.5,
        line=dict(color=C_WARN, width=1.5, dash='dot'),
        row=2, col=col
    )
    fig_go2.add_shape(
        type='line', x0=-0.5, x1=7.5, y0=3.5, y1=3.5,
        line=dict(color=C_WARN, width=1.5, dash='dot'),
        row=2, col=col
    )

fig_go2.update_layout(
    title=dict(
        text='Google Matrix: underlying structure vs after the α correction (α=0.85)',
        font=dict(size=13)
    ),
    plot_bgcolor=C_BG, paper_bgcolor=C_BG,
    height=560, width=720,
    margin=dict(l=50, r=50, t=80, b=40),
)
fig_go2.show()

print("\n  ▸ BA graph: P_link is very sparse (many white cells); in G every cell has a small positive value")
print("  ▸ SBM: dark within communities (high probability), light between communities (low probability)")
print("  ▸ The orange dotted lines mark the boundary between the two SBM communities")
print("  ▸ The essence of the Google correction: turning a sparse matrix into an all-positive matrix guarantees a unique stationary distribution")

# ============================================================
print_header("§4 Google Matrix | the rank-one correction and the lower bound on the spectral gap")
# ============================================================

N_g = 8

# --- Step 1: eigenvalues of the original chain vs the Google Matrix ---
print_step(1, f"Compare the original chain with the Google Matrix (N={N_g}, α=0.85)")
P_link = make_cycle(N=N_g, r=0.0)  # use the cycle chain as P_link
P_google = make_google(N=N_g, alpha=0.85)

eigs_link   = np.sort(np.abs(np.linalg.eigvals(P_link)))[::-1]
eigs_google = np.sort(np.abs(np.linalg.eigvals(P_google)))[::-1]
compare_print(
    "|λ_k| original cycle chain",
    np.round(eigs_link, 4),
    "includes the real eigenvalue λ=−1, with |λ|=1 (periodicity)"
)
compare_print(
    "|λ_k| Google Matrix (α=0.85)",
    np.round(eigs_google, 4),
    f"non-dominant eigenvalues ≤ α = 0.85 (compressed by the rank-one correction)"
)

# --- Step 2: sweep α and observe the spectral gap ---
print_step(2, "Sweep α ∈ [0.5, 0.99] and verify gap ≥ 1-α")
alphas = [0.5, 0.7, 0.85, 0.95, 0.99]
print(f"  {'α':>6}  {'|λ₁|':>8}  {'gap':>8}  {'1-α':>8}  {'gap ≥ 1-α?':>12}")
print("  " + "-" * 52)
for a in alphas:
    Pg = make_google(N=N_g, alpha=a)
    pi_g = stationary_dist(Pg)
    l2 = lambda1_abs(Pg)
    gap = 1 - l2
    check = "✓" if gap >= (1 - a) - 1e-6 else "✗"
    print(f"  {a:>6.2f}  {l2:>8.4f}  {gap:>8.4f}  {1-a:>8.4f}  {check:>12}")

§5 Doubly Stochastic Matrices¶

Physical Intuition: A Balanced Logistics Network¶

Suppose there are $N$ cities, with goods flowing among them. An ordinary stochastic matrix only requires that “the total amount of goods leaving each city = 1” (row sums equal 1). A doubly stochastic matrix also requires that “the total amount of goods flowing into each city = 1” (column sums also equal 1)— like a perfectly balanced logistics network in which every city ships out as much as it takes in.

Definition¶

\sum_{j=0}^{N-1} P_{ij} = 1 \quad (\text{each row sums to 1: balanced out-degree})

(15)

\sum_{i=0}^{N-1} P_{ij} = 1 \quad (\text{each column sums to 1: balanced in-degree})

(16)

Note: an ordinary stochastic matrix satisfies only the first condition; a doubly stochastic matrix must satisfy both.

Why Must the Stationary Distribution Be Uniform?¶

An intuitive derivation: let $\boldsymbol{\pi} = (1/N, 1/N, \ldots, 1/N)$ (the uniform distribution). Check the $j$ -th entry of $\boldsymbol{\pi}^\top \mathbf{P}$ :

[\boldsymbol{\pi}^\top \mathbf{P}]_j = \sum_{i=0}^{N-1} \pi_i P_{ij} = \frac{1}{N} \sum_{i=0}^{N-1} P_{ij} = \frac{1}{N} \cdot 1 = \pi_j

(17)

It holds! The reason is that column $j$ sums to 1 (the second condition for a doubly stochastic matrix). Hence, whatever the particular shape of the matrix, the stationary distribution is always uniform.

The Birkhoff–von Neumann Theorem¶

This theorem says: the doubly stochastic matrices are precisely the convex combinations of permutation matrices.

A permutation matrix is a matrix obtained by rearranging the rows of $\mathbf{I}$ , for example:

\begin{pmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix}

(18)

The intuitive meaning of this theorem is that any balanced logistics plan can be decomposed into a weighted average of several “pure permutation plans.”

A Concrete Example (First, What the Matrix Looks Like)¶

# ============================================================
print_header("§5.0 Doubly stochastic matrix | a concrete matrix (N=4, m=0.7)")
# ============================================================

# --- Step 1: construct and print the 4×4 matrix ---
print_step(1, "Construct make_doubly_stochastic(N=4, m=0.7)")
P_ds4 = make_doubly_stochastic(N=4, m=0.7)
print("\n  m = 0.7 (mixing strength)\n")
print("          s0      s1      s2      s3")
for i, row in enumerate(P_ds4):
    row_str = '  '.join(f'{v:.4f}' for v in row)
    print(f"  s{i}      {row_str}")

# --- Step 2: check double stochasticity ---
print_step(2, "Check: each row sums to 1 and each column sums to 1")
row_sums = P_ds4.sum(axis=1)
col_sums = P_ds4.sum(axis=0)
print("  Row sums:   ", np.round(row_sums, 6))
print("  Column sums:", np.round(col_sums, 6))
ok_row = np.allclose(row_sums, 1)
ok_col = np.allclose(col_sums, 1)
print(f"  All row sums are 1: {'✓' if ok_row else '✗'}  All column sums are 1: {'✓' if ok_col else '✗'}")
print(f"  ▸ An ordinary stochastic matrix guarantees only row sums = 1; a doubly stochastic matrix also requires column sums = 1")

# --- Step 3: comparison with an ordinary stochastic matrix ---
print_step(3, "Comparison: the column sums of an ordinary stochastic matrix (birth-death chain) are not 1")
P_bd4 = make_birth_death(S=3, p=0.4)  # 4×4 birth-death chain
print("  Column sums of the birth-death chain (S=3, p=0.4):",
      np.round(P_bd4.sum(axis=0), 4))
print("  ▸ The column sums are not 1 (the birth-death chain is not doubly stochastic)")

# --- Step 4: heat maps (m=0.1 vs m=0.7 vs m=1.0) ---
print_step(4, "Heat maps: how the matrix looks for different mixing strengths m")
colorscale_ex = [[0.0,"#FFFFFF"],[0.5,"#AB82C5"],[1.0,"#57068C"]]
fig_ds = make_subplots(rows=1, cols=3,
    subplot_titles=('m=0.1 (nearly uniform)', 'm=0.7', 'm=1.0 (close to a permutation matrix)'))
for col_idx, m in enumerate([0.1, 0.7, 1.0], 1):
    P_show = make_doubly_stochastic(N=4, m=m)
    fig_ds.add_trace(go.Heatmap(
        z=P_show,
        colorscale=colorscale_ex, zmin=0, zmax=1,
        text=[[f'{v:.3f}' for v in row] for row in P_show],
        texttemplate='%{text}', textfont=dict(size=11),
        xgap=2, ygap=2, showscale=(col_idx==3)
    ), row=1, col=col_idx)
    fig_ds.update_yaxes(autorange='reversed', row=1, col=col_idx)
fig_ds.update_layout(
    title='Doubly stochastic matrix (N=4): different mixing strengths m',
    plot_bgcolor=C_BG, paper_bgcolor=C_BG,
    height=300, width=720,
    margin=dict(l=40, r=40, t=60, b=40)
)
fig_ds.show()
print("  ▸ As m→1.0 the matrix approaches a permutation matrix (the extreme case of the Birkhoff theorem)")
print("  ▸ However m changes, the stationary distribution is always uniform")

# ============================================================
print_header("§5 Doubly stochastic matrix | uniform stationary distribution and the Birkhoff theorem")
# ============================================================

N_ds = 6

# --- Step 1: check the row and column sums ---
print_step(1, f"Row and column sums of the doubly stochastic matrix P (N={N_ds}, m=0.6)")
P_ds = make_doubly_stochastic(N=N_ds, m=0.6)
compare_print(
    "Row sums",
    np.round(P_ds.sum(axis=1), 6),
    "should all be 1.0"
)
compare_print(
    "Column sums",
    np.round(P_ds.sum(axis=0), 6),
    "should all be 1.0 (definition of a doubly stochastic matrix)"
)

# --- Step 2: check the stationary distribution ---
print_step(2, "Stationary distribution: should be the uniform distribution π = (1/N)1")
pi_ds = stationary_dist(P_ds)
compare_print(
    "Numerical stationary π",
    np.round(pi_ds, 6),
    f"uniform distribution (1/{N_ds})·1 = {1/N_ds:.6f}"
)

# --- Step 3: changing m does not affect the stationary distribution ---
print_step(3, "Change the mixing strength m: the stationary distribution stays uniform")
for m in [0.1, 0.3, 0.6, 0.9, 1.0]:
    P_tmp = make_doubly_stochastic(N=N_ds, m=m)
    pi_tmp = stationary_dist(P_tmp)
    max_dev = np.max(np.abs(pi_tmp - 1/N_ds))
    lam_tmp = lambda1_abs(P_tmp)
    print(f"  m={m:.1f}: max|π_i - 1/{N_ds}| = {max_dev:.2e}, |λ₁| = {lam_tmp:.4f}")

§6 The Hypercube Walk¶

Physical Intuition: Randomly Flipping the Bits of a Password¶

Suppose you have an $n$ -bit binary password; for example, when $n=3$ the possible passwords are $000, 001, 010, \ldots, 111$ , $2^n = 8$ of them in all.

At each step, one bit is chosen at random and flipped ( $0 \to 1$ or $1 \to 0$ ). For example, starting from 011, flipping bit 0 gives 010, flipping bit 1 gives 001, and flipping bit 2 gives 111—each of the three options has probability $1/3$ .

This is the hypercube walk: there are $2^n$ vertices, and two vertices are joined by an edge only when their Hamming distance is 1 (they differ in exactly one bit).

Matrix Structure: The Kronecker Product¶

The transition matrix can be expressed with the Kronecker product structure of Chapter 5:

\mathbf{P}_n = \frac{1}{n} \sum_{k=0}^{n-1} \underbrace{\mathbf{I}_{2^k} \otimes \mathbf{S} \otimes \mathbf{I}_{2^{n-k-1}}}_{\text{flip bit } k}, \quad \mathbf{S} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}

(19)

An intuitive reading: $\mathbf{S}$ is the $2 \times 2$ matrix that “flips one bit,” $\otimes$ “embeds” it at the position of bit $k$ , and $\mathbf{I}$ leaves the other bits unchanged.

Binomial Multiplicities of the Eigenvalues¶

This Kronecker product structure allows the eigenvalues to be computed exactly:

\lambda_j = 1 - \frac{2j}{n}, \quad j = 0, 1, \ldots, n

(20)

with multiplicity $\binom{n}{j}$ (a binomial coefficient!). For example, for $n=3$ :

$j$	$\lambda_j = 1 - 2j/3$	Multiplicity $\binom{3}{j}$
0	1	1
1	$1/3 \approx 0.333$	3
2	$-1/3 \approx -0.333$	3
3	-1	1

$\lambda_n = -1$ means that the hypercube walk is periodic (for every value of $n$ )— adding a lazy probability $r > 0$ removes the periodicity.

The Cutoff Phenomenon¶

The hypercube walk has a peculiar property: mixing does not happen by “gradual convergence” but by “sudden convergence.”

For a long time the TV distance stays close to 1 (mixing has barely begun); then, around $t \approx \frac{n}{2}\ln n$ , it drops sharply to nearly 0 (for the standard version with lazy probability $r = 1/2$ ).

This “cliff-like” convergence is called the cutoff phenomenon, and the mixing time is $t_{\text{mix}} \sim \frac{n}{2}\ln n$ . Each time the dimension increases by 1, the mixing time increases by roughly $\frac{1}{2}(\ln n + 1)$ steps.

A Concrete Example (First, What the Matrix Looks Like)¶

# ============================================================
print_header("§6.0 Hypercube | concrete matrices (n=2 and n=3)")
# ============================================================

# --- Step 1: the 4×4 matrix for n=2 (fully readable) ---
print_step(1, "n=2: the four vertices of {0,1}², a 4×4 matrix")
P_hc2 = make_hypercube(n=2, r=0.0)
print("\n  States: 00=0, 01=1, 10=2, 11=3")
print("  At each step one bit, chosen uniformly, is flipped (2 bits in all)\n")
print("          s00   s01   s10   s11")
state_labels = ['00','01','10','11']
for i, row in enumerate(P_hc2):
    row_str = '  '.join(f'{v:.2f}' for v in row)
    print(f"  s{state_labels[i]}     {row_str}")
print()

# --- Step 2: entry-by-entry check (which states can be reached from which in one step) ---
print_step(2, "Check: P[i,j] = 1/2 if and only if i⊕j is a power of 2")
checks_hc = [
    ("P[00,01]=1/2 (flip bit 0, 0→1)",    P_hc2[0,1], 0.5),
    ("P[00,10]=1/2 (flip bit 1, 0→1)",    P_hc2[0,2], 0.5),
    ("P[00,11]=0   (needs 2 flips, unreachable)", P_hc2[0,3], 0.0),
    ("P[01,00]=1/2 (flip bit 0, 1→0)",    P_hc2[1,0], 0.5),
    ("P[11,01]=1/2 (flip bit 1, 1→0)",    P_hc2[3,1], 0.5),
    ("P[00,00]=0   (never stays in place)",           P_hc2[0,0], 0.0),
]
all_ok = True
for label, val, expected in checks_hc:
    ok = abs(val - expected) < 1e-10
    all_ok = all_ok and ok
    print(f"  {label:<42}  computed={val:.2f}  expected={expected:.2f}  {'✓' if ok else '✗'}")
print(f"\n  All checks: {'✓ passed' if all_ok else '✗ failed'}")

# --- Step 3: the 8×8 matrix for n=3 (printed with labels) ---
print_step(3, "n=3: the eight vertices of {0,1}³, an 8×8 matrix")
P_hc3 = make_hypercube(n=3, r=0.0)
print("\n  States: 000=0, 001=1, 010=2, ..., 111=7")
print("  At each step one bit, chosen uniformly, is flipped (3 bits in all); P[i,j] = 1/3 or 0\n")
labels3 = [f'{i:03b}' for i in range(8)]
print("       " + '  '.join(f'{l:>5}' for l in labels3))
for i, row in enumerate(P_hc3):
    row_str = ''.join(f'  {v:.2f}' if v > 0 else '  0.00' for v in row)
    print(f"  {labels3[i]}  {row_str}")

# --- Step 4: showing the eigenvalue multiplicities ---
print_step(4, "Eigenvalue multiplicities: λ_j = 1-2j/n, multiplicity = C(n,j)")
from math import comb
print(f"  {'j':>3}  {'λ_j = 1-2j/3':>16}  {'C(3,j)':>8}  {'actual':>8}")
print("  " + "-" * 42)
eigs_hc3 = np.linalg.eigvals(P_hc3).real
for j in range(4):
    lam_j = 1 - 2*j/3
    theory = comb(3, j)
    actual = int(np.sum(np.abs(eigs_hc3 - lam_j) < 0.01))
    print(f"  {j:>3}  {lam_j:>16.4f}  {theory:>8}  {actual:>8}")

# --- Step 5: heat maps (n=2 vs n=3) ---
print_step(5, "Heat maps: n=2 (4×4) vs n=3 (8×8)")
colorscale_ex = [[0.0,"#FFFFFF"],[0.5,"#AB82C5"],[1.0,"#57068C"]]
fig_hc = make_subplots(rows=1, cols=2,
    subplot_titles=(
        'n=2 (4×4): exactly two 1/2s in each row',
        'n=3 (8×8): exactly three 1/3s in each row'
    ))
for col_idx, (P_show, labs) in enumerate([
    (P_hc2, ['00','01','10','11']),
    (P_hc3, [f'{i:03b}' for i in range(8)])
], 1):
    fig_hc.add_trace(go.Heatmap(
        z=P_show,
        colorscale=colorscale_ex, zmin=0, zmax=0.6,
        text=[[f'{v:.2f}' if v>0 else '' for v in row] for row in P_show],
        texttemplate='%{text}', textfont=dict(size=10),
        xgap=1, ygap=1, showscale=(col_idx==2),
        x=labs, y=labs
    ), row=1, col=col_idx)
    fig_hc.update_yaxes(autorange='reversed', row=1, col=col_idx)
fig_hc.update_layout(
    title='Transition matrices of the hypercube walk: n=2 and n=3',
    plot_bgcolor=C_BG, paper_bgcolor=C_BG,
    height=400, width=720,
    margin=dict(l=60, r=40, t=60, b=40)
)
fig_hc.show()
print("  ▸ White = 0 (Hamming distance > 1, unreachable)")
print("  ▸ Each row has exactly n nonzero entries (the n adjacent vertices), each equal to 1/n")
print("  ▸ The matrix is very sparse: fraction of nonzero entries = n/2ⁿ")

# ============================================================
print_header("§6 Hypercube walk | eigenvalue multiplicities and cutoff")
# ============================================================

from math import comb

for n_hc in [2, 3, 4]:
    print(f"\n▶ Hypercube dimension n = {n_hc} (number of states 2ⁿ = {1<<n_hc})")
    print("-" * 40)
    P_hc = make_hypercube(n=n_hc, r=0.0)
    eigs_hc = np.linalg.eigvals(P_hc).real

    # --- Step 1: eigenvalues vs the formula ---
    eigs_formula = {1 - 2*j/n_hc: comb(n_hc, j) for j in range(n_hc+1)}
    print(f"  Eigenvalues (formula) vs numerical computation:")
    print(f"  {'λ_j = 1-2j/n':>16}  {'theory C(n,j)':>16}  {'numerical':>10}")
    print(f"  " + "-" * 46)
    for lam_val, mult in sorted(eigs_formula.items(), reverse=True):
        num_mult = int(np.sum(np.abs(eigs_hc - lam_val) < 0.01))
        print(f"  {lam_val:>16.4f}  {mult:>16}  {num_mult:>10}")

    # --- Step 2: mixing time vs theory ---
    pi_hc = stationary_dist(P_hc)
    lam_hc = lambda1_abs(P_hc)
    tmix_hc = mixing_time(lam_hc)
    tmix_theory = (n_hc / 2) * np.log(n_hc) if n_hc > 1 else 1
    compare_print(
        f"t_mix (n={n_hc})",
        tmix_hc,
        f"theory (n/2)·ln n = {tmix_theory:.1f}"
    )

§7 The Perturbation Experiment: What If Every Step Carries a Little Noise?¶

Motivation: The Real World Is Not Perfect¶

The Markov chains studied so far are all “perfect”—their transition probabilities are exact and never change. In reality, however, weather forecasts are wrong, links break, and transmissions are noisy.

If a little random perturbation is added at each transition, how does the convergence behavior change?

$\varepsilon$ -Uniform Mixing Perturbation¶

We use the simplest kind of perturbation:

\tilde{\mathbf{P}} = (1-\varepsilon)\,\mathbf{P} + \varepsilon\,\mathbf{U}

(21)

$\varepsilon = 0$ : no perturbation at all; this is the original chain
$\varepsilon = 1$ : completely random—where each step goes is decided by $\mathbf{U}$ , independently of $\mathbf{P}$
$0 < \varepsilon < 1$ : follow the original rule with “confidence” $(1-\varepsilon)$ , and jump around at random with “probability” $\varepsilon$

Each row of $\mathbf{U}$ is a random probability vector sampled uniformly from the simplex, which guarantees that $\tilde{\mathbf{P}}$ is still a valid stochastic matrix (row sums equal 1).

Two Different Perturbation Models¶

There is a subtle but important distinction here:

Model 1 (a single perturbation): sample $\mathbf{U}$ once, keep it fixed, and then compute $\tilde{\mathbf{P}}^k$ .

Effect: this is just a new Markov chain whose transition matrix happens to be $\tilde{\mathbf{P}}$ rather than $\mathbf{P}$ . Because $\tilde{\mathbf{P}}$ contains a uniform component, its spectral gap is larger, and it actually converges faster.

Model 2 (injection at every step): at every step, sample a new $\mathbf{U}^{(k)}$ :

\boldsymbol{\mu}^{(k+1)} = \boldsymbol{\mu}^{(k)} \mathbf{P}^{(k)}, \quad \mathbf{P}^{(k)} = (1-\varepsilon)\mathbf{P} + \varepsilon \mathbf{U}^{(k)}

(22)

Effect: this is a random matrix product. Every path is different, and the TV distance curve “jitters.” In expectation, however, $\mathbb{E}[\mathbf{P}^{(k)}] = (1-\varepsilon)\mathbf{P} + \varepsilon (1/N)\mathbf{1}\mathbf{1}^\top$ , so the expected rate of convergence is still determined by $(1-\varepsilon)|\lambda_1|$ .

How to Read the Heat-Map Comparison¶

Top row (Model 1): orderly convergence, like a deterministic chain
Bottom row (Model 2): every resampling gives a different path, and the heat maps look “jittery”
Click the “Resample U ↻” button: Model 1 gets a new $\mathbf{U}$ ; Model 2 reruns the whole path

# ============================================================
# §7 Perturbation experiment: interactive dashboard
# ============================================================

# --- Step 1: uniformly random stochastic matrix (each row ~ Dir(1,...,1)) ---
def rand_stochastic(n):
    """Generate an n×n random stochastic matrix, each row uniformly distributed on the simplex"""
    G = -np.log(np.random.rand(n, n) + 1e-15)
    return G / G.sum(axis=1, keepdims=True)

def perturb_once(P, eps):
    """Model 1: a single perturbation"""
    U = rand_stochastic(P.shape[0])
    return (1 - eps) * P + eps * U

def simulate_model2(P, eps, max_k, n_paths=5):
    """Model 2: injection at every step; returns n_paths TV distance trajectories"""
    pi = stationary_dist(P)
    N  = P.shape[0]
    paths = []
    for _ in range(n_paths):
        mu = np.zeros(N); mu[0] = 1.0
        tv_path = []
        for _ in range(max_k + 1):
            tv_path.append(tv_dist(mu, pi))
            Pk = (1 - eps) * P + eps * rand_stochastic(N)
            mu = mu @ Pk
        paths.append(tv_path)
    return np.array(paths)

# --- Step 2: build the perturbation dashboard ---
def build_perturbation_dashboard(P, chain_name, eps, color):
    pi    = stationary_dist(P)
    lam0  = lambda1_abs(P)
    tmix0 = mixing_time(lam0)
    N     = P.shape[0]
    _tm0 = tmix0 if tmix0 is not None else 50
    max_k = min(_tm0 * 5 + 30, 300)
    _ts = tmix0 if tmix0 is not None else 30
    k_steps = [1, max(2, _ts//3), _ts, min(_ts*3, 400)]

    # Model 1
    P1    = perturb_once(P, eps)
    pi1   = stationary_dist(P1)
    lam1  = lambda1_abs(P1)
    tv1   = tv_decay_series(P1, pi, max_k)   # measured against the original π
    tv0   = tv_decay_series(P,  pi, max_k)

    # Model 2 (5 paths)
    tv2_paths = simulate_model2(P, eps, max_k, n_paths=5)
    tv2_mean  = tv2_paths.mean(axis=0)

    # heat maps
    colorscale = [[0.0,"#FFFFFF"],[0.5,"#AB82C5"],[1.0,"#57068C"]]

    fig = make_subplots(
        rows=3, cols=4,
        row_heights=[0.35, 0.35, 0.30],
        subplot_titles=(
            f"Model 1 P̃^k: k={k_steps[0]}", f"k={k_steps[1]}",
            f"k={k_steps[2]}",              f"k={k_steps[3]}",
            f"Model 2 cumulative product: k={k_steps[0]}", f"k={k_steps[1]}",
            f"k={k_steps[2]}",              f"k={k_steps[3]}",
            "TV distance: base vs Model 1", "TV distance: base vs Model 2", "", ""
        ),
        specs=[
            [{"type":"heatmap"},{"type":"heatmap"},{"type":"heatmap"},{"type":"heatmap"}],
            [{"type":"heatmap"},{"type":"heatmap"},{"type":"heatmap"},{"type":"heatmap"}],
            [{"type":"scatter"},{"type":"scatter"},{"type":"scatter"},{"type":"scatter"}],
        ]
    )

    # Model 1 heat maps (row=1)
    for idx, k in enumerate(k_steps):
        Pk1 = mat_pow(P1, k)
        fig.add_trace(go.Heatmap(
            z=Pk1, colorscale=colorscale, zmin=0, zmax=1,
            showscale=False, xgap=0.5, ygap=0.5
        ), row=1, col=idx+1)

    # Model 2 heat maps (row=2): cumulative product step by step
    prod = np.eye(N)
    k_prod_cache = {0: prod.copy()}
    for step in range(1, max(k_steps)+1):
        Pk_step = (1 - eps) * P + eps * rand_stochastic(N)
        prod = prod @ Pk_step
        if step in k_steps:
            k_prod_cache[step] = prod.copy()
    for idx, k in enumerate(k_steps):
        fig.add_trace(go.Heatmap(
            z=k_prod_cache.get(k, np.eye(N)),
            colorscale=colorscale, zmin=0, zmax=1,
            showscale=False, xgap=0.5, ygap=0.5
        ), row=2, col=idx+1)

    # TV curves (row=3, col=1): base vs Model 1
    ks = np.arange(max_k + 1)
    fig.add_trace(go.Scatter(
        x=ks, y=tv0, mode='lines',
        line=dict(color=C_V1, width=2),
        name='Base chain', showlegend=True
    ), row=3, col=1)
    fig.add_trace(go.Scatter(
        x=ks, y=tv1, mode='lines',
        line=dict(color=C_V2, width=2, dash='dash'),
        name='Model 1 P̃', showlegend=True
    ), row=3, col=1)
    fig.add_trace(go.Scatter(
        x=[0, max_k], y=[0.25, 0.25], mode='lines',
        line=dict(color=C_GRID, width=1, dash='dot'),
        showlegend=False
    ), row=3, col=1)

    # TV curves (row=3, col=2): base vs Model 2
    path_colors = ['rgba(255,93,71,0.5)','rgba(255,93,71,0.35)',
                   'rgba(255,93,71,0.25)','rgba(255,93,71,0.2)',
                   'rgba(255,93,71,0.15)']
    for i, path in enumerate(tv2_paths):
        fig.add_trace(go.Scatter(
            x=ks[:len(path)], y=path, mode='lines',
            line=dict(color=path_colors[i], width=1),
            showlegend=(i==0), name='Model 2 paths'
        ), row=3, col=2)
    fig.add_trace(go.Scatter(
        x=ks[:len(tv2_mean)], y=tv2_mean, mode='lines',
        line=dict(color=C_AXIS, width=1.5, dash='dash'),
        name='Model 2 expectation', showlegend=True
    ), row=3, col=2)
    fig.add_trace(go.Scatter(
        x=ks, y=tv0, mode='lines',
        line=dict(color=C_V1, width=2),
        name='Base chain', showlegend=False
    ), row=3, col=2)
    fig.add_trace(go.Scatter(
        x=[0, max_k], y=[0.25, 0.25], mode='lines',
        line=dict(color=C_GRID, width=1, dash='dot'),
        showlegend=False
    ), row=3, col=2)

    fig.update_layout(
        title=dict(
            text=(
                f"<b>Perturbation experiment: {chain_name}</b>&nbsp;&nbsp;ε = {eps:.2f}&nbsp;&nbsp;"
                f"|λ₁| base={lam0:.3f}&nbsp;&nbsp;Model 1={lam1:.3f}"
            ),
            font=dict(size=13, color=C_AXIS)
        ),
        plot_bgcolor=C_BG, paper_bgcolor=C_BG,
        height=680, width=900,
        margin=dict(l=40, r=20, t=80, b=40),
        legend=dict(x=0.5, y=-0.02, orientation='h')
    )
    for r_idx in [1, 2]:
        for c_idx in range(1, 5):
            fig.update_yaxes(autorange='reversed', row=r_idx, col=c_idx)
    for c_idx in [1, 2]:
        fig.update_xaxes(title_text='Step k', row=3, col=c_idx, gridcolor=C_GRID)
        fig.update_yaxes(title_text='TV distance', range=[0, 0.55],
                         row=3, col=c_idx, gridcolor=C_GRID)
    return fig

print_header("§7 Perturbation experiment | functions defined")

# ============================================================
# §7 Interactive controls for the perturbation experiment
# ============================================================

pert_chain_selector = widgets.ToggleButtons(
    options=list(CHAIN_CONFIGS.keys()),
    value='Birth-Death',
    description='Markov chain:',
    style={'description_width': '80px', 'button_width': '110px'},
)

eps_slider = widgets.FloatSlider(
    value=0.15, min=0.0, max=0.5, step=0.01,
    description='Perturbation ε:',
    style={'description_width': '110px'},
    readout_format='.2f',
    continuous_update=False
)

resample_btn = widgets.Button(
    description='Resample U ↻',
    button_style='',
    layout=widgets.Layout(width='160px')
)

pert_param_box   = widgets.VBox([])
pert_output      = widgets.Output()
pert_metrics_html = widgets.HTML(value='')

def update_pert_params(_=None):
    cfg = CHAIN_CONFIGS[pert_chain_selector.value]
    pert_param_box.children = cfg['params']
    for sl in cfg['params']:
        sl.observe(update_pert_figure, names='value')
    update_pert_figure()

def update_pert_figure(_=None):
    cfg  = CHAIN_CONFIGS[pert_chain_selector.value]
    ps   = [sl.value for sl in cfg['params']]
    P    = cfg['builder'](ps)
    eps  = eps_slider.value
    pi   = stationary_dist(P)
    lam0 = lambda1_abs(P)
    tmix = mixing_time(lam0)
    pert_metrics_html.value = (
        f"<div style='font-family:monospace; font-size:13px; "
        f"background:#FFF0ED; padding:8px 16px; border-radius:6px; margin:4px 0;'>"
        f"Matrix size: <b>{P.shape[0]}×{P.shape[0]}</b>&emsp;"
        f"base |λ₁| = <b>{lam0:.4f}</b>&emsp;"
        f"base t_mix = <b>{tmix}</b> steps&emsp;"
        f"ε = <b>{eps:.2f}</b>"
        f"</div>"
    )
    fig = build_perturbation_dashboard(
        P, pert_chain_selector.value, eps, cfg['color']
    )
    with pert_output:
        pert_output.clear_output(wait=True)
        fig.show()

pert_chain_selector.observe(update_pert_params, names='value')
eps_slider.observe(update_pert_figure, names='value')
resample_btn.on_click(update_pert_figure)

update_pert_params()

display(widgets.VBox([
    pert_chain_selector,
    pert_param_box,
    widgets.HBox([eps_slider, resample_btn]),
    pert_metrics_html,
    pert_output,
]))

Questions for Exploration¶

§2 Birth-Death Chain¶

Change $p$ from 0.5 to 0.9. How does the mixing time change? Intuitively, when $p = 0.9$ the system is “biased to one side,” and the stationary distribution is concentrated at the right end; why does this instead make the convergence faster?
Look at the eigenvalue plane for $p = 0.5$ : all the eigenvalues lie on the real axis. Move $p$ away from 0.5 (so that $\mathbf{P}$ is no longer symmetric): do the eigenvalues still lie on the real axis? What does this tell us, compared with the cycle chain (a symmetric matrix, whose eigenvalues are all real)?

§3 Cycle Chain¶

Set $N = 8$ and $r = 0$ , and look at the TV distance curve: does it oscillate? Gradually increase $r$ . When does the oscillation disappear? To the change of which eigenvalue does this correspond?
Fix $N = 8$ and compare the mixing times for $r = 0$ and $r = 0.1$ — adding a little “laziness” actually makes the convergence faster. Can you explain why?

§4 Google Matrix¶

Slowly change the damping factor $\alpha$ from 0.99 to 0.5 and observe:
- How do the heat maps change from sparse to uniform?
- How does PageRank (the stationary distribution) go from “nonuniform” toward “uniform”?
- How does the mixing time change?
If $\alpha = 1$ (no teleportation at all), the Google Matrix degenerates into the pure BA graph. Does a stationary distribution still exist in this case? What anomaly can you see in the eigenvalue plane?

§5 Doubly Stochastic Matrix¶

No matter how the mixing strength $m$ is adjusted, the stationary distribution is always uniform— give a one-line algebraic argument using the property that “the column sums equal 1.”
Set $m$ to 1.0; the matrix approaches a permutation matrix. How does the mixing time change? How can this be explained by the Birkhoff–von Neumann theorem?

§6 Hypercube¶

Compare the mixing times for $n = 2, 3, 4$ , compute their ratios, and check them against the theoretical prediction $t_{\text{mix}} \sim (n/2)\ln n$ . Do the numbers agree?
The hypercube has only $n+1$ distinct eigenvalues, yet the matrix is $2^n \times 2^n$ . What symmetry of this “highly degenerate” eigenvalue distribution can you see in the heat maps?

§7 Perturbation Experiment¶

Set $\varepsilon = 0$ . Do the TV curves of the two models coincide? Why?
Set $\varepsilon = 0.5$ . How does the “jitter amplitude” among the paths of Model 2 compare with that for $\varepsilon = 0.1$ ? What about the position of the expectation curve (the black dashed line)?
The TV curve of Model 1 usually converges faster than that of the base chain (violet)— because the spectral gap of $\tilde{\mathbf{P}}$ is larger. Can you compute a theoretical lower bound for the new spectral gap?

Chapter 8 Experiment 8: Mixing Times of Markov Chains and Matrix Structure

How to Use This Experiment¶

Structure of the Experiment¶

§0 Common Theoretical Framework¶

§0.1 What Is a Transition Matrix?¶

§0.2 What Is a Stationary Distribution?¶

§0.3 Eigenvalues and the Speed of Convergence¶

§0.4 The Spectral Gap and the Mixing Time¶

§0.5 Mathematical Utility Functions¶

§0.6 Constructing the Five Transition Matrices¶

§1 Unified Comparison Dashboard¶

How to Read the Seven Panels¶

The Metrics Bar¶

§2 The Birth-Death Chain¶

Physical Intuition: Hospital Beds¶

Matrix Structure¶

Detailed Balance: Why Are All the Eigenvalues Real?¶

A Concrete Example (First, What the Matrix Looks Like)¶

§3 The Cycle Chain¶

Physical Intuition: A Circular Track¶

Why Does the DFT Diagonalize It?¶

The Periodicity Problem of Even-Length Rings¶

A Concrete Example (First, What the Matrix Looks Like)¶

§4 Google Matrix (PageRank)¶

Motivation: How Do You Judge How Important a Web Page Is?¶

The Underlying Structure: An Irregular Directed Graph¶

Two Technical Problems and Google’s Solutions¶

The Mathematical Effect: A Rank-One Correction¶

The Trade-off between Speed and Quality¶

A Concrete Example (First, What the Matrix Looks Like)¶

§5 Doubly Stochastic Matrices¶

Physical Intuition: A Balanced Logistics Network¶

Definition¶

Why Must the Stationary Distribution Be Uniform?¶

The Birkhoff–von Neumann Theorem¶

A Concrete Example (First, What the Matrix Looks Like)¶

§6 The Hypercube Walk¶

Physical Intuition: Randomly Flipping the Bits of a Password¶

Matrix Structure: The Kronecker Product¶

Binomial Multiplicities of the Eigenvalues¶

The Cutoff Phenomenon¶

A Concrete Example (First, What the Matrix Looks Like)¶

§7 The Perturbation Experiment: What If Every Step Carries a Little Noise?¶

Motivation: The Real World Is Not Perfect¶

ε\varepsilonε-Uniform Mixing Perturbation¶

Two Different Perturbation Models¶

How to Read the Heat-Map Comparison¶

Questions for Exploration¶

§2 Birth-Death Chain¶

§3 Cycle Chain¶

§4 Google Matrix¶

§5 Doubly Stochastic Matrix¶

§6 Hypercube¶

§7 Perturbation Experiment¶

$\varepsilon$ -Uniform Mixing Perturbation¶