# CAMA : ma13 Système matriciel -- Exercices
# Cours du 04/05

## Programmation vectorielle
:::danger
Le but des exercices est
* d'avoir un programme qui donne la bonne réponse
* qui soit le plus rapide possible (et pour cela on utilise massivement Numpy)
:::
:::warning
En règle général si vous avez des `for` imbriqués c'est mauvais signe.
:::

## Méthode du pivot de Gauss partiel
L'ennoncé est dans le cours

:::spoiler Solution

def solve_gauss_partial(A, b):   # on prend le max dans la colonne i parmi les lignes en dessous (plus facile)
    for i in range(len(A)-1):
        pivot = np.argmax(np.abs(A[i:, i]))  # il n'y a que 3 lignes à ajouter pour échanger les lignes
        A[[i, pivot]] = A[[pivot, i]]
        b[[i, pivot]] = b[[pivot, i]]
        E = np.diag(np.array([1.,] * len(A), dtype=A.dtype))
        coefs = - A[i+1:,i] / A[i,i]
        E[i+1:,i] = coefs
        A[i:, i:] = E[i:,i:] @ A[i:,i:]
        b[i+1:] += coefs * b[i]   # multiplication terme à terme
    # A est maintenant triangulaire surpérieur
    res = np.zeros(len(b), dtype=b.dtype)
    res[-1] = b[-1] / A[-1,-1]
    for i in range(len(A)-1)[::-1]:
        res[i] = (b[i] - A[i,i+1:] @ res[i+1:]) / A[i,i]
    return res

e = 1E-6
A = np.array([[e, 1], [1, 2]], dtype='float32')
b = np.array([1., 3.], dtype='float32')
print(f"A\n {A} \nb\n {b}\n")
x = solve_gauss_partial(A, b)
print('solution : ',x)
print('vérification\n', A@x)

A
 [[0.000001 1.      ]
 [1.       2.      ]] 
b
 [1. 3.]

solution :  [1.0000019  0.99999905]
vérification
 [3.       0.999997]

:::

## Factorisation de Choleski
Il s'agit de décomposer A en $A = B\, B^T$ avec B une matrice triangulaire inférieure. Cela n'est possible
que si la matrice A est symétrique et définie positive (c'est d'ailleurs un facon de vérifier qu'une
matrice est définie positive).

Écrire l'algorithme de Choleski qui prend A et retourne B (pour deviner l'algorithme, essayez de trouver les
coefficients de B à partir des coefficients de A sur une matrice A 4x4).

:::spoiler Solution
$$
A = B\, B^T =
\begin{bmatrix}
b_{11} & 0 & \dots & 0\\
b_{21} & b_{22} & \dots & 0\\
& \vdots&\\
b_{n1} & b_{n2} & \dots& b_{n,n}
\end{bmatrix}
\begin{bmatrix}
b_{11} & b_{21} & \dots & b_{n1}\\
0 & b_{22} & \dots & b_{n2}\\
& \vdots&\\
b_{n1} & b_{n2} & \dots& b_{n,n}
\end{bmatrix}=
\begin{bmatrix}
b_{11}^2 & b_{11}b_{21} & \dots & b_{11}b_{n1}\\
x & \sum_{i=1}^2b_{2i}^2 & \dots & \sum_{i=1}^2b_{2i}b_{ni}\\
& & \vdots&\\
x & x & \dots& \sum_{i=1}^2b_{n,i}^2
\end{bmatrix}
$$
avec $x$ la même valeur que de l'autre coté de la diagonale

On voit que $b_{11} = \sqrt{a_{11}}$ et maintenant qu'on a $b_{11}$ on peut trouver toute la première ligne de $B^T$ : $b_{j1}=a_{1j}/b_{11}$.

Une fois qu'on connait la première ligne de $B^T$ , on s'attaque à la deuxième en commencant par trouver $b_{22}$ puis ensuite tous les autres éléments de la ligne comme on a fait pour la première ligne.

On a donc dans le cas général :
* $b_{ii} = \sqrt{a_{ii} - \sum_{k=1}^{i-1}b_{ik}^2}$
* $b_{ji} = a_{ij} - \sum_{k=1}^{i-1}b_{ik}b_{jk}/b_{ii} = a_{ij} - \sum_{k=1}^{i-1}b_{ik}b_{kj}^T/b_{ii} \space\forall j\gt i$

def Choleski(A):
    B = np.zeros(A.shape)
    for i in range(len(A)):
        B[i,i] = np.sqrt(A[i,i] - np.sum(np.square(B[i, :i])))         # garanti ok car A est def positive
        B[i+1:, i] = (A[i, i+1:] - B[i, :i] @ B.T[:i, i+1:]) / B[i,i]  # les ∑ sous forme de prod. scalaire
    return B

:::

## Matrice symetrique
Rappel : pas de boucles `for` imbriquées mais des opérations vectorielles et matricielles (sur des sous-matrices).

Créer une matrice symétrique définie positive est vérifier que votre programme marche bien.

:::spoiler Solution

A = np.random.randint(10, size=(4,4))
A = A + A.T                             # symmétrique
A = A + np.diag(A.sum(axis=0))          # diagonale dominante
print('A:\n', A)
B = Choleski(A)
print('B\n', B)
print('vérification\n', B @ B.T)

A:
 [[55  8 18  5]
 [ 8 33  7 10]
 [18  7 54  9]
 [ 5 10  9 28]]
B
 [[7.4161984871 0.           0.           0.          ]
 [1.0787197799 5.6423721639 0.           0.          ]
 [2.4271195049 0.7765914857 6.8924593995 0.          ]
 [0.6741998625 1.6434093681 0.8831939788 4.9055711788]]
vérification
 [[55.  8. 18.  5.]
 [ 8. 33.  7. 10.]
 [18.  7. 54.  9.]
 [ 5. 10.  9. 28.]]

:::

## Améliorer Jacobi
Lorsqu'on écrit une itération de la méthode de Jacobi avec l'ensemble des coefficients, on constate que
si on calcule la nouvelle valeur de **x** élément par élement alors lorsqu'on veut mettre à jour `x[1]`,
on connait déjà `x[0]`. Idem lorsqu'on met à jour `x[2]` on connait déjà `x[0]` et `x[1]`, etc.

L'idée de la version amméliorée de Jacobi est d'utiliser la nouvelle valeur de `x[0]` et non pas l'ancienne
comme c'est le cas dans l'algorithme du cours. Ainsi en utilisant des valeurs mise à jour on peut espérer
converger plus vite.

Dans ce chaque itération demande un calcul ligne par ligne et donc une boucle `for` dans la boucle `for` sur
les itérations.

#### Test d'arrêt

On ajoutera un argument `error` à la fonction qui indique la précision désirée du résultat. Par
défaut sa valeur est de `1E-6` et pour offrir une bonne garantie on arrête l'algorithme lorsque
$||x_{t+1} - x_t|| < \texttt{error}\, / \, 10$.

:::spoiler Solution

def Jacobi(A, b, error=1E-6, verbose=False):
    L = np.tril(A) 
    U = -np.triu(A, k=1) 
    if verbose:
        print(f"L:\n {L}\nU\n {U}\n")
    previous_x = np.zeros(len(b))
    x = np.random.random(len(b))
    err = (error / 10) ** 2
    while np.sum(np.square(x - previous_x))  > err:
        previous_x = x.copy()
        if verbose:
            print(f"x = {x}")
        # on résoud  L x = U x + b  avec L matrice triangulaire inférieure
        y = U @ x + b
        x[0] = y[0] / L[0,0]
        for i in range(1,len(L)):
            x[i] = (y[i] - L[i,:i] @ x[:i]) / L[i,i]
    return x

A = np.random.randint(10, size=(4,4))
A = A + np.diag(A.sum(axis=0))
b = A.sum(axis=1)                     # ainsi la solution est [1,1,1,1]
print('A:\n', A, "\nb:\n", b, "\n")

Jacobi(A,b, verbose=True)

A:
 [[24  2  1  7]
 [ 5 19  4  6]
 [ 9  2 20  9]
 [ 2  9  9 32]] 
b:
 [34 34 40 52] 

L:
 [[24  0  0  0]
 [ 5 19  0  0]
 [ 9  2 20  0]
 [ 2  9  9 32]]
U
 [[ 0 -2 -1 -7]
 [ 0  0 -4 -6]
 [ 0  0  0 -9]
 [ 0  0  0  0]]

x = [0.8870874823 0.8448958895 0.2146829205 0.8281640711]
x = [1.0957657001 1.194392389  1.0147923641 0.9351814319]
x = [1.0020897014 1.0168049182 1.0265474982 0.9876765266]
x = [1.0010877908 0.9980164155 1.0052544156 0.9990120918]
x = [1.0002345046 0.9991440665 1.000424625  1.000106649 ]
x = [1.0000225291 0.9998709979 0.9999547701 1.0000475947]
x = [0.999998753  0.9999948204 0.9999796615 1.0000072549]
x = [0.9999991631 1.000002211  0.9999968908 1.0000003049]
x = [0.9999998564 1.0000005961 0.9999998678 0.9999998785]
x = [0.9999999913 1.0000000685 1.0000000518 0.9999999667]

array([1.0000000018, 0.9999999991, 1.0000000142, 0.9999999961])

:::