Résolution numérique de problèmes aux valeurs propres

Cours ENPC — Pratique du calcul scientifique

Objectif

But de ce chapitre

Étudier les méthodes numériques pour trouver les couples propres $(\mathbf{\boldsymbol{v}}_i, \mathbf{\boldsymbol{\lambda}}_i)\in \mathbf C^n \times \mathbf C$ vérifiant \[ \mathsf A \mathbf{\boldsymbol{v}}_i= \lambda_i \mathbf{\boldsymbol{v}}_i, \] où $\mathsf A \in \mathbf C^{n \times n}$.

$~$

Hypothèse simplificatrice tout au long de ce chapitre : la matrice $\mathsf A$ est diagonalisable.

Notation : valeurs propres $|\lambda_1| \geqslant\dotsc \geqslant|\lambda_n|$, vecteurs propres normalisés $\mathbf{\boldsymbol{v}}_1, \dotsc, \mathbf{\boldsymbol{v}}_n$. \[ \mathsf A \mathsf V = \mathsf V \mathsf D, \qquad \mathsf V = \begin{pmatrix} \mathbf{\boldsymbol{v}}_1 & \dots & \mathbf{\boldsymbol{v}}_n \end{pmatrix}, \qquad \mathsf D = {\rm diag}(\lambda_1, \dotsc, \lambda_n). \]

Calculer toutes les valeurs propres d’une matrice est coûteux, même avec les meilleures méthodes :

using LinearAlgebra
A = rand(3000, 3000)
@elapsed(eigen(A))

8.625416 seconds

${\color{green} \rightarrow}$ On s’intéresse d’abord aux méthodes pour trouver un ou quelques couples propres seulement.

Programme du jour

Itérations vectorielles simples

calcule un seul couple propre
Itération simultanée

calcule plusieurs couples propres simultanément
Application : PageRank

trouver la distribution invariante d’une chaîne de Markov

Remarque : peut-on calculer les racines du polynôme caractéristique ?

Le coût de calcul du polynôme caractéristique est en ${\color{red}n!}$

${\color{green} \rightarrow}$ pas une approche viable…

La méthode des puissances

Méthode des puissances

La méthode des puissances est une méthode pour trouver $\color{blue}(\lambda_1, v_1)$ : \[ \fcolorbox{blue}{lightgreen}{$ \mathbf{\boldsymbol{x}}_{k+1} = \frac{\mathsf A \mathbf{\boldsymbol{x}}_{k}}{{\lVert {\mathsf A \mathbf{\boldsymbol{x}}_{k}} \rVert}} $} \] Étant donné $\mathbf{\boldsymbol{x}}_k \approx \mathbf{\boldsymbol{v}}_1$, la valeur propre est calculée à partir du quotient de Rayleigh : \[ \fcolorbox{blue}{lightgreen}{$ \lambda_{1} \approx \frac{\mathbf{\boldsymbol{x}}_{k}^* \mathsf A \mathbf{\boldsymbol{x}}_k}{\mathbf{\boldsymbol{x}}_k^* \mathbf{\boldsymbol{x}}_k} $} \]

La normalisation est nécessaire pour éviter les Inf numériques.

A = [2 1; 1 2]
x = rand(2)
for i in 1:20
    x = A*x
    x /= √(x'x)
end
println("Vecteur propre : $x, valeur propre : $(x'A*x)")

Vecteur propre : [0.7071067810178944, 0.7071067813552007], valeur propre : 2.9999999999999996

Convergence de la méthode des puissances

Notation : notons $\angle$ l’angle aigu entre $\mathbf{\boldsymbol{x}}$ et $\mathbf{\boldsymbol{y}}$ : \[ \angle(\mathbf{\boldsymbol{x}}, \mathbf{\boldsymbol{y}}) = \arccos\left( \frac{{\lvert {\mathbf{\boldsymbol{x}}^* \mathbf{\boldsymbol{y}}} \rvert}}{\sqrt{\mathbf{\boldsymbol{x}}^* \mathbf{\boldsymbol{x}}} \sqrt{\mathbf{\boldsymbol{y}}^* \mathbf{\boldsymbol{y}}}}\right). \]

Puisque les vecteurs propres de $\mathsf A$ forment une base de $\mathbf C^n$, tout vecteur $\mathbf{\boldsymbol{x}}_0$ peut se décomposer sous la forme \[ \mathbf{\boldsymbol{x}}_0 = \alpha_1 \mathbf{\boldsymbol{v}}_1 + \dotsb + \alpha_n \mathbf{\boldsymbol{v}}_n. \]

Proposition : convergence de la méthode des puissances démarrée depuis $\mathbf{\boldsymbol{x}}_0$

Supposons que ${\lvert {\lambda_1} \rvert} > {\lvert {\lambda_2} \rvert}$ et $\alpha_1 \neq 0$. Alors $(\mathbf{\boldsymbol{x}}_k)_{k \geqslant 0}$ vérifie \[ \lim_{k \to \infty} \angle(\mathbf{\boldsymbol{x}}_k, \mathbf{\boldsymbol{v}}_1) = 0. \] La suite $(\mathbf{\boldsymbol{x}}_k)$ est dite converger essentiellement vers $\mathbf{\boldsymbol{v}}_1$. $~$

Preuve

Par construction \[ \mathbf{\boldsymbol{x}}_k = \frac{\lambda_1^k \alpha_1 \mathbf{\boldsymbol{v}}_1 + \dotsb + \lambda_n^k \alpha_n \mathbf{\boldsymbol{v}}_n}{{\lVert {\lambda_1^k \alpha_1 \mathbf{\boldsymbol{v}}_1 + \dotsb + \lambda_n^k \alpha_n \mathbf{\boldsymbol{v}}_n} \rVert}} = \frac{\lambda_1^k \alpha_1}{|\lambda_1^k \alpha_1|} \frac{\mathbf{\boldsymbol{v}}_1 + \frac{\lambda_2^k \alpha_2}{\lambda_1^k \alpha_1} \mathbf{\boldsymbol{v}}_2 + \dotsb + \frac{\lambda_n^k \alpha_2}{\lambda_1^k \alpha_1} \mathbf{\boldsymbol{v}}_n }{\left\lVert\mathbf{\boldsymbol{v}}_1 + \frac{\lambda_2^k \alpha_2}{\lambda_1^k \alpha_1} \mathbf{\boldsymbol{v}}_2 + \dotsb + \frac{\lambda_n^k \alpha_n}{\lambda_1^k \alpha_1} \mathbf{\boldsymbol{v}}_n\right\rVert}. \] Donc \[ \mathbf{\boldsymbol{x}}_k \left(\frac{|\lambda_1^k \alpha_1|}{\lambda_1^k \alpha_1}\right) \xrightarrow[k \to \infty]{} \mathbf{\boldsymbol{v}}_1. \] (Remarquer que le facteur entre parenthèses dans le membre de gauche est un scalaire de module 1.)

Le terme d’erreur dominant est en $\color{red} |\lambda_2/\lambda_1|^k$.

Itération inverse et quotient de Rayleigh

Itération inverse : pour trouver $\lambda_J$, la valeur propre la plus proche de $\mu$

Soit $\mathsf B = (\mathsf A - \mu \mathsf I)^{-1}$. L’itération inverse autour de $\mu \in \mathbf C$ est donnée par \[ \fcolorbox{blue}{lightgreen}{$ \mathbf{\boldsymbol{x}}_{k+1} = \frac{\mathsf B \mathbf{\boldsymbol{x}}_{k}}{{\lVert {\mathsf B \mathbf{\boldsymbol{x}}_{k}} \rVert}}. $} \] Étant donné $\mathbf{\boldsymbol{x}}_k \approx \mathbf{\boldsymbol{v}}_J$, la valeur propre est calculée à partir du quotient de Rayleigh : \[ \fcolorbox{blue}{lightgreen}{$ \lambda_{J} \approx \frac{\mathbf{\boldsymbol{x}}_{k}^* \mathsf A \mathbf{\boldsymbol{x}}_k}{\mathbf{\boldsymbol{x}}_k^* \mathbf{\boldsymbol{x}}_k} $} \]

La valeur propre dominante de $\mathsf B$ est donnée par $(\lambda_J - \mu)^{-1}$.
Le vecteur propre associé est $\mathbf{\boldsymbol{v}}_J$.
Utiliser l’opérateur \ plutôt que de calculer l’inverse :

A = [2 1; 1 2]
x = rand(2)
μ = 0
B = factorize(A - μ*I)
for i in 1:10
    x = B\x
    x /= √(x'x)
end
println("Vecteur propre : $x, valeur propre : $(x'A*x)")

Vecteur propre : [-0.7067802249705825, 0.7074331866618451], valeur propre : 1.0000004263589704

Peut-on adapter $\mu$ au cours de l’itération pour accélérer la convergence ?

L’erreur de l’itération inverse est en $\color{red} \left|\frac{\lambda_J - \mu} {\lambda_K - \mu}\right|^k$.

${\color{green} \rightarrow}$ suggère de prendre \[ \fcolorbox{blue}{lightgreen}{$ \mu_k = \frac{\mathbf{\boldsymbol{x}}_k^* \mathsf A \mathbf{\boldsymbol{x}}_k}{\mathbf{\boldsymbol{x}}_k^* \mathbf{\boldsymbol{x}}_k}. $} \]

Cela conduit à l’itération du quotient de Rayleigh.

Si la convergence a lieu, la valeur propre associée converge cubiquement :

\[ |\lambda^{(k+1)} - \lambda_J| \leqslant C |\lambda^{(k)} - \lambda_J|^3. \]

${\color{green} \rightarrow}$ Le nombre de chiffres significatifs est environ triplé à chaque itération.

A = [2 1; 1 2]
setprecision(500)
x = BigFloat.([1.2; 0.9])
μ = 0
for i in 1:6
    x = (A - μ*I)\x
    x /= √(x'x)
    μ = x'A*x
    println(μ)
end

2.689655172413793214338566074331922284628319727818488081860106227371262853189736760002654726539617400119969038165484095334893321873413803868251681795864
2.9876835222760986149331492272467355588185696109799282648595173771965397075426195063225678761434065353075142377089820027363840168362403721691362285304077
2.9999995241744513496497052218657690803207923292634754813048708145726507175594474194389978335452546331521437999768268695077710762632721314300139734193988
2.9999999999999999999730670707803381925912042248826047017208189162200256216460941091027939474791840614585910098581290345505474157692691833334644060848503
2.9999999999999999999999999999999999999999999999999999999999951158299301653110614230820165226784762374074665384291286744975203198864866424226623297764997
3.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000012

Critère d’arrêt (pour information seulement)

Un bon critère d’arrêt est \[ \begin{equation} \tag{critère d'arrêt} {\lVert {\mathsf A \widehat {\mathbf{\boldsymbol{v}}} - \widehat \lambda \widehat {\mathbf{\boldsymbol{v}}}} \rVert} \leqslant\varepsilon {\lVert {\widehat {\mathbf{\boldsymbol{v}}}} \rVert} \label{stopping} \end{equation} \]

Proposition : estimation d’erreur directe

Supposons $\mathsf A$ hermitienne et $\eqref{stopping}$ vérifié. Alors il existe $\lambda \in \sigma(\mathsf A)$ tel que $|\lambda - \widehat \lambda| \leqslant\varepsilon$.

Preuve

Soit $\mathbf{\boldsymbol{r}} := \mathsf A \widehat {\mathbf{\boldsymbol{v}}} - \widehat \lambda \widehat {\mathbf{\boldsymbol{v}}}$. Alors \[ \frac{1}{\min_{\lambda \in \sigma(\mathsf A)} |\lambda - \widehat \lambda|} = {\lVert {(\mathsf A - \widehat \lambda \mathsf I)^{-1}} \rVert} \geqslant\frac{{\lVert {(\mathsf A - \widehat \lambda \mathsf I)^{-1}\mathbf{\boldsymbol{r}}} \rVert}}{{\lVert {\mathbf{\boldsymbol{r}}} \rVert}} \geqslant\frac{1}{\varepsilon}, \] En prenant le réciproque, on déduit que \[ \min_{\lambda \in \sigma(\mathsf A)} |\lambda - \widehat \lambda| \leqslant\varepsilon, \] ce qui conclut la preuve.

Proposition : estimation d’erreur rétrograde

Supposons que $\eqref{stopping}$ soit vérifié et soit $\mathcal E = \Bigl\{\mathsf E \in \mathbf C^{n \times n}: (\mathsf A + \mathsf E) \widehat {\mathbf{\boldsymbol{v}}} = \widehat \lambda \widehat {\mathbf{\boldsymbol{v}}} \Bigr\}$. Alors \[ \min_{\mathsf E \in \mathcal E} {\lVert {\mathsf E} \rVert} = \frac{{\lVert {\mathsf A \widehat {\mathbf{\boldsymbol{v}}} - \widehat \lambda \widehat {\mathbf{\boldsymbol{v}}}} \rVert}}{{\lVert {\widehat {\mathbf{\boldsymbol{v}}}} \rVert}} \quad(\leqslant\varepsilon). \]

Définition de l’erreur rétrograde par l’analyste numérique Nick Higham :

L’erreur rétrograde est une mesure de l’erreur associée à une solution approchée d’un problème. Tandis que l’erreur directe est la distance entre la solution approchée et la solution exacte, l’erreur rétrograde mesure de combien les données doivent être perturbées pour produire la solution approchée.

Esquisse de preuve

Montrer d’abord que tout $\mathsf E \in \mathcal E$ vérifie $\mathsf E \widehat {\mathbf{\boldsymbol{v}}} = - \mathbf{\boldsymbol{r}}$, où $\mathbf{\boldsymbol{r}} := \mathsf A \widehat {\mathbf{\boldsymbol{v}}} - \widehat \lambda \widehat {\mathbf{\boldsymbol{v}}}$.
Déduire du premier point que \[ \inf_{\mathsf E \in \mathcal E} {\lVert {\mathsf E} \rVert} \geqslant\frac{{\lVert {\mathbf{\boldsymbol{r}}} \rVert}}{{\lVert {\widehat {\mathbf{\boldsymbol{v}}}} \rVert}}. \]
Trouver une matrice de rang 1 $\mathsf E_*$ telle que ${\lVert {\mathsf E_*} \rVert} = \frac{{\lVert {\mathbf{\boldsymbol{r}}} \rVert}}{{\lVert {\widehat {\mathbf{\boldsymbol{v}}}} \rVert}}$, puis conclure.

(Rappeler que toute matrice de rang 1 est de la forme $\mathsf E_* = \mathbf{\boldsymbol{u}} \mathbf{\boldsymbol{w}}^*$, de norme ${\lVert {\mathbf{\boldsymbol{u}}} \rVert} {\lVert {\mathbf{\boldsymbol{w}}} \rVert}$.)

Calcul simultané de plusieurs couples propres

Préliminaire : décomposition QR

Toute matrice $\mathsf A \in \mathbf C^{n \times p}$ (avec $p ⩽ n$) de rang $p$ peut se décomposer sous la forme $\mathsf A = \mathsf Q \mathsf R$ où $\mathsf Q\in \mathbf C^{n \times p}$ est une matrice orthogonale (${\mathsf Q}^* \mathsf Q = \mathsf I_{p \times p}$) et $\mathsf R\in \mathbf C^{p \times p}$ est une matrice triangulaire supérieure.

Remarques sur les matrices orthogonales rectangulaires :

Si ${(\mathbf q_j)}_{1⩽j⩽p} \in \mathbf C^{n}$ désignent les $p$ vecteurs colonnes orthonormaux de $\mathsf Q$ et ${(\mathbf e_j)}_{1⩽j⩽p} \in \mathbf C^{p}$ les $p$ vecteurs canoniques de $\mathbf C^{p}$, alors $\mathsf Q = \sum_{j=1}^p {\mathbf q}_j {\mathbf e}_j^*$
Il s’ensuit que ${\mathsf Q}^* \mathsf Q = \sum_{j=1}^p {\mathbf e}_j {\mathbf e}_j^* = \mathsf I_{p \times p}$ mais $\mathsf Q {\mathsf Q}^* = \sum_{j=1}^p {\mathbf q}_j {\mathbf q}_j^* \in \mathbf C^{n \times n}$ est le projecteur orthogonal sur les colonnes de $\mathsf Q$ et est différent de $\mathsf I_{n \times n}$ dès que $p<n$.

Preuve par récurrence sur $p$

$p=1$: $\quad \mathsf A = \begin{pmatrix}\mathbf a \end{pmatrix} = \begin{pmatrix}\frac{\mathbf a}{\lVert \mathbf a \rVert} \end{pmatrix}\begin{pmatrix} \lVert \mathbf a \rVert \end{pmatrix}$
Supposons $\mathsf A \in \mathbf C^{n \times p}$ de la forme $\mathsf A = \left( \widetilde {\mathsf A} \; \mathbf a \right)$ et supposons que $\widetilde {\mathsf A} = \widetilde {\mathsf Q} \widetilde {\mathsf R}$, avec $\widetilde {\mathsf Q}$ orthogonale et $\widetilde {\mathsf R}$ triangulaire supérieure.

Alors \[\mathsf A = \begin{pmatrix} \widetilde {\mathsf Q} & \mathbf a \end{pmatrix} \begin{pmatrix} \widetilde {\mathsf R} & \mathbf 0 \\ \mathbf 0^T & 1 \end{pmatrix} = \begin{pmatrix} \widetilde {\mathsf Q} & \left(\mathbf a - \widetilde {\mathsf Q}\widetilde {\mathsf Q}^* \mathbf a\right) \end{pmatrix} \begin{pmatrix} \widetilde {\mathsf R} & \widetilde {\mathsf Q}^* \mathbf a \\ \mathbf 0^T & 1 \\ \end{pmatrix} \quad\textrm{motivé par}\quad \widetilde {\mathsf Q}^*\left(\mathbf a - \widetilde {\mathsf Q}\widetilde {\mathsf Q}^* \mathbf a\right)=\mathbf 0. \]

Posant $\mathbf q = \left(\mathbf a - \widetilde {\mathsf Q}\widetilde {\mathsf Q}^* \mathbf a\right) / \left\lVert \mathbf a - \widetilde {\mathsf Q}\widetilde {\mathsf Q}^* \mathbf a \right\rVert$, on obtient \[ \mathsf A = \underbrace{ \begin{pmatrix} \widetilde {\mathsf Q} & \mathbf q \end{pmatrix}}_{\text{orthogonale}} \underbrace{ \begin{pmatrix} \widetilde {\mathsf R} & \widetilde {\mathsf Q}^* \mathbf a \\ \mathbf 0 & {\left\lVert \mathbf a - \widetilde {\mathsf Q}\widetilde {\mathsf Q}^* \mathbf a \right\rVert} \end{pmatrix}}_{\text{triangulaire supérieure}} \quad\textrm{avec}\quad \widetilde {\mathsf Q}^* \mathbf q=\mathbf 0 \quad\textrm{et}\quad \mathsf Q {\mathsf Q}^*=\widetilde {\mathsf Q}\widetilde {\mathsf Q}^* + {\mathbf q} {\mathbf q}^*. \]

Itération simultanée

Supposons que $\mathsf X_0 \in \mathbf C^{n \times p}$ avec $p < n$ et considérons l’itération

\[ \fcolorbox{blue}{lightgreen}{$ \mathsf X_{k+1} = {\color{green}\text{normalisation}} (\mathsf A \mathsf X_{k}) $} \]

Pour la normalisation, appliquer Gram-Schmidt aux colonnes ;

En pratique, cela peut se faire par décomposition QR : \[ \mathsf X_{k+1} = \mathsf Q_{k+1}, \qquad \mathsf Q_{k+1} \mathsf R_{k+1} = \mathsf A \mathsf X_k \]
- $\mathsf Q_{k+1} \in \mathbf C^{n\times p}$ est une matrice unitaire : $\mathsf Q_{k+1}^* \mathsf Q_{k+1} = \mathsf I_{p \times p}$ ;
- $\mathsf R_{k+1} \in \mathbf C^{p\times p}$ est triangulaire supérieure.

Remarques

La première colonne effectue une simple itération des puissances ;
Mathématiquement, il est équivalent d’appliquer la décomposition QR à la fin. En effet \[ \mathsf X_{k} \underbrace{\mathsf R_{k} \mathsf R_{k-1} \dotsc \mathsf R_1}_{\text{tri. sup.}} = \mathsf A \mathsf X_{k-1} \mathsf R_{k-1} \mathsf R_{k-2} \dotsc \mathsf R_{1} = \dots = \mathsf A^k \mathsf X_0. \] ${\color{green} \rightarrow}$ en particulier, on a $\mathop{\mathrm{col}}(\mathsf X_k) = \mathop{\mathrm{col}}(\mathsf A^k \mathsf X_0)$.

Théorème : Convergence de l’itération simultanée

Supposons que

$\mathsf A \in \mathbf C^{n \times n}$ est hermitienne
$\mathsf X_0 \in \mathbf C^{n \times p}$ a des colonnes linéairement indépendantes
Le sous-espace engendré par les colonnes de $\mathsf X_0$ vérifie \[ \mathop{\mathrm{col}}(\mathsf X_0) \cap \mathop{\mathrm{Span}}\{ \mathbf{\boldsymbol{v}}_{p+1}, \dotsc, \mathbf{\boldsymbol{v}}_n \} = \varnothing. \]
Les $p+1$ premières valeurs propres sont distinctes : \[ |\lambda_1| > |\lambda_2| > \dotsb > |\lambda_p| > |\lambda_{p+1}| \geqslant|\lambda_{p+2}| \geqslant\dotsc \geqslant|\lambda_n|, \] Alors les colonnes de $\mathsf X_k$ convergent essentiellement vers $\mathbf{\boldsymbol{v}}_1, \dots, \mathbf{\boldsymbol{v}}_p$.

Application : Marche aléatoire sur un graphe

Considérons une marche aléatoire sur le graphe :

Hypothèse : les sauts vers les nœuds adjacents sont équiprobables.

Formulation du problème

La matrice d’adjacence associée $\mathsf A \in \{0, 1\}^{n \times n}$ est donnée par $a_{ij} = \begin{cases} 1 & \text{s'il existe une arête entre $i$ et $j$,} \\ 0 & \text{sinon.} \end{cases}$

Soit $\left(X_k\right)_{k \in \mathbf N}$ une chaîne de Markov à valeurs dans les nœuds telle que la probabilité de transition $i \rightarrow j$ est notée

\[ \fcolorbox{blue}{lightgreen}{$p_{ij} = \mathbf P{(X_{k+1}=j | X_k=i)}$} \quad\textrm{≠0 s'il existe une arête entre $i$ et $j$} \]

En supposant l’équiprobabilité des sauts vers les nœuds adjacents, on a $p_{ij} = \frac{a_{ij}}{n_e(i)}$ où $n_e(i)$ est le nombre d’arêtes connectées au nœud $i$.

${\color{green} \rightarrow}$ noter que $\mathsf A$ est symétrique mais pas $\mathsf P := (p_{ij})$.

La probabilité marginale de $X_{k+1}$ s’écrit par la loi des probabilités totales

\[ \fcolorbox{blue}{lightgreen}{$ \mathbf P{(X_{k+1}=i)}=\sum_{j=1}^n \mathbf P{(X_{k+1}=i | X_k=j)} \mathbf P{(X_{k}=j)} =\sum_{j=1}^n p_{ji} \mathbf P{(X_{k}=j)} $} \]

${\color{green} \rightarrow}$ la question est celle de l’état stationnaire de la distribution, c.-à-d. $\lim_{k\to\infty} \mathbf P{(X_{k}=i)}$

La distribution de probabilité stationnaire $\mathbf{\boldsymbol{v}}_1 \in \mathbf R_{\geqslant 0}^n$ vérifie \[ \mathsf P^T\mathbf{\boldsymbol{v}}_1 = \mathbf{\boldsymbol{v}}_1, \qquad \mathbf{\boldsymbol{1}}^T\mathbf{\boldsymbol{v}}_1 = 1. \] (il est possible de montrer que $1$ est la valeur propre dominante)

${\color{green} \rightarrow}$ problème aux valeurs propres pour la matrice $\mathsf P^T$.

Construction de la matrice d’adjacence

function adj(m)
  x = vcat([collect(0:i) for i in 0:m-1]...)
  y = vcat([fill(m-i, i) for i in 1:m]...)
  row_lengths = 1:m
  row_indices = [0; cumsum(row_lengths)]
  A = zeros(Int, length(x), length(x))
  for i in 1:m-1
      for j in 1:row_lengths[i]
          r1 = row_indices[i]
          r2 = row_indices[i+1]
          A[r2 + j, r2 + j + 1] = 1  # Horizontal
          A[r1 + j, r2 + j]     = 1  # Vertical
          A[r1 + j, r2 + j + 1] = 1  # Diagonal
      end
  end
  x, y, A + A'
end

Ici ... est l’opérateur splat (?... pour plus d’informations)

sum3(a, b, c) = a + b + c
numbers = [1, 2, 3]
sum3(numbers...)

Fonctionne aussi avec les arguments nommés (utiliser ; et NamedTuple)

plotargs = (color=:blue, label="example")
plot(sin; plotargs...)

Exercice : écrire la fonction transition(m) qui construit $\mathsf P$.

Remarquer que $\mathsf A$ est une matrice creuse :

# Matrice d'adjacence
x, y, A = adj(10);
Plots.spy(A, xlabel=L"i", ylabel=L"j")
Plots.plot!(size=(600, 600))

C’est le cas dans de nombreuses applications :

Discrétisation des équations aux dérivées partielles ;
Applications aux graphes.

Calcul de l’état stationnaire

La distribution de probabilité stationnaire vérifie \[ \mathsf P^T\mathbf{\boldsymbol{v}}_1 = \mathbf{\boldsymbol{v}}_1 \]

Illustration avec des unités en ‰, arrondies aux entiers :

using LinearAlgebra
v₁ = eigvecs(P', sortby=real)[:, end] |> real
v₁ = v₁ / sum(v₁)  # Normalisation
plotgraph(v₁)

Stratégie pour trouver $\mathbf{\boldsymbol{v}}_1$ : faire évoluer le vecteur de probabilité jusqu’à l’équilibre : \[ \mathbf{\boldsymbol{x}}^{(k+1)} = \mathsf P^T\mathbf{\boldsymbol{x}}^{(k)} \qquad \Leftrightarrow \qquad x^{(k+1)}_i = \sum_{j=1}^{n_n} p_{ji} x^{(k)}_j \]

$x^{(k)}_i$ est la probabilité de se trouver au nœud $i$ à l’itération $k$ ;
C’est la méthode des puissances ;
Ci-dessous, on prend $\mathbf{\boldsymbol{x}}^{(0)} = \mathbf{\boldsymbol{e}}_1$.

Considérations d’efficacité

Considérons l’implémentation suivante de la méthode des puissances

function power_iteration(M, x, niter)
    for i in 1:niter
        x = M*x
        # La normalisation n'est pas nécessaire ici (pourquoi ?)
    end
    return x
end

# Matrice de transition dense transposée
Pᵗ = transition(20) |> transpose |> Matrix;

# Itéré initial
nn = size(Pᵗ, 1) ; x = zeros(nn) ; x[1] = 1;

$~$

Approche naïve : utiliser une structure de données dense pour $\mathsf P^T$

using CountFlops
flops = @count_ops power_iteration(Pᵗ, x, 10) # only Julia ≤ 1.11
memory = sizeof(Pᵗ)
println(flops, "\n", "Mémoire : $memory octets")

add64: 441000
mul64: 441000

Mémoire : 352800 octets

${\color{green} \rightarrow}$ du temps est perdu sur les 0. Essayer zeros(10_000, 10_000)^2.

Meilleure solution : utiliser une structure de données creuse

import Base: *, transpose
struct coo_matrix
    rows::Array{Int}       # Indices de ligne des ≠0
    cols::Array{Int}       # Indices de colonne des ≠0
    vals::Array{Float64}   # Valeurs des entrées ≠0
    m::Int                 # Nombre de lignes
    n::Int                 # Nombre de colonnes
end

function transpose(A::coo_matrix)
    coo_matrix(A.cols, A.rows, A.vals, A.n, A.m)
end

${\color{green} \rightarrow}$ Exercice : définir les fonctions * et to_coo

Considérations d’efficacité

Considérons l’implémentation suivante de la méthode des puissances

function power_iteration(M, x, niter)
    for i in 1:niter
        x = M*x
        # La normalisation n'est pas nécessaire ici (pourquoi ?)
    end
    return x
end

# Matrice de transition dense transposée
Pᵗ = transition(20) |> transpose |> Matrix;

# Itéré initial
nn = size(Pᵗ, 1) ; x = zeros(nn) ; x[1] = 1;

$~$

Approche naïve : utiliser une structure de données dense pour $\mathsf P^T$

using CountFlops
flops = @count_ops power_iteration(Pᵗ, x, 10) # only Julia ≤ 1.11
memory = sizeof(Pᵗ)
println(flops, "\n", "Mémoire : $memory octets")

add64: 441000
mul64: 441000

Mémoire : 352800 octets

${\color{green} \rightarrow}$ du temps est perdu sur les 0. Essayer `zeros(10_000, 10_000)^2$.

Meilleure solution : utiliser une structure de données creuse

using CountFlops
sparse_Pᵗ = to_coo(Pᵗ)
flops = @count_ops power_iteration(sparse_Pᵗ, x, 10) # only Julia ≤ 1.11
memory = sizeof(sparse_Pᵗ)
println(flops, "\n", "Mémoire : $memory octets")

add64: 11400
mul64: 11400

Mémoire : 40 octets

Formats de matrices creuses

Nous prenons la matrice suivante pour illustration : \[ M = \begin{pmatrix} 5 & . & . \\ . & 6 & 7 \\ 8 & . & 9 \end{pmatrix} \]

Liste de coordonnées (COO), présentée sur la diapositive précédente ;
- rows vaut [1, 2, 2, 3, 3]
- cols vaut [1, 2, 3, 1, 3]
- vals vaut [5, 6, 7, 8, 9]

Lignes compressées (CSR, pour information seulement) :
- Similaire à COO, mais souvent plus économe en mémoire : rows est compressé ;
- Le vecteur rows contient les indices, dans cols et vals, où commence chaque ligne…
- … sauf pour rows[m+1], qui vaut 1 + le nombre d’éléments non nuls ;
- rows vaut [1, 2, 4, 6]
- cols vaut [1, 2, 3, 1, 3]
- vals vaut [5, 6, 7, 8, 9]

Colonnes compressées (CSC, pour information seulement) :
- Identique à CSR avec lignes $\leftrightarrow$ colonnes ;
- Implémenté en Julia par SparseArrays.

using SparseArrays
Z = spzeros(3, 3)

3×3 SparseMatrixCSC{Float64, Int64} with 0 stored entries:
  ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅

$~$

M = [5 0 0; 0 6 7; 8 0 9]
M = sparse(M)

3×3 SparseMatrixCSC{Int64, Int64} with 5 stored entries:
 5  ⋅  ⋅
 ⋅  6  7
 8  ⋅  9

$~$

# Constructeur utilisant le format COO
M = sparse([1, 2, 2, 3, 3],
           [1, 2, 3, 1, 3],
           [5, 6, 7, 8, 9])

3×3 SparseMatrixCSC{Int64, Int64} with 5 stored entries:
 5  ⋅  ⋅
 ⋅  6  7
 8  ⋅  9

$~$

println(M.rowval)
println(M.colptr)
println(M.nzval)

[1, 3, 2, 2, 3]
[1, 3, 4, 6]
[5, 8, 6, 7, 9]

Voir aussi sprand et spdiagm.