feat(LinearAlgebra/Matrix): add row-stochastic and column-stochastic matrices

Define `rowStochastic` and `colStochastic` as submonoids of square matrices,
and establish their basic properties including:

- Membership characterizations via row/column sums
- Closure under multiplication and powers
- Convexity of the sets
- Relationship with `doublyStochastic` (it equals their infimum)
- Transpose duality between row and column stochastic
- Connection to the standard simplex

Also prove the L¹ non-expansiveness property for row-stochastic matrices
(`rowStochastic.nnnorm_vecMul_sub_le`), which is fundamental to Markov chain
theory, and use it to show existence of stationary distributions via Cesàro
averages and compactness of the standard simplex.

Author: dennj

Mathlib.lean

@@ -4749,6 +4749,7 @@ public import Mathlib.LinearAlgebra.Matrix.ProjectiveSpecialLinearGroup
 public import Mathlib.LinearAlgebra.Matrix.Rank
 public import Mathlib.LinearAlgebra.Matrix.Reindex
 public import Mathlib.LinearAlgebra.Matrix.RowCol
+public import Mathlib.LinearAlgebra.Matrix.RowStochastic
 public import Mathlib.LinearAlgebra.Matrix.SchurComplement
 public import Mathlib.LinearAlgebra.Matrix.SemiringInverse
 public import Mathlib.LinearAlgebra.Matrix.SesquilinearForm

Mathlib/LinearAlgebra/Matrix/RowStochastic.lean

@@ -0,0 +1,389 @@
+/-
+Copyright (c) 2026 Dennj Osele. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+Authors: Dennj Osele
+-/
+module
+
+public import Mathlib.Analysis.Convex.DoublyStochasticMatrix
+public import Mathlib.Analysis.Convex.StdSimplex
+public import Mathlib.Analysis.Normed.Lp.PiLp
+public import Mathlib.Order.Filter.AtTopBot.Archimedean
+
+/-!
+# Row-Stochastic and Column-Stochastic Matrices
+
+This file defines row-stochastic and column-stochastic matrices, and establishes their
+basic properties including closure under multiplication and the L¹ non-expansiveness
+property that is fundamental to Markov chain theory.
+
+## Main definitions
+
+* `rowStochastic`: A square matrix is row-stochastic if all entries are nonnegative and
+  each row sums to 1.
+* `colStochastic`: A square matrix is column-stochastic if all entries are nonnegative and
+  each column sums to 1.
+* `IsStationary`: A distribution `μ` is stationary for a matrix `P` if `μ ᵥ* P = μ`.
+* `cesaroAverage`: The Cesàro average of the iterates of a vector under a matrix.
+* `uniformDistribution`: The uniform distribution on a finite nonempty type.
+
+## Main statements
+
+* `rowStochastic.pow_mem`: Powers of row-stochastic matrices are row-stochastic.
+* `rowStochastic.nnnorm_vecMul_sub_le`: Row-stochastic matrices are non-expansive
+  operators on probability vectors in the L¹ norm.
+* `rowStochastic.exists_stationary_distribution`: Every row-stochastic matrix on a finite
+  nonempty state space has a stationary distribution in the standard simplex.
+
+## Implementation notes
+
+Row-stochastic matrices form a submonoid of square matrices. They are the natural
+representation of transition matrices for discrete-time Markov chains, where `P i j`
+represents the probability of transitioning from state `i` to state `j`.
+
+The relationship with `doublyStochastic` is:
+`doublyStochastic R n = rowStochastic R n ⊓ colStochastic R n`
+
+## Tags
+
+stochastic matrix, Markov chain, row-stochastic, column-stochastic, L¹ contraction
+-/
+
+@[expose] public section
+
+open Finset Function Matrix ENNReal Filter
+open scoped BigOperators
+
+variable {R n : Type*} [Fintype n] [DecidableEq n]
+
+section OrderedSemiring
+
+variable [Semiring R] [PartialOrder R] [IsOrderedRing R] {M : Matrix n n R}
+
+/--
+A square matrix is row-stochastic iff all entries are nonnegative and each row sums to 1.
+Equivalently, right multiplication by the all-ones column vector gives the all-ones vector.
+-/
+def rowStochastic (R n : Type*) [Fintype n] [DecidableEq n] [Semiring R] [PartialOrder R]
+    [IsOrderedRing R] : Submonoid (Matrix n n R) where
+  carrier := {M | (∀ i j, 0 ≤ M i j) ∧ M *ᵥ 1 = 1}
+  mul_mem' {M N} hM hN :=
+    ⟨fun i j => sum_nonneg fun k _ => mul_nonneg (hM.1 _ _) (hN.1 _ _),
+     by rw [← mulVec_mulVec, hN.2, hM.2]⟩
+  one_mem' := by simp [zero_le_one_elem]
+
+/--
+A square matrix is column-stochastic iff all entries are nonnegative and each column sums to 1.
+Equivalently, left multiplication by the all-ones row vector gives the all-ones vector.
+-/
+def colStochastic (R n : Type*) [Fintype n] [DecidableEq n] [Semiring R] [PartialOrder R]
+    [IsOrderedRing R] : Submonoid (Matrix n n R) where
+  carrier := {M | (∀ i j, 0 ≤ M i j) ∧ 1 ᵥ* M = 1}
+  mul_mem' {M N} hM hN :=
+    ⟨fun i j => sum_nonneg fun k _ => mul_nonneg (hM.1 _ _) (hN.1 _ _),
+     by rw [← vecMul_vecMul, hM.2, hN.2]⟩
+  one_mem' := by simp [zero_le_one_elem]
+
+lemma mem_rowStochastic : M ∈ rowStochastic R n ↔ (∀ i j, 0 ≤ M i j) ∧ M *ᵥ 1 = 1 := Iff.rfl
+
+lemma mem_colStochastic : M ∈ colStochastic R n ↔ (∀ i j, 0 ≤ M i j) ∧ 1 ᵥ* M = 1 := Iff.rfl
+
+lemma mem_rowStochastic_iff_sum :
+    M ∈ rowStochastic R n ↔ (∀ i j, 0 ≤ M i j) ∧ ∀ i, ∑ j, M i j = 1 := by
+  simp [funext_iff, rowStochastic, mulVec, dotProduct]
+
+lemma mem_colStochastic_iff_sum :
+    M ∈ colStochastic R n ↔ (∀ i j, 0 ≤ M i j) ∧ ∀ j, ∑ i, M i j = 1 := by
+  simp [funext_iff, colStochastic, vecMul, dotProduct]
+
+/-- Every entry of a row-stochastic matrix is nonnegative. -/
+lemma nonneg_of_mem_rowStochastic (hM : M ∈ rowStochastic R n) {i j : n} : 0 ≤ M i j :=
+  hM.1 _ _
+
+/-- Every entry of a column-stochastic matrix is nonnegative. -/
+lemma nonneg_of_mem_colStochastic (hM : M ∈ colStochastic R n) {i j : n} : 0 ≤ M i j :=
+  hM.1 _ _
+
+/-- Each row sum of a row-stochastic matrix is 1. -/
+lemma sum_row_of_mem_rowStochastic (hM : M ∈ rowStochastic R n) (i : n) : ∑ j, M i j = 1 :=
+  (mem_rowStochastic_iff_sum.1 hM).2 _
+
+/-- Each column sum of a column-stochastic matrix is 1. -/
+lemma sum_col_of_mem_colStochastic (hM : M ∈ colStochastic R n) (j : n) : ∑ i, M i j = 1 :=
+  (mem_colStochastic_iff_sum.1 hM).2 _
+
+/-- A row-stochastic matrix multiplied with the all-ones column vector is 1. -/
+lemma mulVec_one_of_mem_rowStochastic (hM : M ∈ rowStochastic R n) : M *ᵥ 1 = 1 :=
+  (mem_rowStochastic.1 hM).2
+
+/-- The all-ones row vector multiplied with a column-stochastic matrix is 1. -/
+lemma one_vecMul_of_mem_colStochastic (hM : M ∈ colStochastic R n) : 1 ᵥ* M = 1 :=
+  (mem_colStochastic.1 hM).2
+
+/-- Every entry of a row-stochastic matrix is at most 1. -/
+lemma le_one_of_mem_rowStochastic (hM : M ∈ rowStochastic R n) {i j : n} : M i j ≤ 1 :=
+  sum_row_of_mem_rowStochastic hM i ▸ single_le_sum (fun k _ => hM.1 _ k) (mem_univ j)
+
+/-- Every entry of a column-stochastic matrix is at most 1. -/
+lemma le_one_of_mem_colStochastic (hM : M ∈ colStochastic R n) {i j : n} : M i j ≤ 1 :=
+  sum_col_of_mem_colStochastic hM j ▸ single_le_sum (fun k _ => hM.1 k _) (mem_univ i)
+
+/-- The set of row-stochastic matrices is convex. -/
+lemma convex_rowStochastic : Convex R (rowStochastic R n : Set (Matrix n n R)) := by
+  intro x hx y hy a b ha hb hab
+  simp only [SetLike.mem_coe, mem_rowStochastic_iff_sum] at hx hy ⊢
+  simp [add_nonneg, ha, hb, mul_nonneg, hx, hy, sum_add_distrib, ← mul_sum, hab]
+
+/-- The set of column-stochastic matrices is convex. -/
+lemma convex_colStochastic : Convex R (colStochastic R n : Set (Matrix n n R)) := by
+  intro x hx y hy a b ha hb hab
+  simp only [SetLike.mem_coe, mem_colStochastic_iff_sum] at hx hy ⊢
+  simp [add_nonneg, ha, hb, mul_nonneg, hx, hy, sum_add_distrib, ← mul_sum, hab]
+
+/-- Any permutation matrix is row-stochastic. -/
+lemma permMatrix_mem_rowStochastic {σ : Equiv.Perm n} : σ.permMatrix R ∈ rowStochastic R n := by
+  rw [mem_rowStochastic_iff_sum]
+  exact ⟨fun i j => by aesop, fun _ => by simp [Equiv.toPEquiv_apply]⟩
+
+/-- Any permutation matrix is column-stochastic. -/
+lemma permMatrix_mem_colStochastic {σ : Equiv.Perm n} : σ.permMatrix R ∈ colStochastic R n := by
+  rw [mem_colStochastic_iff_sum]
+  exact ⟨fun i j => by aesop, fun _ => by simp [Equiv.toPEquiv_apply, ← Equiv.eq_symm_apply]⟩
+
+/-- A matrix is doubly stochastic iff it is both row-stochastic and column-stochastic. -/
+lemma doublyStochastic_eq_inf :
+    doublyStochastic R n = rowStochastic R n ⊓ colStochastic R n := by
+  ext M
+  simp only [Submonoid.mem_inf, mem_doublyStochastic, mem_rowStochastic, mem_colStochastic]
+  exact ⟨fun ⟨h1, h2, h3⟩ => ⟨⟨h1, h2⟩, ⟨h1, h3⟩⟩, fun ⟨⟨h1, h2⟩, ⟨_, h3⟩⟩ => ⟨h1, h2, h3⟩⟩
+
+/-- The transpose of a row-stochastic matrix is column-stochastic. -/
+lemma transpose_mem_colStochastic (hM : M ∈ rowStochastic R n) : Mᵀ ∈ colStochastic R n := by
+  simp only [mem_colStochastic_iff_sum, transpose_apply]
+  exact ⟨fun i j => hM.1 j i, fun j => sum_row_of_mem_rowStochastic hM j⟩
+
+/-- The transpose of a column-stochastic matrix is row-stochastic. -/
+lemma transpose_mem_rowStochastic (hM : M ∈ colStochastic R n) : Mᵀ ∈ rowStochastic R n := by
+  simp only [mem_rowStochastic_iff_sum, transpose_apply]
+  exact ⟨fun i j => hM.1 j i, fun i => sum_col_of_mem_colStochastic hM i⟩
+
+end OrderedSemiring
+
+/-! ### Row-stochastic matrices and the standard simplex -/
+
+section Simplex
+
+variable [Semiring R] [PartialOrder R] [IsOrderedRing R]
+
+/-- Each row of a row-stochastic matrix belongs to the standard simplex. -/
+lemma row_mem_stdSimplex {M : Matrix n n R} (hM : M ∈ rowStochastic R n) (i : n) :
+    M i ∈ stdSimplex R n :=
+  ⟨fun j => hM.1 i j, sum_row_of_mem_rowStochastic hM i⟩
+
+/-- Each column of a column-stochastic matrix belongs to the standard simplex. -/
+lemma col_mem_stdSimplex {M : Matrix n n R} (hM : M ∈ colStochastic R n) (j : n) :
+    (fun i => M i j) ∈ stdSimplex R n :=
+  ⟨fun i => hM.1 i j, sum_col_of_mem_colStochastic hM j⟩
+
+/-- Multiplying a probability vector on the left by a row-stochastic matrix
+preserves membership in the standard simplex. -/
+lemma vecMul_mem_stdSimplex {M : Matrix n n R} (hM : M ∈ rowStochastic R n)
+    {v : n → R} (hv : v ∈ stdSimplex R n) : v ᵥ* M ∈ stdSimplex R n := by
+  refine ⟨fun j => sum_nonneg fun i _ => mul_nonneg (hv.1 i) (hM.1 i j), ?_⟩
+  simp only [vecMul, dotProduct]
+  rw [sum_comm, ← hv.2]
+  exact sum_congr rfl fun i _ => by simp [← mul_sum, sum_row_of_mem_rowStochastic hM]
+
+end Simplex
+
+/-! ### L¹ non-expansiveness for row-stochastic matrices -/
+
+section L1Norm
+
+variable {M : Matrix n n ℝ}
+
+omit [DecidableEq n] in
+/-- The L¹ norm of a function equals the sum of absolute values. -/
+private lemma l1_nnnorm_eq_sum (f : PiLp 1 (fun _ : n => ℝ)) : (‖f‖₊ : ℝ) = ∑ i, |f.ofLp i| := by
+  rw [PiLp.nnnorm_eq_sum one_ne_top]
+  simp only [ENNReal.toReal_one, NNReal.rpow_one, div_one, NNReal.coe_sum, coe_nnnorm,
+    Real.norm_eq_abs]
+
+namespace rowStochastic
+
+/-- Row-stochastic matrices are non-expansive operators on probability vectors in the L¹ norm.
+This is the key contraction property for Markov chains. -/
+theorem nnnorm_vecMul_sub_le (hM : M ∈ rowStochastic ℝ n) (x y : n → ℝ) :
+    ‖WithLp.toLp 1 (x ᵥ* M - y ᵥ* M)‖₊ ≤ ‖WithLp.toLp 1 (x - y)‖₊ := by
+  have hxy : x ᵥ* M - y ᵥ* M = fun j => ∑ i, (x i - y i) * M i j := by
+    ext j; simp only [Pi.sub_apply, vecMul, dotProduct, sub_mul, sum_sub_distrib]
+  have key : ∑ j, |∑ i, (x i - y i) * M i j| ≤ ∑ i, |x i - y i| := calc
+    ∑ j, |∑ i, (x i - y i) * M i j|
+      ≤ ∑ j, ∑ i, |x i - y i| * M i j := sum_le_sum fun j _ => (abs_sum_le_sum_abs _ _).trans
+        (sum_le_sum fun i _ => by rw [abs_mul, abs_of_nonneg (hM.1 i j)])
+    _ = ∑ i, |x i - y i| := by
+        rw [sum_comm]; simp_rw [← mul_sum, sum_row_of_mem_rowStochastic hM, mul_one]
+  have hnorm : (‖WithLp.toLp 1 (x ᵥ* M - y ᵥ* M)‖₊ : ℝ) = ∑ j, |∑ i, (x i - y i) * M i j| := by
+    rw [l1_nnnorm_eq_sum]; exact sum_congr rfl fun j _ => congrArg abs (congrFun hxy j)
+  exact NNReal.coe_le_coe.mp (by rw [hnorm, l1_nnnorm_eq_sum]; exact key)
+
+/-- Row-stochastic matrices are non-expansive in L¹ norm (version with norm instead of nnnorm). -/
+theorem norm_vecMul_sub_le (hM : M ∈ rowStochastic ℝ n) (x y : n → ℝ) :
+    ‖WithLp.toLp 1 (x ᵥ* M - y ᵥ* M)‖ ≤ ‖WithLp.toLp 1 (x - y)‖ :=
+  mod_cast nnnorm_vecMul_sub_le hM x y
+
+end rowStochastic
+
+end L1Norm
+
+/-! ### Powers of stochastic matrices -/
+
+namespace rowStochastic
+
+variable [Semiring R] [PartialOrder R] [IsOrderedRing R]
+
+/-- Powers of row-stochastic matrices are row-stochastic. -/
+theorem pow_mem {M : Matrix n n R} (hM : M ∈ rowStochastic R n) (k : ℕ) :
+    M ^ k ∈ rowStochastic R n :=
+  Submonoid.pow_mem (rowStochastic R n) hM k
+
+end rowStochastic
+
+namespace colStochastic
+
+variable [Semiring R] [PartialOrder R] [IsOrderedRing R]
+
+/-- Powers of column-stochastic matrices are column-stochastic. -/
+theorem pow_mem {M : Matrix n n R} (hM : M ∈ colStochastic R n) (k : ℕ) :
+    M ^ k ∈ colStochastic R n :=
+  Submonoid.pow_mem (colStochastic R n) hM k
+
+end colStochastic
+
+/-! ### Stationary distributions -/
+
+section Stationary
+
+variable {R : Type*} [Semiring R]
+
+/-- A distribution `μ` is stationary for a matrix `P` if `μ ᵥ* P = μ`. -/
+class IsStationary (μ : n → R) (P : Matrix n n R) : Prop where
+  stationary : μ ᵥ* P = μ
+
+/-- Powers of a matrix preserve stationary distributions. -/
+lemma IsStationary.pow (μ : n → R) (P : Matrix n n R) [IsStationary μ P] (k : ℕ) :
+    IsStationary μ (P ^ k) :=
+  ⟨by induction k with
+      | zero => simp [Matrix.vecMul_one]
+      | succ k ih => rw [pow_succ, ← Matrix.vecMul_vecMul, ih, IsStationary.stationary]⟩
+
+end Stationary
+
+/-! ### Cesàro averages and existence of stationary distributions -/
+
+section CesaroAverage
+
+variable {R : Type*} [DivisionRing R] [PartialOrder R] [IsStrictOrderedRing R]
+
+/-- The Cesàro average of the iterates of a vector under a matrix. -/
+def cesaroAverage (x₀ : n → R) (P : Matrix n n R) (k : ℕ) : n → R :=
+  (k + 1 : R)⁻¹ • ∑ i ∈ Finset.range (k + 1), x₀ ᵥ* (P ^ i)
+
+variable [Nonempty n]
+
+/-- The uniform distribution on a finite nonempty type. -/
+@[reducible, nolint unusedArguments]
+def uniformDistribution (R : Type*) (n : Type*) [Fintype n] [DivisionRing R] :
+    n → R := fun _ => 1 / Fintype.card n
+
+end CesaroAverage
+
+/-! ### Cesàro averages on the real simplex -/
+
+section CesaroAverageReal
+
+variable {M : Matrix n n ℝ}
+
+/-- The Cesàro average of a probability vector under a row-stochastic matrix
+belongs to the standard simplex. -/
+lemma cesaroAverage_mem_stdSimplex (hM : M ∈ rowStochastic ℝ n)
+    {x₀ : n → ℝ} (hx₀ : x₀ ∈ stdSimplex ℝ n) (k : ℕ) :
+    cesaroAverage x₀ M k ∈ stdSimplex ℝ n := by
+  have hmem i : x₀ ᵥ* M ^ i ∈ stdSimplex ℝ n :=
+    vecMul_mem_stdSimplex (rowStochastic.pow_mem hM i) hx₀
+  refine ⟨fun j => mul_nonneg (inv_nonneg.mpr (by linarith : (0 : ℝ) ≤ k + 1))
+    (by simp only [Finset.sum_apply]; exact sum_nonneg fun i _ => (hmem i).1 j), ?_⟩
+  simp only [cesaroAverage, Pi.smul_apply, smul_eq_mul, Finset.sum_apply, ← mul_sum,
+    sum_comm (γ := n)]
+  rw [show ∑ i ∈ range (k + 1), ∑ j, (x₀ ᵥ* M ^ i) j = k + 1 from
+    (sum_congr rfl fun i _ => (hmem i).2).trans (by simp), mul_comm]
+  exact mul_inv_cancel₀ (by linarith : (k : ℝ) + 1 ≠ 0)
+
+omit [DecidableEq n] in
+/-- The L¹ norm of a probability vector is 1. -/
+private lemma l1_nnnorm_eq_one {x : n → ℝ} (hx : x ∈ stdSimplex ℝ n) : ‖WithLp.toLp 1 x‖₊ = 1 := by
+  rw [← NNReal.coe_inj, NNReal.coe_one, l1_nnnorm_eq_sum,
+    (sum_congr rfl fun i _ => abs_of_nonneg (hx.1 i) : ∑ i, |x i| = ∑ i, x i), hx.2]
+
+/-- The Cesàro average is almost invariant: applying the matrix changes it by at most `2/(k+1)`. -/
+lemma norm_cesaroAverage_vecMul_sub_le (hM : M ∈ rowStochastic ℝ n)
+    {x₀ : n → ℝ} (hx₀ : x₀ ∈ stdSimplex ℝ n) (k : ℕ) :
+    ‖WithLp.toLp 1 (cesaroAverage x₀ M k ᵥ* M - cesaroAverage x₀ M k)‖ ≤ 2 / (k + 1) := by
+  have hk : (0 : ℝ) < k + 1 := by linarith
+  have hsimp : cesaroAverage x₀ M k ᵥ* M - cesaroAverage x₀ M k =
+      (k + 1 : ℝ)⁻¹ • (x₀ ᵥ* M ^ (k + 1) - x₀) := by
+    unfold cesaroAverage
+    rw [Matrix.smul_vecMul, ← smul_sub, Matrix.sum_vecMul]; congr 1
+    have h1 : ∀ i, (x₀ ᵥ* M ^ i) ᵥ* M = x₀ ᵥ* M ^ (i + 1) := fun i => by
+      rw [Matrix.vecMul_vecMul, ← pow_succ]
+    simp_rw [h1]
+    rw [Finset.sum_range_succ' (fun i => x₀ ᵥ* M ^ i), Finset.sum_range_succ, pow_zero,
+      Matrix.vecMul_one]; abel
+  rw [hsimp, WithLp.toLp_smul, norm_smul, Real.norm_eq_abs, abs_of_pos (inv_pos.mpr hk)]
+  have h1 := vecMul_mem_stdSimplex (rowStochastic.pow_mem hM (k + 1)) hx₀
+  have : ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1) - x₀)‖ ≤ 2 := calc
+    ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1) - x₀)‖₊
+    _ = ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1)) - WithLp.toLp 1 x₀‖₊ := by rw [← WithLp.toLp_sub]
+    _ ≤ ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1))‖₊ + ‖WithLp.toLp 1 x₀‖₊ := nnnorm_sub_le _ _
+    _ = 2 := by rw [l1_nnnorm_eq_one h1, l1_nnnorm_eq_one hx₀]; norm_num
+  calc (k + 1 : ℝ)⁻¹ * ‖WithLp.toLp 1 (x₀ ᵥ* M ^ (k + 1) - x₀)‖
+    _ ≤ (k + 1 : ℝ)⁻¹ * 2 := by gcongr
+    _ = 2 / (k + 1) := by ring
+
+variable [Nonempty n]
+
+omit [DecidableEq n] in
+/-- The uniform distribution belongs to the standard simplex. -/
+lemma uniformDistribution_mem_stdSimplex : uniformDistribution ℝ (n := n) ∈ stdSimplex ℝ n :=
+  ⟨fun _ => by simp only [uniformDistribution, one_div, inv_nonneg]; positivity,
+   by simp [uniformDistribution, Finset.card_univ, nsmul_eq_mul]⟩
+
+namespace rowStochastic
+
+/-- Every row-stochastic matrix on a finite nonempty state space has a stationary distribution. -/
+theorem exists_stationary_distribution (hM : M ∈ rowStochastic ℝ n) :
+    ∃ μ : n → ℝ, μ ∈ stdSimplex ℝ n ∧ IsStationary μ M := by
+  let x₀ := uniformDistribution ℝ (n := n)
+  let xₖ : ℕ → n → ℝ := fun k => cesaroAverage x₀ M k
+  have hxₖ k : xₖ k ∈ stdSimplex ℝ n :=
+    cesaroAverage_mem_stdSimplex hM uniformDistribution_mem_stdSimplex k
+  obtain ⟨μ, hμ_mem, nₖ, hnₖ_mono, hnₖ_lim⟩ := (isCompact_stdSimplex n).tendsto_subseq hxₖ
+  refine ⟨μ, hμ_mem, ⟨?_⟩⟩
+  have h_lim_diff : Tendsto (fun k => xₖ (nₖ k) ᵥ* M - xₖ (nₖ k)) atTop (nhds (μ ᵥ* M - μ)) :=
+    ((Continuous.matrix_vecMul continuous_id continuous_const).continuousAt.tendsto.comp
+      hnₖ_lim).sub hnₖ_lim
+  have h_error_tendsto : Tendsto (fun k => ‖WithLp.toLp 1 (xₖ (nₖ k) ᵥ* M - xₖ (nₖ k))‖)
+      atTop (nhds 0) :=
+    tendsto_of_tendsto_of_tendsto_of_le_of_le tendsto_const_nhds
+      ((tendsto_const_nhds.div_atTop tendsto_id).comp
+        ((tendsto_natCast_atTop_atTop.comp hnₖ_mono.tendsto_atTop).atTop_add tendsto_const_nhds))
+      (fun _ => norm_nonneg _)
+      (fun k => norm_cesaroAverage_vecMul_sub_le hM uniformDistribution_mem_stdSimplex (nₖ k))
+  have h_norm_zero := tendsto_nhds_unique
+    (((PiLp.continuous_toLp 1 (fun _ => ℝ)).continuousAt.tendsto.comp h_lim_diff).norm)
+    h_error_tendsto
+  exact sub_eq_zero.mp ((WithLp.toLp_injective 1).eq_iff.mp (norm_eq_zero.mp h_norm_zero))
+
+end rowStochastic
+
+end CesaroAverageReal