Rappels : normes vectorielles et matricielles (1/3)

Plots.plot(aspect_ratio=:equal, xlims=(-1.1,1.1), ylims=(-1.1, 1.1))
Plots.plot!(framestyle=:origin, grid=true, legend=:outertopright)
Plots.plot!(title="Cercle unité dans différentes normes")
Plots.plot!([0, 1, 0, -1, 0], [-1, 0, 1, 0, -1], label=L"$p = 1$ (norme taxicab)")
Plots.plot!(t -> cos(t), t -> sin(t), 0, 2π, label=L"$p = 2$ (norme euclidienne)")
Plots.plot!([1, 1, -1, -1, 1], [1, -1, -1, 1, 1], label=L"p = \infty")

Rappels : normes vectorielles et matricielles (2/3)

La norme d’opérateur sur les matrices induite par la \(p\)-norme vectorielle est donnée par \[ {\lVert {\mathsf A} \rVert}_{p} := \sup_{{\lVert {\mathbf{\boldsymbol{x}}} \rVert}_{p} \leqslant 1} {\lVert {\mathsf A \mathbf{\boldsymbol{x}}} \rVert}_{p} \] Il découle de la définition que \({\lVert {\mathsf A \mathsf B} \rVert}_p \leqslant{\lVert {\mathsf A} \rVert}_p {\lVert {\mathsf B} \rVert}_p\). En Julia, calculer \({\lVert {\mathsf A} \rVert}_p\) avec opnorm(A, p).

Rappels : normes vectorielles et matricielles (3/3)

À partir des normes matricielles, on définit le nombre de conditionnement d’une matrice par \[ \kappa_p(\mathsf A) = {\lVert {\mathsf A} \rVert}_p {\lVert {\mathsf A^{-1}} \rVert}_p \]

Conditionnement des systèmes linéaires

Illustration numérique

Considérons le système linéaire \[ \mathsf A(α) \mathbf x := \begin{pmatrix} 1 & 1 \\ 1 & 1 - α \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} π \\ π - πα \end{pmatrix} =: \mathbf b(α). \]

On trace sur le même graphe les fonctions qui à \(\alpha\) associent …

• … la quantité \(\kappa_2\bigl(\mathsf A(α)\bigr) \times ε\), où \(ε\) est l’epsilon machine pour le format Float64.

• … l’erreur relative en norme euclidienne, obtenue en résolvant le système linéaire avec l’opérateur \.

Af(α) = [1. 1.; 1. (1-α)]
bf(α) = [π; π - π*α]

# Solution exacte
x_exact = [0.; π]

# Plage de α
αs = 10. .^ (-15:.1:-2)
rel_err(α) = norm(Af(α)\bf(α) - x_exact)/norm(x_exact)
Plots.plot(αs, cond.(Af.(αs)) * eps(), label=L"κ_2(A) \times ε")
Plots.plot!(αs, rel_err.(αs), label=L"erreur $ℓ^2$")
Plots.plot!(xscale=:log10, yscale=:log10, xlabel=L"α", bottom_margin=5Plots.mm)

Décomposition LU par élimination de Gauss

\[ \begin{pmatrix} 2 & 1 & -1 \\ -3 & -1 & 2 \\ -2 & 1 & 2 \end{pmatrix} =: \mathsf A \]

Décomposition LU : coût de calcul

# A est une matrice inversible de taille n x n
L = [r == c ? 1.0 : 0.0 for r in 1:n, c in 1:n]
U = copy(A)
for c in 1:n-1, r in c+1:n
   U[c, c] == 0 && error("L'entrée pivot est nulle !")
   L[r, c] = U[r, c] / U[c, c]
   U[r, c:end] -= U[c, c:end] * L[r, c]
end
# L est triangulaire inférieure unitaire et U est triangulaire supérieure

\(~\)

Comparaison avec d’autres approches directes

using LinearAlgebra
A = [1 1 1; 1 2 3; 1 4 9]
A = factorize(A)
@show A  # Remarque : A peut encore être utilisée comme matrice par la suite !

A = LU{Float64, Matrix{Float64}, Vector{Int64}}([1.0 1.0 1.0; 1.0 3.0 8.0; 1.0 0.3333333333333333 -0.6666666666666665], [1, 3, 3], 0)

LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
 1.0  0.0       0.0
 1.0  1.0       0.0
 1.0  0.333333  1.0
U factor:
3×3 Matrix{Float64}:
 1.0  1.0   1.0
 0.0  3.0   8.0
 0.0  0.0  -0.666667

Méthode directe pour \(\mathsf A\) symétrique définie positive

• Il est plus efficace de calculer la factorisation de Cholesky \(\mathsf A = \mathsf C \mathsf C^*\) avec \(\mathsf C\) triangulaire inférieure à éléments diagonaux positifs.

function cholesky(A)
    n = size(A, 1)
    C = copy(LowerTriangular(A))
    for j ∈ 1:n
        Cⱼₖ = C[j, 1:j-1]
        C[j,j] -= Cⱼₖ'Cⱼₖ
        C[j,j] = √C[j,j]
        for i ∈ j+1:n
            C[i,j] -= Cⱼₖ'C[i, 1:j-1]
            C[i,j] = C[i,j]/C[j,j]
        end
    end
    return C
end

using Polynomials, Statistics
generate_defpos_matrix(n) = begin A = rand(n,n); A'A + I end
generate_defpos_matrix(1000) |> cholesky # Pour forcer la compilation

nb_samples = 10
Plots.plot(xaxis=:log10, yaxis=:log10, xlabel="n", ylabel="Temps CPU", legend=:topleft, size=(900,400),
           bottom_margin=5Plots.mm, left_margin=5Plots.mm, top_margin=5Plots.mm)
tf, tb = Float64[], Float64[]
tn = 2 .^(6:9)
for n in tn
    A = generate_defpos_matrix(n)
    push!(tf, mean(@elapsed cholesky(A) for _ in 1:nb_samples))
end

Pf = fit(log10.(tn), log10.(tf), 1) ; af = round(coeffs(Pf)[2]; digits=2)
Plots.plot!(tn, tf, marker=:o, label=L"Matrice dense $n^{%$af}$")
Plots.xticks!(tn, [L"2^{%$p}" for p in 6:9])

Cholesky pour les matrices bandes

• \(\mathsf A\) est une matrice bande, c.-à-d. qu’il existe \(b\in{\mathbb{{N}}}\) tel que \(a_{ij}=0\) si \(|i-j|>b\).

function cholesky_banded(A, b)
    n = size(A, 1)
    C = copy(LowerTriangular(A))
    for j ∈ 1:n
        m, M = max(1,j-b), min(n,j+b)
        Cⱼₖ = C[j, m:j-1]
        C[j,j] -= Cⱼₖ'Cⱼₖ
        C[j,j] = √C[j,j]
        for i ∈ j+1:M
            C[i,j] -= Cⱼₖ'C[i, m:j-1]
            C[i,j] = C[i,j]/C[j,j]
        end
    end
    return C
end

function generate_banded_matrix(n, b)
    C = [j≤i≤j+b ? rand() : 0.0 for i in 1:n, j in 1:n]
    return C*C'+I
end

b = 4; A = generate_banded_matrix(10, b)
cholesky_banded(A, b) # Pour forcer la compilation

Plots.plot(xaxis=:log10, yaxis=:log10, xlabel="n", ylabel="Temps CPU", legend=:topleft, size=(900,400),
           bottom_margin=5Plots.mm, left_margin=5Plots.mm, top_margin=5Plots.mm)
tf, tb = Float64[], Float64[]
tn = 2 .^(6:9)
for n in tn
    A = generate_banded_matrix(n, b)
    push!(tf, mean(@elapsed cholesky(A) for _ in 1:nb_samples))
    push!(tb, mean(@elapsed cholesky_banded(A,b) for _ in 1:nb_samples))
end

Pf = fit(log10.(tn), log10.(tf), 1) ; af = round(coeffs(Pf)[2]; digits=2)
Plots.plot!(tn, tf, marker=:o, label=L"Matrice dense $n^{%$af}$")

ntn = 2 .^(10:12)
for n in ntn
    A = generate_banded_matrix(n, b)
    push!(tb, mean(@elapsed cholesky_banded(A,b) for _ in 1:nb_samples))
end
append!(tn,ntn)
Pb = fit(log10.(tn),log10.(tb),1) ; ab = round(coeffs(Pb)[2]; digits=2)
Plots.plot!(tn, tb, marker=:diamond, label="Matrice bande "*latexstring("n^{$(ab)}"))
Plots.xticks!(ntn, [L"2^{%$p}" for p in 6:12])

Une première méthode itérative

Rayon spectral

Remarques :

• \(\rho(A)\) n’est pas une norme ;

• \(\rho(A) \leqslant{\lVert {\mathsf A} \rVert}\) pour toute norme matricielle induite.

• 🔎 Formule de Gelfand : Pour tout \(\mathsf A \in \mathbf R^{n \times n}\) et toute norme matricielle \({\lVert {\cdot} \rVert}\), \[ \lim_{k \to \infty} {\lVert {\mathsf A^k} \rVert}^{1/k} = \rho(\mathsf A). \]

Analyse de la méthode de Richardson (1/2)

\[\begin{align*} \mathbf{\boldsymbol{x}}^{(k + 1)} = \mathbf{\boldsymbol{x}}^{(k)} + \omega (\mathbf{\boldsymbol{b}} - \mathsf A \mathbf{\boldsymbol{x}}^{(k)}) \end{align*}\]

Lien avec l’optimisation

\(~\)

Illustration de la méthode de Richardson

Illustration de la méthode de Richardson : \(\omega = .2\)

Illustration de la méthode de Richardson : \(\omega = .6\)

Illustration de la méthode de Richardson : \(\omega\) optimal \(= .5\)

🔎 Méthodes itératives de base plus générales

Considérons la décomposition \[ \mathsf A = {\color{green}\mathsf M} - {\color{blue} \mathsf N} \] où \(\mathsf M\) est inversible et facile à inverser. Alors \[ \mathsf A \mathbf{\boldsymbol{x}}_* = \mathbf{\boldsymbol{b}} \quad \Leftrightarrow \quad {\color{green}\mathsf M} \mathbf{\boldsymbol{x}}_* = {\color{blue} \mathsf N} \mathbf{\boldsymbol{x}}_* + \mathbf{\boldsymbol{b}} \] ce qui suggère la méthode itérative \[ {\color{green}\mathsf M} \mathbf{\boldsymbol{x}}^{(k+1)} = {\color{blue} \mathsf N} \mathbf{\boldsymbol{x}}^{(k)} + \mathbf{\boldsymbol{b}} \]

🔎 Convergence de la méthode itérative de base générale : solution

🔎 Garanties de convergence concrètes

Les situations où \(\rho(\mathsf M^{-1} \mathsf N) < 1\) sont identifiées au cas par cas :

• Itération de Richardson, condition suffisante : \(\mathsf A\) symétrique définie positive ;

• Itération de Jacobi, condition suffisante : \(\mathsf A\) strictement à diagonale dominante par lignes ou par colonnes :

\[ \lvert a_{ii} \rvert > \sum_{j \neq i} \lvert a_{ij} \rvert \quad \forall i \qquad \text{ ou } \qquad \lvert a_{jj} \rvert > \sum_{i \neq j} \lvert a_{ij} \rvert \quad \forall j. \]

• Itération de Gauss-Seidel, condition suffisante : \(\mathsf A\) strictement à diagonale dominante par lignes ou par colonnes ;

• Méthode de relaxation, condition suffisante : \(\mathsf A\) symétrique définie positive et \(\omega \in (0, 2)\) ;

• Méthode de relaxation, condition nécessaire : \(\omega \in (0, 2)\) ;

• …

Améliorer la méthode de Richardson

Rappel : jusqu’à la fin de ce cours, la matrice \(\mathsf A\) est supposée symétrique définie positive

On rappelle l’itération de Richardson : \[\begin{align*} \mathbf{\boldsymbol{x}}^{(k + 1)} = \mathbf{\boldsymbol{x}}^{(k)} - \omega \nabla f(\mathbf{\boldsymbol{x}}^{(k)}), \qquad f(\mathbf{\boldsymbol{x}}) := \frac{1}{2} \mathbf{\boldsymbol{x}}^T\mathsf A \mathbf{\boldsymbol{x}} - \mathbf{\boldsymbol{b}}^T\mathbf{\boldsymbol{x}}. \end{align*}\]

On peut améliorer la méthode de Richardson en prenant \[ \mathbf{\boldsymbol{x}}^{(k + 1)} = \mathbf{\boldsymbol{x}}^{(k)} - \omega_{\color{red}k} \nabla f(\mathbf{\boldsymbol{x}}^{(k)}), \qquad \omega_{\color{red}k} := \frac{\mathbf{\boldsymbol{d}}_k^T\mathbf{\boldsymbol{d}}_k}{\mathbf{\boldsymbol{d}}_k^T\mathsf A \mathbf{\boldsymbol{d}}_k}, \qquad \mathbf{\boldsymbol{d}}_k := \mathsf A \mathbf{\boldsymbol{x}}^{(k)} - \mathbf{\boldsymbol{b}}. \]

C’est la méthode de descente de gradient avec pas optimal.

Illustration de la descente de gradient

Amélioration supplémentaire : la méthode des directions conjuguées

\[ f(\mathbf{\boldsymbol{x}}) := \frac{1}{2} \mathbf{\boldsymbol{x}}^T\mathsf A \mathbf{\boldsymbol{x}} - \mathbf{\boldsymbol{b}}^T\mathbf{\boldsymbol{x}} = \frac{1}{2} {\lVert {\mathbf{\boldsymbol{x}} - \mathbf{\boldsymbol{x}}_*} \rVert}_{\mathsf A}^2 + \text{constante}, \qquad {\lVert {\mathbf{\boldsymbol{y}}} \rVert}_{\mathsf A} := \sqrt{{\langle {\mathbf{\boldsymbol{y}}, \mathbf{\boldsymbol{y}}} \rangle}_{\mathsf A}}, \qquad {\langle {\mathbf{\boldsymbol{x}}, \mathbf{\boldsymbol{y}}} \rangle}_{\mathsf A} := \mathbf{\boldsymbol{x}}^T\mathsf A \mathbf{\boldsymbol{y}}. \] Étant donné \(n\) directions conjuguées \(\{\mathbf{\boldsymbol{d}}_0, \dotsc, \mathbf{\boldsymbol{d}}_{n-1}\}\) telles que \({\langle {\mathbf{\boldsymbol{d}}_i, \mathbf{\boldsymbol{d}}_j} \rangle}_{\mathsf A} = \delta_{ij}\), on a \[ \mathbf{\boldsymbol{x}}_* - \mathbf{\boldsymbol{x}}^{(0)} = \sum_{i=0}^{n-1} \mathbf{\boldsymbol{d}}_i \mathbf{\boldsymbol{d}}_i^T\mathsf A (\mathbf{\boldsymbol{x}}_* - \mathbf{\boldsymbol{x}}^{(0)}) = \sum_{i=0}^{n-1} \mathbf{\boldsymbol{d}}_i \mathbf{\boldsymbol{d}}_i^T({\color{green}\mathbf{\boldsymbol{b}}} - \mathsf A \mathbf{\boldsymbol{x}}^{(0)}) \]

# On suppose que `ds` est une liste de n directions conjuguées
x = zeros(n)
for k in 1:n
    d = ds[k]
    r = A*x - b
    ω = d'r / (d'A*d)  # Dénominateur = 1 si `d` est normalisé
    x -= ω * d
end

Directions conjuguées : exemple

• Considérons à nouveau le système linéaire \[ \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} 7 \\ 8 \end{pmatrix} \]

• Directions conjuguées (non normalisées) calculées par Gram-Schmidt : \[ \mathbf{\boldsymbol{d}}_0 = \begin{pmatrix} 1 \\ 0 \end{pmatrix}, \qquad \mathbf{\boldsymbol{d}}_1 = \begin{pmatrix} -1 \\ 2 \end{pmatrix}. \]

using LinearAlgebra
n = 2
A = [2.0 1.0; 1.0 2.0]
b = [7.0; 8.0]
x = zeros(n)
ds = [[1.0; 0.0], [-1.0; 2.0]]
for i in 1:n
    d = ds[i]
    r = A*x - b
    ω = d'r / (d'A*d)  # Dénominateur utile si `d` n'est pas normalisé
    x -= ω*d
    println("Résidu : $(norm(A*x - b))")
end

Résidu : 4.5
Résidu : 0.0

Directions conjuguées : illustration

Illustration de la méthode du gradient conjugué