Direct sums, Kronecker products and the block structure of product affine systems

Definition

SASAlgebra

by olivier · Sep 11, 2026 · Mathlib 0df444a (Lean v4.33.1)

kronecker-productlinear-algebrareservoir-computingstate-affine-system

This module fixes the algebraic apparatus behind the closure properties of state-affine reservoir systems: direct sums and Kronecker products of polynomial coefficient families, and the block structure of the product system.

Polynomials with matrix coefficients. A family of coefficients is evaluated as $p(z) = \sum_j z^j P_j$ . The direct sum of two such families is taken coefficient by coefficient, giving $\left(p_1 \oplus p_2\right)(z) = p_1(z) \oplus p_2(z)$ ; the Kronecker product is indexed by pairs of coefficient indices, giving $\sum_{i,j} z^{i+j} \left(A^1_i \otimes A^2_j\right) = p_1(z) \otimes p_2(z)$ , which is equation (3.20) of the source.

The product system. Given two affine recursions $x_1 = p_1 x_1' + q_1$ and $x_2 = p_2 x_2' + q_2$ , the tensor product $x_1 \otimes x_2$ obeys

x_1 \otimes x_2 \;=\; (p_1\otimes p_2)(x_1'\otimes x_2') \;+\; (p_1\otimes q_2)\,x_1' \;+\; (q_1\otimes p_2)\,x_2' \;+\; q_1\otimes q_2 .

The two cross terms act on $x_1'$ and $x_2'$ separately, and the last is constant. An augmented state carrying $x_1$ , $x_2$ and $x_1 \otimes x_2$ is therefore needed, with a block-triangular matrix; the constant term $q_1 \otimes q_2$ is what makes the non-homogeneous form necessary, and is why the corresponding closure fails for purely linear reservoirs.

Formalization Note The augmented state is indexed by the disjoint union $(\iota_1 \oplus \iota_2) \oplus (\iota_1 \times \iota_2)$ ; the block picture is a rendering of that index type and fixes no ordering or basis. The two cross blocks are given as explicit matrices rather than through a column-replication construction. Coefficient families are indexed by Fin r, so both polynomials carry the same number of coefficients; padding by zero is how a shorter one is accommodated. The recursion runs from the larger index to the smaller, so that index $k$ denotes the instant $k$ steps into the past. Degenerate dimensions are admitted and make the statements vacuous rather than false.

Definition code

import Mathlib

set_option autoImplicit false
open Matrix

namespace SASAlgebra

variable {N₁ N₂ r : ℕ}

/-- Evaluation d'un polynome a coefficients matriciels, indexee par un type quelconque. -/
noncomputable def polyEval {ι : Type*} [Fintype ι] [DecidableEq ι]
    (P : Fin r → Matrix ι ι ℝ) (z : ℝ) : Matrix ι ι ℝ :=
  ∑ j : Fin r, (z ^ (j : ℕ)) • P j

/-- Somme directe de deux polynomes matriciels, coefficient par coefficient. -/
noncomputable def polyAdd (P₁ : Fin r → Matrix (Fin N₁) (Fin N₁) ℝ)
    (P₂ : Fin r → Matrix (Fin N₂) (Fin N₂) ℝ) :
    Fin r → Matrix (Fin N₁ ⊕ Fin N₂) (Fin N₁ ⊕ Fin N₂) ℝ :=
  fun j => Matrix.fromBlocks (P₁ j) 0 0 (P₂ j)

/-- L'evaluation commute avec la somme directe : `(p₁ ⊕ p₂)(z) = p₁(z) ⊕ p₂(z)`. -/
lemma polyEval_polyAdd (P₁ : Fin r → Matrix (Fin N₁) (Fin N₁) ℝ)
    (P₂ : Fin r → Matrix (Fin N₂) (Fin N₂) ℝ) (z : ℝ) :
    polyEval (polyAdd P₁ P₂) z
      = Matrix.fromBlocks (polyEval P₁ z) 0 0 (polyEval P₂ z) := by
  simp only [polyEval, polyAdd]
  ext i j
  rcases i with i | i <;> rcases j with j | j <;>
    simp [Matrix.fromBlocks, Matrix.sum_apply, Finset.sum_apply]

/-- Somme directe de deux polynomes a coefficients vectoriels. -/
noncomputable def vpolyAdd (Q₁ : Fin r → (Fin N₁ → ℝ)) (Q₂ : Fin r → (Fin N₂ → ℝ)) :
    Fin r → (Fin N₁ ⊕ Fin N₂ → ℝ) :=
  fun j => Sum.elim (Q₁ j) (Q₂ j)

noncomputable def vpolyEval {ι : Type*} (Q : Fin r → (ι → ℝ)) (z : ℝ) : ι → ℝ :=
  ∑ j : Fin r, (z ^ (j : ℕ)) • Q j

lemma vpolyEval_vpolyAdd (Q₁ : Fin r → (Fin N₁ → ℝ)) (Q₂ : Fin r → (Fin N₂ → ℝ)) (z : ℝ) :
    vpolyEval (vpolyAdd Q₁ Q₂) z = Sum.elim (vpolyEval Q₁ z) (vpolyEval Q₂ z) := by
  funext i
  rcases i with i | i <;>
    simp [vpolyEval, vpolyAdd, Finset.sum_apply]

/-- Le systeme d'etat affine, pour un type d'index quelconque. -/
def IsSAS {ι : Type*} [Fintype ι] [DecidableEq ι]
    (P : Fin r → Matrix ι ι ℝ) (Q : Fin r → (ι → ℝ))
    (z : ℕ → ℝ) (x : ℕ → (ι → ℝ)) : Prop :=
  ∀ k, x k = polyEval P (z k) *ᵥ x (k + 1) + vpolyEval Q (z k)

/-- **Proposition 3.10(i), moitie « etat ».** Le couple de deux solutions resout le
systeme somme directe : c'est ce qui fait des combinaisons lineaires de fonctionnelles
SAS des fonctionnelles SAS. -/
theorem isSAS_polyAdd
    (P₁ : Fin r → Matrix (Fin N₁) (Fin N₁) ℝ) (Q₁ : Fin r → (Fin N₁ → ℝ))
    (P₂ : Fin r → Matrix (Fin N₂) (Fin N₂) ℝ) (Q₂ : Fin r → (Fin N₂ → ℝ))
    (z : ℕ → ℝ) (x₁ : ℕ → (Fin N₁ → ℝ)) (x₂ : ℕ → (Fin N₂ → ℝ))
    (h₁ : IsSAS P₁ Q₁ z x₁) (h₂ : IsSAS P₂ Q₂ z x₂) :
    IsSAS (polyAdd P₁ P₂) (vpolyAdd Q₁ Q₂) z (fun k => Sum.elim (x₁ k) (x₂ k)) := by
  intro k
  rw [polyEval_polyAdd, vpolyEval_vpolyAdd, Matrix.fromBlocks_mulVec]
  funext i
  rcases i with i | i
  · have := congrFun (h₁ k) i
    simpa using this
  · have := congrFun (h₂ k) i
    simpa using this

/-- **Proposition 3.10(i), moitie « lecture ».** Le readout `W₁ ⊕ λW₂` sur l'etat couple
rend la combinaison lineaire des deux sorties. -/
theorem readout_polyAdd (W₁ : Fin N₁ → ℝ) (W₂ : Fin N₂ → ℝ) (lam : ℝ)
    (u₁ : Fin N₁ → ℝ) (u₂ : Fin N₂ → ℝ) :
    Sum.elim W₁ (lam • W₂) ⬝ᵥ Sum.elim u₁ u₂ = W₁ ⬝ᵥ u₁ + lam * (W₂ ⬝ᵥ u₂) := by
  simp only [dotProduct, Fintype.sum_sum_type, Sum.elim_inl, Sum.elim_inr, Pi.smul_apply,
    smul_eq_mul]
  rw [Finset.mul_sum]
  congr 1
  exact Finset.sum_congr rfl (fun i _ => by ring)

/-! ### Partie (ii) : le produit

Le produit de deux fonctionnelles SAS demande de porter, outre les deux etats, leur
produit tensoriel : si `x₁` et `x₂` resolvent leurs systemes, alors `x₁ ⊗ x₂` resout un
systeme dont la matrice est `p₁(z) ⊗ p₂(z)` et le terme affine fait intervenir
`q₁ ⊗ q₂`. C'est ce qui force la forme non homogene. -/

/-- Produit de Kronecker de deux polynomes matriciels, a degres egaux :
`(p₁ ⊗ p₂)(z) = p₁(z) ⊗ p₂(z)` n'est pas un polynome de meme degre, d'ou l'indexation
par les couples d'indices de coefficients. -/
noncomputable def polyKron (P₁ : Fin r → Matrix (Fin N₁) (Fin N₁) ℝ)
    (P₂ : Fin r → Matrix (Fin N₂) (Fin N₂) ℝ) :
    Fin r × Fin r → Matrix (Fin N₁ × Fin N₂) (Fin N₁ × Fin N₂) ℝ :=
  fun ij => Matrix.kroneckerMap (· * ·) (P₁ ij.1) (P₂ ij.2)

/-- L'evaluation du produit tensoriel : `Σ_{i,j} z^{i+j} (A¹ᵢ ⊗ A²ⱼ) = p₁(z) ⊗ p₂(z)`.
C'est l'equation (3.20) de la source. -/
lemma kron_eval (P₁ : Fin r → Matrix (Fin N₁) (Fin N₁) ℝ)
    (P₂ : Fin r → Matrix (Fin N₂) (Fin N₂) ℝ) (z : ℝ) :
    ∑ ij : Fin r × Fin r, (z ^ ((ij.1 : ℕ) + (ij.2 : ℕ))) • polyKron P₁ P₂ ij
      = Matrix.kroneckerMap (· * ·) (polyEval P₁ z) (polyEval P₂ z) := by
  ext a b
  simp only [polyKron, polyEval, Matrix.kroneckerMap_apply, Matrix.sum_apply,
    Finset.sum_apply, Matrix.smul_apply, smul_eq_mul, Fintype.sum_prod_type]
  rw [Finset.sum_mul_sum]
  exact Finset.sum_congr rfl (fun i _ => Finset.sum_congr rfl (fun j _ => by
    rw [pow_add]; ring))

/-- Produit tensoriel de deux vecteurs. -/
def vkron {ι κ : Type*} (u : ι → ℝ) (v : κ → ℝ) : ι × κ → ℝ := fun ij => u ij.1 * v ij.2

/-- `(A u) ⊗ (B v) = (A ⊗ B) (u ⊗ v)` : le produit de Kronecker transporte l'action. -/
lemma vkron_mulVec {m n p q : Type*} [Fintype n] [Fintype q] [DecidableEq n] [DecidableEq q]
    (A : Matrix m n ℝ) (B : Matrix p q ℝ) (u : n → ℝ) (v : q → ℝ) :
    vkron (A *ᵥ u) (B *ᵥ v) = (Matrix.kroneckerMap (· * ·) A B) *ᵥ (vkron u v) := by
  funext ij
  simp only [vkron, Matrix.mulVec, dotProduct, Matrix.kroneckerMap_apply,
    Fintype.sum_prod_type]
  rw [Finset.sum_mul_sum]
  exact Finset.sum_congr rfl (fun i _ => Finset.sum_congr rfl (fun j _ => by ring))

/-- Les deux matrices croisees : `p₁ ⊗ q₂` et `q₁ ⊗ p₂`, vues comme applications
lineaires de `x₁` (resp. `x₂`) vers l'espace produit. -/
def kronRight {ι κ ν : Type*} (A : Matrix ι ν ℝ) (v : κ → ℝ) : Matrix (ι × κ) ν ℝ :=
  fun ij a => A ij.1 a * v ij.2

def kronLeft {ι κ ν : Type*} (u : ι → ℝ) (B : Matrix κ ν ℝ) : Matrix (ι × κ) ν ℝ :=
  fun ij a => u ij.1 * B ij.2 a

/-- **Le decoupage du produit.** Si `x₁` et `x₂` resolvent leurs systemes, le produit
tensoriel verifie une recurrence affine dont la partie lineaire porte sur `x₁`, `x₂` et
`x₁ ⊗ x₂` simultanement. C'est pourquoi l'etat du systeme produit doit porter les trois,
et pourquoi le terme affine `q₁ ⊗ q₂` apparait : la forme homogene ne suffit pas. -/
theorem vkron_expand
    (p₁ : Matrix (Fin N₁) (Fin N₁) ℝ) (q₁ : Fin N₁ → ℝ)
    (p₂ : Matrix (Fin N₂) (Fin N₂) ℝ) (q₂ : Fin N₂ → ℝ)
    (u₁ : Fin N₁ → ℝ) (u₂ : Fin N₂ → ℝ) :
    vkron (p₁ *ᵥ u₁ + q₁) (p₂ *ᵥ u₂ + q₂)
      = (Matrix.kroneckerMap (· * ·) p₁ p₂) *ᵥ (vkron u₁ u₂)
        + (kronRight p₁ q₂) *ᵥ u₁
        + (kronLeft q₁ p₂) *ᵥ u₂
        + vkron q₁ q₂ := by
  have hk := vkron_mulVec p₁ p₂ u₁ u₂
  funext ij
  have hk' := congrFun hk ij
  show ((p₁ *ᵥ u₁) ij.1 + q₁ ij.1) * ((p₂ *ᵥ u₂) ij.2 + q₂ ij.2) = _
  rw [Pi.add_apply, Pi.add_apply, Pi.add_apply, ← hk']
  show _ = vkron (p₁ *ᵥ u₁) (p₂ *ᵥ u₂) ij + _ + _ + _
  simp only [vkron, kronRight, kronLeft, Matrix.mulVec, dotProduct]
  have e1 : ∑ x, p₁ ij.1 x * q₂ ij.2 * u₁ x = (∑ x, p₁ ij.1 x * u₁ x) * q₂ ij.2 := by
    rw [Finset.sum_mul]
    exact Finset.sum_congr rfl (fun x _ => by ring)
  have e2 : ∑ x, q₁ ij.1 * p₂ ij.2 x * u₂ x = q₁ ij.1 * (∑ x, p₂ ij.2 x * u₂ x) := by
    rw [Finset.mul_sum]
    exact Finset.sum_congr rfl (fun x _ => by ring)
  rw [e1, e2]
  ring

/-! ### Le systeme produit

L'etat du systeme produit porte les trois composantes `x₁`, `x₂` et `x₁ ⊗ x₂`, indexees
par `(ι₁ ⊕ ι₂) ⊕ (ι₁ × ι₂)`. Sa matrice est triangulaire par blocs : les deux premieres
composantes evoluent librement, la troisieme recoit les deux termes croises. -/

/-- La ligne croisee : de `x₁` et `x₂` vers la composante tensorielle. -/
noncomputable def crossRow (p₁ : Matrix (Fin N₁) (Fin N₁) ℝ) (q₁ : Fin N₁ → ℝ)
    (p₂ : Matrix (Fin N₂) (Fin N₂) ℝ) (q₂ : Fin N₂ → ℝ) :
    Matrix (Fin N₁ × Fin N₂) (Fin N₁ ⊕ Fin N₂) ℝ :=
  Matrix.of fun ij a => Sum.elim (fun a₁ => kronRight p₁ q₂ ij a₁)
    (fun a₂ => kronLeft q₁ p₂ ij a₂) a

/-- La matrice du systeme produit, par blocs. -/
noncomputable def prodMat (p₁ : Matrix (Fin N₁) (Fin N₁) ℝ) (q₁ : Fin N₁ → ℝ)
    (p₂ : Matrix (Fin N₂) (Fin N₂) ℝ) (q₂ : Fin N₂ → ℝ) :
    Matrix ((Fin N₁ ⊕ Fin N₂) ⊕ (Fin N₁ × Fin N₂))
           ((Fin N₁ ⊕ Fin N₂) ⊕ (Fin N₁ × Fin N₂)) ℝ :=
  Matrix.fromBlocks (Matrix.fromBlocks p₁ 0 0 p₂) 0
    (crossRow p₁ q₁ p₂ q₂) (Matrix.kroneckerMap (· * ·) p₁ p₂)

/-- Le terme affine du systeme produit. -/
def prodVec (q₁ : Fin N₁ → ℝ) (q₂ : Fin N₂ → ℝ) :
    (Fin N₁ ⊕ Fin N₂) ⊕ (Fin N₁ × Fin N₂) → ℝ :=
  Sum.elim (Sum.elim q₁ q₂) (vkron q₁ q₂)

/-- L'etat augmente : les deux etats et leur produit tensoriel. -/
def prodState (u₁ : Fin N₁ → ℝ) (u₂ : Fin N₂ → ℝ) :
    (Fin N₁ ⊕ Fin N₂) ⊕ (Fin N₁ × Fin N₂) → ℝ :=
  Sum.elim (Sum.elim u₁ u₂) (vkron u₁ u₂)

/-- **Proposition 3.10(ii), le pas central.** L'etat augmente verifie la recurrence
affine du systeme produit. Les deux termes croises de `vkron_expand` sont exactement
les deux blocs de `crossRow`, et `q₁ ⊗ q₂` le bloc affine correspondant. -/
theorem prodState_step
    (p₁ : Matrix (Fin N₁) (Fin N₁) ℝ) (q₁ : Fin N₁ → ℝ)
    (p₂ : Matrix (Fin N₂) (Fin N₂) ℝ) (q₂ : Fin N₂ → ℝ)
    (u₁ : Fin N₁ → ℝ) (u₂ : Fin N₂ → ℝ) :
    prodState (p₁ *ᵥ u₁ + q₁) (p₂ *ᵥ u₂ + q₂)
      = prodMat p₁ q₁ p₂ q₂ *ᵥ prodState u₁ u₂ + prodVec q₁ q₂ := by
  funext i
  rcases i with i | i
  · -- les deux premieres composantes : evolution libre
    rcases i with i | i <;>
      simp [prodState, prodMat, prodVec, Matrix.fromBlocks_mulVec, Matrix.mulVec, dotProduct,
        Fintype.sum_sum_type, Fintype.sum_prod_type]
  · -- la composante tensorielle : c'est `vkron_expand`
    have hx := congrFun (vkron_expand p₁ q₁ p₂ q₂ u₁ u₂) i
    simp only [prodState, prodMat, prodVec, Sum.elim_inr, Matrix.fromBlocks_mulVec,
      Pi.add_apply] at *
    rw [hx]
    simp only [crossRow, Matrix.mulVec, dotProduct, Matrix.of_apply, Fintype.sum_sum_type,
      Sum.elim_inl, Sum.elim_inr, Function.comp_apply]
    ring

/-- Le readout du systeme produit : nul sur les deux premieres composantes,
`W₁ ⊗ W₂` sur la composante tensorielle. -/
def prodReadout (W₁ : Fin N₁ → ℝ) (W₂ : Fin N₂ → ℝ) :
    (Fin N₁ ⊕ Fin N₂) ⊕ (Fin N₁ × Fin N₂) → ℝ :=
  Sum.elim (Sum.elim 0 0) (vkron W₁ W₂)

/-- **Proposition 3.10(ii), moitie « lecture ».** Le readout `0 ⊕ 0 ⊕ (W₁ ⊗ W₂)` sur
l'etat augmente rend le produit des deux sorties. -/
theorem readout_prod (W₁ : Fin N₁ → ℝ) (W₂ : Fin N₂ → ℝ)
    (u₁ : Fin N₁ → ℝ) (u₂ : Fin N₂ → ℝ) :
    prodReadout W₁ W₂ ⬝ᵥ prodState u₁ u₂ = (W₁ ⬝ᵥ u₁) * (W₂ ⬝ᵥ u₂) := by
  simp only [prodReadout, prodState, dotProduct, Fintype.sum_sum_type, Sum.elim_inl,
    Sum.elim_inr, Pi.zero_apply, zero_mul, Finset.sum_const_zero, vkron,
    Fintype.sum_prod_type, zero_add]
  rw [Finset.sum_mul_sum]
  exact Finset.sum_congr rfl (fun i _ => Finset.sum_congr rfl (fun j _ => by ring))


end SASAlgebra

Source

L. Grigoryeva, J.-P. Ortega, Universal discrete-time reservoir computers with stochastic inputs and linear readouts using non-homogeneous state-affine systems, Journal of Machine Learning Research 19(24) (2018), 1-40, https://arxiv.org/abs/1712.00754, p. 12, Proposition 3.10 and equations (3.19)-(3.22).

View graph