Décomposition de Cholesky

La décomposition de Cholesky est une décomposition matricielle d’une matrice symétrique définie positive $A$ sous la forme d’un produit $T \cdot T^t$ où $T$ est une matrice triangulaire inférieure, dont les coefficients sont obtenus à l’aide des équations suivantes :

$$ t_{i,i}^2 = a_{i,i} - \sum_{k=1}^{i-1} t_{i,k}^2 $$

$$ t_{j,i} = \frac{a_{i,j} - \sum_{k=1}^{i-1} t_{i,k}t_{j,k}}{t_{i,i}} ; \forall j\geq i $$

Le résultat de la factorisation est la matrice $T$. Ceci permet de voir que toute matrice symétrique définie positive peut se mettre sous la forme d’un produit $T \cdot T^t$.

Factorisation complète

  • Questions :
    1. Écrire l’algorithme de factorisation dense de Cholesky et le tester. Quelle est sa complexité ?
    2. En utilisant cet algorithme, combien coûte une résolution de système linéaire dense $A.x = b$ ?

Factorisation incomplète

La factorisation dense de Cholesky, bien qu’il s’agisse d’un algorithme privilégié pour la résolution de systèmes linéaires symétriques (du fait de sa stabilité numérique), est encore trop coûteuse lorsque la matrice est creuse (i.e ne contient que peu de valeurs non nulles). Il existe un algorithme de factorisation de Cholesky dans le cas creux, mais il est complexe à implémenter (en particulier pour gérer les termes de remplissage). Alternativement, la factorisation de Cholesky incomplète consiste à ne calculer qu’un sous-ensemble des éléments de la matrice $T$.

Dans ce projet, la factorisation incomplète correspond à la factorisation classique, dans laquelle on ne calcule pas les éléments $T_{i,j}$ lorsque $A_{i,j}$ est nul. Ces éléments restent alors à zéro et n’influent pas sur la suite du calcul de la factorisation.

  • Questions :
    1. Écrire un algorithme permettant de générer des matrices symétriques définies positives creuses avec un nombre de termes extra-diagonaux non nuls réglable.
    2. Écrire l’algorithme de factorisation de Cholesky incomplète et le tester. Quels tests réaliser ?.
    3. Quelle est sa complexité ? Penser à choisir de manière appropriée en fonction de quels paramètres établir cette complexité et expliquer en quoi elle est meilleure que celle de l’algorithme précédent.
    4. Un préconditionneur $M$ est une matrice facile à inverser telle que $\mathrm{cond}(M^{-1}A)$ soit inférieur à $\mathrm{cond}(A)$, où $\mathrm{cond}(A)$ est le conditionnement de la matrice $A$. En général, on choisit $M$ comme un inverse approché de $A$ ($M^{-1} \approx A^{-1}$) relativement facile à calculer. Évaluer si la matrice ${T^{-1}}^t \cdot T^{-1}$ obtenue dans les deux cas précédents (factorisations complète et incomplète) est un préconditionneur de bonne ou de mauvaise qualité. Vvous pouvez à utiliser np.linalg.cond et, exceptionnement pour cette question, vous pouvez calculer l’inverse de $A$ avec np.inv(A).
  • Remarque 1 : Pour les évaluations de complexité, on pourra considérer que le produit matrices creuse / vecteur, ainsi que la résolution d’un système linéaire triangulaire supérieure creux sont de complexité linéaire en le nombre d’éléments non nuls de la matrice.
  • Remarque 2 : La matrice $A$ étant factorisée (de manière complète ou incomplète) sous la forme $T \cdot T^t$, pour résoudre un système $A.x =b$ :
    • on a alors $T . T^t . x = b$, on fait une résolution triangulaire : $T . y = b$
    • puis une seconde résolution triangulaire : $T^t . x = y$

D’une manière générale, on a rarement besoin de calculer l’inverse d’une matrice $A^{-1}$ !