Skip to content

Commit 3517c37

Browse files
dennjdennj
committed
feat(Probability/Markov): stationary distributions for stochastic matrices
This PR proves that every row-stochastic matrix on a finite nonempty state space has a stationary distribution in the standard simplex. Main additions to `Mathlib/Probability/Markov/Stationary.lean`: - `IsStationary`: A distribution μ is stationary for matrix P if μ ᵥ* P = μ - `cesaroAverage`: Cesàro average of iterates of a vector under a matrix - `Matrix.rowStochastic.exists_stationary_distribution`: existence theorem The proof uses Cesàro averaging: start with uniform distribution, form averages, extract convergent subsequence by compactness, show limit is stationary via L¹ non-expansiveness. Also adds `vecMul_mem_stdSimplex` to `Stochastic.lean`: multiplying a probability vector by a row-stochastic matrix preserves simplex membership.
1 parent c7a9524 commit 3517c37

3 files changed

Lines changed: 216 additions & 0 deletions

File tree

‎Mathlib.lean‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6151,6 +6151,7 @@ public import Mathlib.Probability.Kernel.RadonNikodym
61516151
public import Mathlib.Probability.Kernel.Representation
61526152
public import Mathlib.Probability.Kernel.SetIntegral
61536153
public import Mathlib.Probability.Kernel.WithDensity
6154+
public import Mathlib.Probability.Markov.Stationary
61546155
public import Mathlib.Probability.Martingale.Basic
61556156
public import Mathlib.Probability.Martingale.BorelCantelli
61566157
public import Mathlib.Probability.Martingale.Centering

‎Mathlib/LinearAlgebra/Matrix/Stochastic.lean‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@ module
88
public import Mathlib.Data.Matrix.Basic
99
public import Mathlib.Data.Matrix.Mul
1010
public import Mathlib.Analysis.Convex.Basic
11+
public import Mathlib.Analysis.Convex.StdSimplex
1112
public import Mathlib.LinearAlgebra.Matrix.Permutation
1213

1314
/-!
@@ -258,4 +259,13 @@ lemma reindex_mem_colStochastic_iff {m : Type*} [Fintype m] [DecidableEq m] {M :
258259
@[aesop safe apply]
259260
alias ⟨_, reindex_mem_colStochastic⟩ := reindex_mem_colStochastic_iff
260261

262+
/-- Multiplying a probability vector on the right by a row-stochastic matrix
263+
preserves membership in the standard simplex. -/
264+
lemma vecMul_mem_stdSimplex (hM : M ∈ rowStochastic R n)
265+
{v : n → R} (hv : v ∈ stdSimplex R n) : v ᵥ* M ∈ stdSimplex R n := by
266+
refine ⟨fun j => Finset.sum_nonneg fun i _ => mul_nonneg (hv.1 i) (hM.1 i j), ?_⟩
267+
simp only [vecMul, dotProduct]
268+
rw [Finset.sum_comm, ← hv.2]
269+
exact Finset.sum_congr rfl fun i _ => by simp [← mul_sum, sum_row_of_mem_rowStochastic hM]
270+
261271
end Matrix

‎Mathlib/Probability/Markov/Stationary.lean‎

Lines changed: 205 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,205 @@
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+
* `uniformDistribution`: the uniform distribution on a finite nonempty type.
24+
25+
## Main results
26+
27+
* `Matrix.rowStochastic.exists_stationary_distribution`: every row-stochastic matrix on a finite
28+
nonempty state space has a stationary distribution in the standard simplex.
29+
30+
## Tags
31+
32+
stochastic matrix, Markov chain, stationary distribution, Cesàro average
33+
-/
34+
35+
@[expose] public section
36+
37+
open Finset Function Matrix ENNReal Filter
38+
open scoped BigOperators
39+
40+
variable {n : Type*} [Fintype n] [DecidableEq n]
41+
42+
/-! ### Stationary distributions -/
43+
44+
section Stationary
45+
46+
variable {R : Type*} [Semiring R]
47+
48+
/-- A distribution `μ` is stationary for a matrix `P` if `μ ᵥ* P = μ`. -/
49+
def IsStationary (μ : n → R) (P : Matrix n n R) : Prop := μ ᵥ* P = μ
50+
51+
/-- A stationary distribution for `P` is stationary for every power `P ^ k`. -/
52+
lemma IsStationary.pow {μ : n → R} {P : Matrix n n R} (h : IsStationary μ P) (k : ℕ) :
53+
IsStationary μ (P ^ k) := by
54+
change μ ᵥ* P ^ k = μ
55+
induction k with
56+
| zero => simp
57+
| succ k ih => rw [pow_succ, ← Matrix.vecMul_vecMul, ih, h]
58+
59+
end Stationary
60+
61+
/-! ### Cesàro averages -/
62+
63+
section CesaroAverage
64+
65+
variable {R : Type*} [DivisionRing R] [PartialOrder R] [IsStrictOrderedRing R]
66+
67+
/-- The Cesàro average of the iterates of a vector under a matrix. -/
68+
def cesaroAverage (x₀ : n → R) (P : Matrix n n R) (k : ℕ) : n → R :=
69+
(k + 1 : R)⁻¹ • ∑ i ∈ Finset.range (k + 1), x₀ ᵥ* (P ^ i)
70+
71+
variable [Nonempty n]
72+
73+
/-- The uniform distribution on a finite nonempty type. -/
74+
@[reducible, nolint unusedArguments]
75+
def uniformDistribution (R : Type*) (n : Type*) [Fintype n] [DivisionRing R] :
76+
n → R := fun _ => 1 / Fintype.card n
77+
78+
end CesaroAverage
79+
80+
/-! ### L¹ non-expansiveness for row-stochastic matrices -/
81+
82+
section L1Norm
83+
84+
variable {M : Matrix n n ℝ}
85+
86+
omit [DecidableEq n] in
87+
private lemma l1_nnnorm_eq_sum (f : PiLp 1 (fun _ : n => ℝ)) : (‖f‖₊ : ℝ) = ∑ i, |f.ofLp i| := by
88+
rw [PiLp.nnnorm_eq_sum one_ne_top]
89+
simp only [ENNReal.toReal_one, NNReal.rpow_one, div_one, NNReal.coe_sum, coe_nnnorm,
90+
Real.norm_eq_abs]
91+
92+
omit [DecidableEq n] in
93+
private lemma l1_nnnorm_eq_one {x : n → ℝ} (hx : x ∈ stdSimplex ℝ n) : ‖WithLp.toLp 1 x‖₊ = 1 := by
94+
rw [← NNReal.coe_inj, NNReal.coe_one, l1_nnnorm_eq_sum,
95+
(sum_congr rfl fun i _ => abs_of_nonneg (hx.1 i) : ∑ i, |x i| = ∑ i, x i), hx.2]
96+
97+
namespace Matrix.rowStochastic
98+
99+
/-- Row-stochastic matrices are non-expansive on `ℝⁿ` in the L¹ norm. -/
100+
theorem nnnorm_vecMul_sub_le (hM : M ∈ rowStochastic ℝ n) (x y : n → ℝ) :
101+
‖WithLp.toLp 1 (x ᵥ* M - y ᵥ* M)‖₊ ≤ ‖WithLp.toLp 1 (x - y)‖₊ := by
102+
have hxy : x ᵥ* M - y ᵥ* M = fun j => ∑ i, (x i - y i) * M i j := by
103+
ext j; simp only [Pi.sub_apply, vecMul, dotProduct, sub_mul, sum_sub_distrib]
104+
have key : ∑ j, |∑ i, (x i - y i) * M i j| ≤ ∑ i, |x i - y i| := calc
105+
∑ j, |∑ i, (x i - y i) * M i j|
106+
≤ ∑ j, ∑ i, |x i - y i| * M i j := sum_le_sum fun j _ => (abs_sum_le_sum_abs _ _).trans
107+
(sum_le_sum fun i _ => by rw [abs_mul, abs_of_nonneg (hM.1 i j)])
108+
_ = ∑ i, |x i - y i| := by
109+
rw [sum_comm]; simp_rw [← mul_sum, sum_row_of_mem_rowStochastic hM, mul_one]
110+
have hnorm : (‖WithLp.toLp 1 (x ᵥ* M - y ᵥ* M)‖₊ : ℝ) = ∑ j, |∑ i, (x i - y i) * M i j| := by
111+
rw [l1_nnnorm_eq_sum]; exact sum_congr rfl fun j _ => congrArg abs (congrFun hxy j)
112+
exact NNReal.coe_le_coe.mp (by rw [hnorm, l1_nnnorm_eq_sum]; exact key)
113+
114+
/-- Row-stochastic matrices are non-expansive on `ℝⁿ` in the L¹ norm. -/
115+
theorem norm_vecMul_sub_le (hM : M ∈ rowStochastic ℝ n) (x y : n → ℝ) :
116+
‖WithLp.toLp 1 (x ᵥ* M - y ᵥ* M)‖ ≤ ‖WithLp.toLp 1 (x - y)‖ :=
117+
mod_cast nnnorm_vecMul_sub_le hM x y
118+
119+
end Matrix.rowStochastic
120+
121+
end L1Norm
122+
123+
/-! ### Existence of stationary distributions -/
124+
125+
section Existence
126+
127+
variable {M : Matrix n n ℝ}
128+
129+
/-- The Cesàro average of a probability vector under a row-stochastic matrix
130+
belongs to the standard simplex. -/
131+
lemma cesaroAverage_mem_stdSimplex (hM : M ∈ Matrix.rowStochastic ℝ n)
132+
{x₀ : n → ℝ} (hx₀ : x₀ ∈ stdSimplex ℝ n) (k : ℕ) :
133+
cesaroAverage x₀ M k ∈ stdSimplex ℝ n := by
134+
have hmem i : x₀ ᵥ* M ^ i ∈ stdSimplex ℝ n :=
135+
vecMul_mem_stdSimplex (Submonoid.pow_mem _ hM i) hx₀
136+
refine ⟨fun j => mul_nonneg (inv_nonneg.mpr (by linarith : (0 : ℝ) ≤ k + 1))
137+
(by simp only [Finset.sum_apply]; exact sum_nonneg fun i _ => (hmem i).1 j), ?_⟩
138+
simp only [cesaroAverage, Pi.smul_apply, smul_eq_mul, Finset.sum_apply, ← mul_sum,
139+
sum_comm (γ := n)]
140+
rw [show ∑ i ∈ range (k + 1), ∑ j, (x₀ ᵥ* M ^ i) j = k + 1 from
141+
(sum_congr rfl fun i _ => (hmem i).2).trans (by simp), mul_comm]
142+
exact mul_inv_cancel₀ (by linarith : (k : ℝ) + 1 ≠ 0)
143+
144+
/-- The Cesàro average is almost invariant: applying the matrix changes it by at most `2/(k+1)`. -/
145+
lemma norm_cesaroAverage_vecMul_sub_le (hM : M ∈ Matrix.rowStochastic ℝ n)
146+
{x₀ : n → ℝ} (hx₀ : x₀ ∈ stdSimplex ℝ n) (k : ℕ) :
147+
‖WithLp.toLp 1 (cesaroAverage x₀ M k ᵥ* M - cesaroAverage x₀ M k)‖ ≤ 2 / (k + 1) := by
148+
have hk : (0 : ℝ) < k + 1 := by linarith
149+
have hsimp : cesaroAverage x₀ M k ᵥ* M - cesaroAverage x₀ M k =
150+
(k + 1 : ℝ)⁻¹ • (x₀ ᵥ* M ^ (k + 1) - x₀) := by
151+
unfold cesaroAverage
152+
rw [Matrix.smul_vecMul, ← smul_sub, Matrix.sum_vecMul]; congr 1
153+
have h1 : ∀ i, (x₀ ᵥ* M ^ i) ᵥ* M = x₀ ᵥ* M ^ (i + 1) := fun i => by
154+
rw [Matrix.vecMul_vecMul, ← pow_succ]
155+
simp_rw [h1]
156+
rw [Finset.sum_range_succ' (fun i => x₀ ᵥ* M ^ i), Finset.sum_range_succ, pow_zero,
157+
Matrix.vecMul_one]; abel
158+
rw [hsimp, WithLp.toLp_smul, norm_smul, Real.norm_eq_abs, abs_of_pos (inv_pos.mpr hk)]
159+
have h1 := vecMul_mem_stdSimplex (Submonoid.pow_mem _ hM (k + 1)) hx₀
160+
have : ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1) - x₀)‖ ≤ 2 := calc
161+
‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1) - x₀)‖₊
162+
_ = ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1)) - WithLp.toLp 1 x₀‖₊ := by rw [← WithLp.toLp_sub]
163+
_ ≤ ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1))‖₊ + ‖WithLp.toLp 1 x₀‖₊ := nnnorm_sub_le _ _
164+
_ = 2 := by rw [l1_nnnorm_eq_one h1, l1_nnnorm_eq_one hx₀]; norm_num
165+
calc (k + 1 : ℝ)⁻¹ * ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1) - x₀)‖
166+
_ ≤ (k + 1 : ℝ)⁻¹ * 2 := by gcongr
167+
_ = 2 / (k + 1) := by ring
168+
169+
variable [Nonempty n]
170+
171+
omit [DecidableEq n] in
172+
/-- The uniform distribution belongs to the standard simplex. -/
173+
lemma uniformDistribution_mem_stdSimplex : uniformDistribution ℝ (n := n) ∈ stdSimplex ℝ n :=
174+
⟨fun _ => by simp only [uniformDistribution, one_div, inv_nonneg]; positivity,
175+
by simp [uniformDistribution, Finset.card_univ, nsmul_eq_mul]⟩
176+
177+
namespace Matrix.rowStochastic
178+
179+
/-- Every row-stochastic matrix on a finite nonempty state space has a stationary distribution. -/
180+
theorem exists_stationary_distribution (hM : M ∈ rowStochastic ℝ n) :
181+
∃ μ : n → ℝ, μ ∈ stdSimplex ℝ n ∧ IsStationary μ M := by
182+
let x₀ := uniformDistribution ℝ (n := n)
183+
let xₖ : ℕ → n → ℝ := fun k => cesaroAverage x₀ M k
184+
have hxₖ k : xₖ k ∈ stdSimplex ℝ n :=
185+
cesaroAverage_mem_stdSimplex hM uniformDistribution_mem_stdSimplex k
186+
obtain ⟨μ, hμ_mem, nₖ, hnₖ_mono, hnₖ_lim⟩ := (isCompact_stdSimplex ℝ n).tendsto_subseq hxₖ
187+
refine ⟨μ, hμ_mem, ?_⟩
188+
have h_lim_diff : Tendsto (fun k => xₖ (nₖ k) ᵥ* M - xₖ (nₖ k)) atTop (nhds (μ ᵥ* M - μ)) :=
189+
((Continuous.matrix_vecMul continuous_id continuous_const).continuousAt.tendsto.comp
190+
hnₖ_lim).sub hnₖ_lim
191+
have h_error_tendsto : Tendsto (fun k => ‖WithLp.toLp 1 (xₖ (nₖ k) ᵥ* M - xₖ (nₖ k))‖)
192+
atTop (nhds 0) :=
193+
tendsto_of_tendsto_of_tendsto_of_le_of_le tendsto_const_nhds
194+
((tendsto_const_nhds.div_atTop tendsto_id).comp
195+
((tendsto_natCast_atTop_atTop.comp hnₖ_mono.tendsto_atTop).atTop_add tendsto_const_nhds))
196+
(fun _ => norm_nonneg _)
197+
(fun k => norm_cesaroAverage_vecMul_sub_le hM uniformDistribution_mem_stdSimplex (nₖ k))
198+
have h_norm_zero := tendsto_nhds_unique
199+
(((PiLp.continuous_toLp 1 (fun _ => ℝ)).continuousAt.tendsto.comp h_lim_diff).norm)
200+
h_error_tendsto
201+
exact sub_eq_zero.mp ((WithLp.toLp_injective 1).eq_iff.mp (norm_eq_zero.mp h_norm_zero))
202+
203+
end Matrix.rowStochastic
204+
205+
end Existence

0 commit comments

Comments
 (0)