--- title: VISI601_CMI (5a) Représentation des matrices creuses type: slide slideOptions: transition: slide progress: true slideNumber: true --- # Représentation des matrices creuses > [name=Jacques-Olivier Lachaud][time=January 2022][color=#907bf7] > (Les images peuvent être soumises à des droits d'auteur. Elles sont utilisées ici exclusivement dans un but pédagogique) ###### tags: `visi601` Retour à [VISI601_CMI (Main) Algorithmique numérique](https://codimd.math.cnrs.fr/s/IWTaBkA9m) --- ## Taille des matrices et coûts de calcul 1. Le vecteur b représente des données issues * des images donc $n=10^6$ pour une image $1000 \times 1000$ * des maillages de centaines de milliers voire millions de sommets, $n \gg 10^5$ 2. La solution cherchée x est un vecteur de même taille 3. La matrice A correspond à des opérateurs liant les x * sa taille est $n \times n$ soit $\gg 10^{10}$ coefficients * coût mémoire $\gg$ 10 gigaoctets * coût résolution $\gg 10^{30}$ flops * même supercalculateur à $10^{18}$ flops/s prend $\gg$ 30 000 ans :::danger Mais les matrices $A$ représentent des *relations locales* $\Rightarrow$ elles ont bcp de coefficients à 0, on parle de **matrices creuses** ::: --- ## Exemple: système masse ressort (I) On fixe $x_0$ et $x_n$ à un mur, et on met $n$ ressorts entre. $\bullet\leftrightsquigarrow\bullet\leftrightsquigarrow\bullet\leftrightsquigarrow \cdots \bullet\leftrightsquigarrow\bullet$ $x_0 \quad x_1 \quad x_2 \quad \quad x_{n-1} \quad x_n$ Si $k$ raideur, alors la force en $x_i$ est $F_i=-k(x_{i+1}-x_i)-k(x_{i-1}-x_i)$ sauf aux bords: $F_0=-k(x_{1}-x_0)$ et $F_n=-k(x_{n-1}-x_n)$ On applique la loi de Newton: somme des forces = masse fois accélération $$ \mathbf{f}:= \underbrace{\begin{bmatrix} 2k & -k & 0 & \cdots & 0 \\ -k & 2k & -k & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -k \\ 0 & \cdots & 0 & -k & 2k \end{bmatrix}}_{kA} \underbrace{\begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_{n-2}\\ x_{n-1} \end{bmatrix}}_{\mathbf{x}}+\underbrace{\begin{bmatrix}-kx_0 \\0\\ \vdots \\ 0 \\ -kx_n \end{bmatrix}}_{k\mathbf{b}} $$ --- ## Exemple: système masse ressort (II) ==Position de repos== du système : $\mathbf{f}=0 \quad\Leftrightarrow\quad A\mathbf{x}=\mathbf{b}$ On voit que $A$ est **creux**. Pour $n=6$, $x_0=-5$ et $x_6=5$, on trouve $\mathbf{x}^T=(-10/3, -5/3, 0, 5/3, 10/3)$ :::danger $A$ creux, mais $A^{-1}$ n'est en général pas creux ! ::: $A^{-1}=\left(\begin{array}{rrrrr} \frac{5}{6} & \frac{2}{3} & \frac{1}{2} & \frac{1}{3} & \frac{1}{6} \\ \frac{2}{3} & \frac{4}{3} & 1 & \frac{2}{3} & \frac{1}{3} \\ \frac{1}{2} & 1 & \frac{3}{2} & 1 & \frac{1}{2} \\ \frac{1}{3} & \frac{2}{3} & 1 & \frac{4}{3} & \frac{2}{3} \\ \frac{1}{6} & \frac{1}{3} & \frac{1}{2} & \frac{2}{3} & \frac{5}{6} \end{array}\right)\quad$ donc ne jamais inverser $A$ ! --- ## Représentations creuses de matrices Soit $n_{nz}$ le nombre de coefficients non nuls de $A$ Il faut trouver des structures de données qui permettent: * d'économiser la mémoire utilisée, i.e. proportionnel à $n_{nz}$ * d'accéder efficacement à un élément * d'accéder efficacement à tous les éléments d'une ligne / d'une colonne * de modifier les éléments d'une ligne / d'une colonne * d'extraire une ligne / une colonne * faire des multiplications matrice / vecteur ou matrice / matrice efficaces ==Pas de structures de données efficace sur tout== --- ## Format COO (COOrdinates sparse matrix) COO: format triplet ou coordonnées. On stocke des triplets (valeur, colonne, ligne) * sert à la construction d'une matrice creuse. * (+) permet d'avoir des entrées dupliquées (qui sont additionnées) * (+) conversion rapide vers les autres formats (CSR/CSC) * (-) sinon, inefficace. --- [Python, Scipy.sparse.coo_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.coo_matrix.html) ```python= # Constructing a matrix with duplicate indices >>> row = np.array([0, 0, 1, 3, 1, 0, 0]) >>> col = np.array([0, 2, 1, 3, 1, 0, 0]) >>> data = np.array([1, 1, 1, 1, 1, 1, 1]) >>> coo = coo_matrix((data, (row, col)), shape=(4, 4)) >>> # Duplicate indices are maintained until implicitly or explicitly summed >>> np.max(coo.data) 1 >>> coo.toarray() array([[3, 0, 1, 0], [0, 2, 0, 0], [0, 0, 0, 0], [0, 0, 0, 1]]) ``` --- ## Format CSR (Compressed Sparse Row Matrix) * les lignes sont compressées, donc très rapides à parcourir * `data` : tableau des **valeurs** de taille $n_{nz}$, * `indices` : tableau des **colonnes** de chaque valeur, de taille $n_{nz}$ * `indptr` : tableau du numéro du 1er élement de la **ligne**, de taille $m+1$ * (+) opérations arithmétiques efficaces: CSR + CSR, CSR * CSR, etc. * (+) extraire une ligne efficace * (+) produit matrice vecteur efficace * (-) extraire une colonne inefficace (utiliser CSC) * (-) modifications assez inefficaces (voir LIL ou DOK) --- [Python Scipy scipy.sparse.csr_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csr_matrix.html) ```python= >>> data = np.array([1, 2, 3, 4, 5, 6]) >>> indices = np.array([0, 2, 2, 0, 1, 2]) >>> indptr = np.array([0, 2, 3, 6]) >>> csr_matrix((data, indices, indptr), shape=(3, 3)).toarray() array([[1, 0, 2], [0, 0, 3], [4, 5, 6]]) ``` --- ## Format CSC (Compressed Sparse Column Matrix) * les colonnes sont compressées, donc très rapides à parcourir * `data` : tableau des **valeurs** de taille $n_{nz}$, * `indices` : tableau des **lignes** de chaque valeur, de taille $n_{nz}$ * `indptr` : tableau du numéro du 1er élement de la **colonne**, de taille $m+1$ * (+) opérations arithmétiques efficaces: CSC + CSC, CSC * CSC, etc. * (+) extraire une colonne efficace * (+) produit matrice vecteur efficace (mais CSR plus rapide) * (-) extraire une ligne inefficace (utiliser CSR) * (-) modifications assez inefficaces (voir LIL ou DOK) --- [Python Scipy scipy.sparse.csc_matrix](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.csc_matrix.html) ```python= >>> data = np.array([1, 2, 3, 4, 5, 6]) >>> indices = np.array([0, 2, 2, 0, 1, 2]) >>> indptr = np.array([0, 2, 3, 6]) >>> csc_matrix((data, indices, indptr), shape=(3, 3)).toarray() array([[1, 0, 4], [0, 0, 5], [2, 3, 6]]) ``` --- ## Exercices 1. Ecrire une fonction `sommeLigne(A)` pour sommer tous les coefficients non nuls de la ligne $i$ d'une matrice $A$ codée en CSR 2. Ecrire une fonction `at(A,i,j)` qui retourne la valeur du coefficient $A_{ij}$ ($A$ CSR) 3. Ecrire la fonction `produit(A,v)`, qui réalise le produit $A\mathbf{v}$, où $A$ est une matrice CSR, $\mathbf{v}$ un vecteur de taille $n$. 4. Dans l'exemple du système masse-ressort, quel serait le coût de stockage de la matrice $A$ pleine ? Même question si elle est stockée sous forme CSR ?