|
| 1 | +/- |
| 2 | +Copyright (c) 2026 Dennj Osele. All rights reserved. |
| 3 | +Released under Apache 2.0 license as described in the file LICENSE. |
| 4 | +Authors: Dennj Osele |
| 5 | +-/ |
| 6 | +module |
| 7 | + |
| 8 | +public import Mathlib.LinearAlgebra.Matrix.Stochastic |
| 9 | +public import Mathlib.Analysis.Convex.StdSimplex |
| 10 | +public import Mathlib.Analysis.Normed.Lp.PiLp |
| 11 | +public import Mathlib.Order.Filter.AtTopBot.Archimedean |
| 12 | + |
| 13 | +/-! |
| 14 | +# Stationary distributions for stochastic matrices |
| 15 | +
|
| 16 | +This file proves that every row-stochastic matrix on a finite nonempty state space has a |
| 17 | +stationary distribution in the standard simplex, via Cesàro averaging. |
| 18 | +
|
| 19 | +## Main definitions |
| 20 | +
|
| 21 | +* `IsStationary`: a distribution `μ` is stationary for `P` if `μ ᵥ* P = μ`. |
| 22 | +* `cesaroAverage`: the Cesàro average of the iterates of a vector under a matrix. |
| 23 | +
|
| 24 | +## Main results |
| 25 | +
|
| 26 | +* `Matrix.rowStochastic.exists_stationary_distribution`: every row-stochastic matrix on a finite |
| 27 | + nonempty state space has a stationary distribution in the standard simplex. |
| 28 | +
|
| 29 | +## Tags |
| 30 | +
|
| 31 | +stochastic matrix, Markov chain, stationary distribution, Cesàro average |
| 32 | +-/ |
| 33 | + |
| 34 | +@[expose] public section |
| 35 | + |
| 36 | +open Finset Matrix ENNReal Filter |
| 37 | + |
| 38 | +variable {n : Type*} [Fintype n] [DecidableEq n] |
| 39 | + |
| 40 | +/-! ### Stationary distributions -/ |
| 41 | + |
| 42 | +section Stationary |
| 43 | + |
| 44 | +variable {R : Type*} [Semiring R] |
| 45 | + |
| 46 | +/-- A distribution `μ` is stationary for a matrix `P` if `μ ᵥ* P = μ`. -/ |
| 47 | +def IsStationary (μ : n → R) (P : Matrix n n R) : Prop := μ ᵥ* P = μ |
| 48 | + |
| 49 | +/-- A stationary distribution for `P` is stationary for every power `P ^ k`. -/ |
| 50 | +lemma IsStationary.pow {μ : n → R} {P : Matrix n n R} (h : IsStationary μ P) (k : ℕ) : |
| 51 | + IsStationary μ (P ^ k) := by |
| 52 | + change μ ᵥ* P ^ k = μ |
| 53 | + induction k with |
| 54 | + | zero => simp |
| 55 | + | succ k ih => rw [pow_succ, ← Matrix.vecMul_vecMul, ih, h] |
| 56 | + |
| 57 | +end Stationary |
| 58 | + |
| 59 | +/-! ### Cesàro averages -/ |
| 60 | + |
| 61 | +section CesaroAverage |
| 62 | + |
| 63 | +variable {R : Type*} [DivisionRing R] |
| 64 | + |
| 65 | +/-- The Cesàro average of the iterates of a vector under a matrix. -/ |
| 66 | +def cesaroAverage (x₀ : n → R) (P : Matrix n n R) (k : ℕ) : n → R := |
| 67 | + (k + 1 : R)⁻¹ • ∑ i ∈ Finset.range (k + 1), x₀ ᵥ* (P ^ i) |
| 68 | + |
| 69 | +end CesaroAverage |
| 70 | + |
| 71 | +/-! ### L¹ non-expansiveness for row-stochastic matrices -/ |
| 72 | + |
| 73 | +section L1Norm |
| 74 | + |
| 75 | +variable {M : Matrix n n ℝ} |
| 76 | + |
| 77 | +namespace Matrix.rowStochastic |
| 78 | + |
| 79 | +/-- Row-stochastic matrices are non-expansive on `ℝⁿ` in the L¹ norm. -/ |
| 80 | +theorem nnnorm_vecMul_le (hM : M ∈ rowStochastic ℝ n) (z : n → ℝ) : |
| 81 | + ‖WithLp.toLp 1 (z ᵥ* M)‖₊ ≤ ‖WithLp.toLp 1 z‖₊ := by |
| 82 | + rw [← NNReal.coe_le_coe] |
| 83 | + simp only [coe_nnnorm, PiLp.norm_eq_of_L1, Real.norm_eq_abs] |
| 84 | + calc ∑ j, |(WithLp.toLp 1 (z ᵥ* M)) j| |
| 85 | + = ∑ j, |∑ i, z i * M i j| := by simp [vecMul, dotProduct] |
| 86 | + _ ≤ ∑ j, ∑ i, |z i| * M i j := sum_le_sum fun j _ ↦ |
| 87 | + (abs_sum_le_sum_abs _ _).trans <| sum_le_sum fun i _ ↦ by |
| 88 | + rw [abs_mul, abs_of_nonneg (hM.1 i j)] |
| 89 | + _ = ∑ i, |(WithLp.toLp 1 z) i| := by |
| 90 | + simp [sum_comm, ← mul_sum, sum_row_of_mem_rowStochastic hM] |
| 91 | + |
| 92 | +end Matrix.rowStochastic |
| 93 | + |
| 94 | +end L1Norm |
| 95 | + |
| 96 | +/-! ### Existence of stationary distributions -/ |
| 97 | + |
| 98 | +section Existence |
| 99 | + |
| 100 | +variable {M : Matrix n n ℝ} |
| 101 | + |
| 102 | +/-- The Cesàro average of a probability vector under a row-stochastic matrix |
| 103 | +belongs to the standard simplex. -/ |
| 104 | +lemma cesaroAverage_mem_stdSimplex (hM : M ∈ Matrix.rowStochastic ℝ n) |
| 105 | + {x₀ : n → ℝ} (hx₀ : x₀ ∈ stdSimplex ℝ n) (k : ℕ) : |
| 106 | + cesaroAverage x₀ M k ∈ stdSimplex ℝ n := by |
| 107 | + rw [cesaroAverage, Finset.smul_sum] |
| 108 | + refine (convex_stdSimplex ℝ n).sum_mem (fun _ _ ↦ by positivity) |
| 109 | + (by simp [mul_inv_cancel₀ (by positivity : (k + 1 : ℝ) ≠ 0)]) |
| 110 | + (fun i _ ↦ vecMul_mem_stdSimplex (Submonoid.pow_mem _ hM i) hx₀) |
| 111 | + |
| 112 | +/-- The Cesàro average is almost invariant: applying the matrix changes it by at most `2/(k+1)`. -/ |
| 113 | +lemma norm_cesaroAverage_vecMul_sub_le (hM : M ∈ Matrix.rowStochastic ℝ n) |
| 114 | + {x₀ : n → ℝ} (hx₀ : x₀ ∈ stdSimplex ℝ n) (k : ℕ) : |
| 115 | + ‖WithLp.toLp 1 (cesaroAverage x₀ M k ᵥ* M - cesaroAverage x₀ M k)‖ ≤ 2 / (k + 1) := by |
| 116 | + have hk : (0 : ℝ) < k + 1 := by positivity |
| 117 | + have hsimp : cesaroAverage x₀ M k ᵥ* M - cesaroAverage x₀ M k = |
| 118 | + (k + 1 : ℝ)⁻¹ • (x₀ ᵥ* M ^ (k + 1) - x₀) := by |
| 119 | + unfold cesaroAverage |
| 120 | + rw [Matrix.smul_vecMul, ← smul_sub, Matrix.sum_vecMul, ← Finset.sum_sub_distrib] |
| 121 | + simp_rw [Matrix.vecMul_vecMul, ← pow_succ] |
| 122 | + rw [Finset.sum_range_sub (fun i ↦ x₀ ᵥ* M ^ i), pow_zero, Matrix.vecMul_one] |
| 123 | + have nn1 : ∀ {x : n → ℝ}, x ∈ stdSimplex ℝ n → ‖WithLp.toLp 1 x‖₊ = 1 := fun hx ↦ by |
| 124 | + rw [← NNReal.coe_inj, NNReal.coe_one, coe_nnnorm, PiLp.norm_eq_of_L1] |
| 125 | + simp [Real.norm_eq_abs, Finset.sum_congr rfl fun i _ ↦ abs_of_nonneg (hx.1 i), hx.2] |
| 126 | + have hb : ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1) - x₀)‖₊ ≤ 2 := by |
| 127 | + rw [WithLp.toLp_sub]; refine (nnnorm_sub_le _ _).trans_eq ?_ |
| 128 | + rw [nn1 (vecMul_mem_stdSimplex (Submonoid.pow_mem _ hM (k + 1)) hx₀), nn1 hx₀]; norm_num |
| 129 | + rw [hsimp, WithLp.toLp_smul, norm_smul, Real.norm_eq_abs, abs_of_pos (inv_pos.mpr hk), |
| 130 | + div_eq_inv_mul] |
| 131 | + gcongr; exact_mod_cast hb |
| 132 | + |
| 133 | +variable [Nonempty n] |
| 134 | + |
| 135 | +namespace Matrix.rowStochastic |
| 136 | + |
| 137 | +/-- Every row-stochastic matrix on a finite nonempty state space has a stationary distribution. -/ |
| 138 | +theorem exists_stationary_distribution (hM : M ∈ rowStochastic ℝ n) : |
| 139 | + ∃ μ : n → ℝ, μ ∈ stdSimplex ℝ n ∧ IsStationary μ M := by |
| 140 | + let x₀ : n → ℝ := fun _ ↦ (Fintype.card n : ℝ)⁻¹ |
| 141 | + have hx₀ : x₀ ∈ stdSimplex ℝ n := |
| 142 | + ⟨fun _ ↦ by positivity, by simp [x₀, Finset.card_univ, nsmul_eq_mul]⟩ |
| 143 | + obtain ⟨μ, hμ, nₖ, hmono, hlim⟩ := (isCompact_stdSimplex ℝ n).tendsto_subseq |
| 144 | + fun k ↦ cesaroAverage_mem_stdSimplex hM hx₀ k |
| 145 | + have hcont : Continuous fun y : n → ℝ ↦ ‖WithLp.toLp 1 (y ᵥ* M - y)‖ := by fun_prop |
| 146 | + refine ⟨μ, hμ, sub_eq_zero.mp <| (WithLp.toLp_injective 1).eq_iff.mp <| norm_eq_zero.mp <| |
| 147 | + tendsto_nhds_unique ((hcont.tendsto μ).comp hlim) |
| 148 | + (squeeze_zero (fun _ ↦ norm_nonneg _) |
| 149 | + (fun k ↦ norm_cesaroAverage_vecMul_sub_le hM hx₀ (nₖ k)) |
| 150 | + ((tendsto_const_nhds.div_atTop tendsto_id).comp |
| 151 | + ((tendsto_natCast_atTop_atTop.comp hmono.tendsto_atTop).atTop_add tendsto_const_nhds)))⟩ |
| 152 | + |
| 153 | +end Matrix.rowStochastic |
| 154 | + |
| 155 | +end Existence |
0 commit comments