Semaine 4 : Problèmes aux valeurs propres

using LinearAlgebra
using TestImages
using LaTeXStrings
using Plots
import Random

Valeur propre unique

[Exercice 1] Itération inverse et itération de Rayleigh

1. Étant donné une matrice $N\times N$ $A$ et un nombre complexe $\mu \notin \sigma(A)$, implémenter l'itération inverse pour approcher la valeur propre de $A$ la plus proche de $\mu$ et un vecteur propre associé, en utilisant comme critère d'arrêt que $$\|A\hat v - \hat\lambda\hat v\| \leq \varepsilon\|\hat v\|.$$

   La méthode prendra comme arguments la matrice $A$, le nombre $\mu$, le vecteur initial $x_0$, et comme arguments nommés la tolérance $\varepsilon$ et un nombre maximal d'itérations avant abandon. Renvoyer un triplet $(\lambda_n,v_n,n)$ avec la paire propre approximative (normalisée) et le nombre $n$ d'itérations effectuées.

   Indice (cliquer pour afficher)

   • Par calcul fonctionnel, $\lambda_j = \underset{\lambda\in \sigma(A)}{\mathrm{argmin}} |\lambda-\mu|$ si et seulement si $(\lambda_j-\mu)^{-1}$ est la valeur propre dominante de $(A-\mu I_N)^{-1}$
   • Il n'est pas nécessaire (ni efficace) de calculer l'inverse de $A-\mu I_N$ (lire la documentation de LinearAlgebra.factorize)
   • Il n'est pas nécessaire de créer explicitement la matrice identité ; pour construire $A - \mu I_N$, il suffit d'écrire A - μ*I (en supposant que LinearAlgebra a été importé avec using)

function inverse_iteration(A, μ, x₀; tol=1e-12, maxiter=100)
    x = copy(x₀)
    ### BEGIN SOLUTION
    M = factorize(A-μ*I)
    niter = 0
    while niter < maxiter
        niter += 1
        x .= M\x
        norm2_x = x'x
        norm_x = √norm2_x
        x /= norm_x
        λ = x'A*x
        (norm(A*x - λ*x) ≤ tol) && return (λ,x,niter)
    end
    ### END SOLUTION
    return nothing
end

inverse_iteration (generic function with 1 method)

2. Implémenter l'itération du quotient de Rayleigh, prenant comme arguments la matrice $A$, la valeur initiale $\mu_0$, le vecteur initial $x_0$, et comme arguments nommés la tolérance $\varepsilon$ et un nombre maximal d'itérations. Comme ci-dessus, renvoyer un triplet avec la paire propre et le nombre d'itérations effectuées, en utilisant le même critère d'arrêt.

   Indice (cliquer pour afficher)

   • Cette méthode consiste à fixer $\mu_k$ à l'estimation courante de la valeur propre cible, qui est le quotient de Rayleigh $\frac{v_k^* A v_k}{v_k^*v_k}$.
   • Là encore, il n'est pas nécessaire d'inverser $A-\mu_k I_N$. Faut-il utiliser LinearAlgebra.factorize ?

function rayleigh_iteration(A, μ₀, x₀; tol=1e-12, maxiter=100)
    x = copy(x₀)
    μ = μ₀
    niter = 0
    ### BEGIN SOLUTION
    while niter < maxiter
        niter += 1
        x .= (A - μ*I) \ x
        norm2_x = x'x
        norm_x = √norm2_x
        x /= norm_x
        μ = x'A*x
        (norm(A*x - μ*x) ≤ tol) && return (μ,x,niter)
    end
    ### END SOLUTION
    return nothing
end

rayleigh_iteration (generic function with 1 method)

Faire varier la valeur de $N$ dans le code ci-dessous (garder $N$ en dessous de 1000). Il n'y a aucune garantie que la méthode d'itération de Rayleigh converge vers la valeur propre la plus proche de $\mu_0$, donc il peut être nécessaire d'exécuter la cellule plusieurs fois pour que cela se produise. La méthode du quotient de Rayleigh est-elle plus rapide en nombre d'itérations ? En temps d'exécution ? Pourquoi ?

# Test code
Random.seed!(2025)
N = 1000
A = randn(N,N); A = A*A'/N
x₀ = randn(N); x₀ /= √(x₀'x₀)
μ = 2.0
ε = 1e-12

@time result = inverse_iteration(A,μ,x₀;tol=ε)
result === nothing ? println("Inverse iteration did not converge.") : begin λ,x,n = result ; println("Inverse iteration: eigenvalue $(λ) found in $n iterations.") end
@time result = rayleigh_iteration(A,μ,x₀;tol=ε)
result === nothing ? println("Rayleigh iteration did not converge.") : begin λ_r,x_r,n_r = result ; println("Rayleigh quotient iteration: eigenvalue $(λ_r) found in $n_r iterations.") end

### BEGIN HIDDEN TESTS
if abs(λ-λ_r) < 1/N^2 # A has about N eigenvalues in [0,4]
    @assert (abs(λ-λ_r) ≤ 2ε)
    @assert all(abs.((x'A*x_r)/(x'x_r) .- λ) .≤ ε)
end

@assert n_r ≤ n
### END HIDDEN TESTS

  1.481179 seconds (5.26 M allocations: 270.944 MiB, 2.99% gc time, 93.76% compilation time)
Inverse iteration: eigenvalue 1.998378274695557 found in 27 iterations.

  0.306550 seconds (25.95 k allocations: 108.427 MiB, 44.54% gc time, 17.74% compilation time)

Rayleigh quotient iteration: eigenvalue 1.9873152321450804 found in 7 iterations.

Valeurs propres multiples

Dorénavant, $A$ sera une matrice hermitienne.

[Exercice 2] Itération en sous-espaces et SVD

La méthode d'itération en sous-espaces est définie, étant donné une condition initiale $X_0\in\mathbb{C}^{n\times p}$, par la récurrence : $$ X_{n+1} R_{n+1} = A X_n,$$ où $X_{n+1} R_{n+1}$ est la décomposition QR partielle de $AX_n \in \mathbb{C}^{n\times p}$. On peut montrer, sous des conditions appropriées, que les colonnes de $X_n$ convergent vers les $p$ vecteurs propres dominants de $A$ (voir les notes de cours).

1. Implémenter la décomposition QR partielle, c'est-à-dire une fonction myQR prenant comme argument une matrice $n\times p$ $M$ et renvoyant $Q,R$, où $Q$ est une matrice orthogonale $n\times p$ et $R$ est une matrice triangulaire supérieure $p\times p$ avec des entrées diagonales positives.

   Indice (cliquer pour afficher)

   • Procéder par récurrence sur $p$, en écrivant $$ Q = \begin{pmatrix} \widetilde Q &\mathbf{q}\end{pmatrix},\qquad R = \begin{pmatrix} \widetilde R & \mathbf{r} \\ 0 & \alpha\end{pmatrix}, \qquad M=\begin{pmatrix} \widetilde{M} & \mathbf m\end{pmatrix},$$ avec $\widetilde Q,\widetilde M \in \mathbb{C}^{n\times (p-1)}$, $\widetilde R\in\mathbb{C}^{(p-1)\times (p-1)}$, $\mathbf r \in \mathbb{C}^{p-1}$, $\alpha\in \mathbb{R}_+$ et $\mathbf m\in \mathbb{C}^n$, en supposant que l'on connaît la décomposition $\widetilde Q \widetilde R = \widetilde M$ de la matrice constituée des $(p-1)$ premières colonnes de $M$ (le cas $p=1$ est trivial).

• Pour éviter des allocations mémoire inutiles, il est souvent utile en Julia d'utiliser view(A,ix...) au lieu de A[ix...] lors d'opérations sur des tableaux. Pour comprendre pourquoi, rappelons que Julia alloue par défaut de la mémoire lors du découpage d'un tableau en tant que r-value B = A[ix...] (B est une copie du découpage de A) vs le découpage en tant que l-value. Par conséquent, en Julia :

      A = randn(20000,20000)
      A[1:20000,1] === A[1:20000,1]
      # false
  

  Pour éviter cela, la méthode view(A,ix...) renvoie une référence au découpage de A, qui se comporte par ailleurs comme un tableau. Par exemple :

      @time A[1:20000,1]'A[1,1:20000]
      # 0.000655 seconds (6 allocations: 312.625 KiB)
      @time view(A,1:20000,1)'view(A,1,1:20000)
      # 0.000369 seconds (4 allocations: 208 bytes)
  

  Pour plus de commodité, Julia fournit également la macro @views qui transforme toute expression de la forme A[ix...] en une expression de la forme view(A,ix...) dans le code sur lequel elle opère :

      @views A[1:20000,1] === A[1:20000,1] # transforme une ligne de code
      # true
  
      @time @views A[1:20000,1]'A[1,1:20000]
      # 0.000289 seconds (4 allocations: 208 bytes)
  
      @views function f(M) # ou un bloc de code complet
              return M[1:20000,1]'M[1,1:20000]
          end
      @time f(A)
      # 0.000314 seconds (1 allocation: 16 bytes)
  

@views function myQR(M)
    n,p = size(M)

    @assert p <= n "Error: p > n"

    Q = zero(M)
    R = zeros(eltype(M),p,p)

    ### BEGIN SOLUTION

    for k=1:p
        m = view(M, :, k)
        Q̃ = view(Q, :, 1:k-1)
        r = view(R, 1:k-1, k)
        q = view(Q, :, k)
        
        r .= Q̃'m
        q .= m - Q̃*r
        α = sqrt(q'q)
        q ./= α
        R[k,k] = α
    end

    ### END SOLUTION

    return Q,R
end

myQR (generic function with 1 method)

# Test code
M = randn(ComplexF64,100,50)
M = big.(M)
Q,R = myQR(M)
@show norm(Q*R-M)

### BEGIN HIDDEN TESTS
@assert norm(Q*R-M) < 1e-50
@assert norm(Q'Q-I) < 1e-50
@assert all(isreal.(diag(R)))
@assert all(real.(diag(R)) .>= 0)
### END HIDDEN TESTS

norm(Q * R - M) = 1.036416049574850419425422173720991419736827080356812444774911952379079269252268e-75

2. Implémenter l'itération en sous-espaces, prenant comme arguments la matrice $B$, le nombre de valeurs propres $p$ et un nombre fixe d'itérations $n = $niter. Renvoyer les paires propres approximatives $(\boldsymbol{\lambda}_n,X_n)$ après $n$ itérations.

function myEigen(B, p, niter)

    ### BEGIN SOLUTION
    n = size(B, 1)
    X = randn(n, p)
    R = zeros(p,p)
    for i in 1:niter
        X, R = myQR(B*X)
    end
    λs = diag(X'*B*X)
    ### END SOLUTION

    return λs, X
end

myEigen (generic function with 1 method)

@show myEigen([1. 2.; 2. 1.], 2, 100)[1]
@show myEigen([1. 2.; 2. 1.], 2, 100)[2]

### BEGIN HIDDEN TESTS
@assert begin n = 2; A = randn(n, n); myEigen(A'A, n, 100)[1] ≈ reverse(eigen(A'A).values) end
@assert begin n = 4; A = randn(n, n); myEigen(A'A, n, 100)[1] ≈ reverse(eigen(A'A).values) end
@assert begin A = randn(5, 5); q, r = qr(A); B = q*Diagonal(1:5)*q'; myEigen(B, 3, 100)[1] ≈ [5; 4; 3] end
### END HIDDEN TESTS

(myEigen([1.0 2.0; 2.0 1.0], 2, 100))[1] = [3.0, -1.0]

(myEigen([1.0 2.0; 2.0 1.0], 2, 100))[2] = [-0.7071067811865475 0.7071067811865476; -0.7071067811865475 -0.7071067811865476]

Rappeler la décomposition en valeurs singulières pour $B \in \mathbb{C}^{m\times n}$ $$B = U \Sigma V^*, \quad U\in \mathbb{C}^{m\times m},\quad \Sigma \in \mathbb{R}^{m\times n},\quad V\in \mathbb{C}^{n\times n},$$ où $U U^* = U^* U = I_m$, $V V^* = V^* V = I_n$, et $\Sigma \in \mathbb{R}^{m \times n}$ est une matrice diagonale rectangulaire avec des entrées réelles non négatives sur la diagonale. Les colonnes de $U$ (resp. $V$) sont appelées vecteurs singuliers gauches (resp. droits) de $B$, et les entrées diagonales de $\Sigma$ sont les valeurs singulières (non négatives).

3. Écrire une fonction mySVD(B, p, niter) qui renvoie les p valeurs singulières dominantes de la matrice carrée B (dans un vecteur σs), ainsi que les vecteurs singuliers gauches et droits associés (dans des matrices Up et Vp).

   Indice (cliquer pour afficher)

   • Remarquer que si $A$ est une matrice carrée $$ A A^* = U \Sigma^2 U^*, \qquad A^* A = V \Sigma^2 V^*. $$ Ainsi, les vecteurs singuliers gauches de $A$ sont les vecteurs propres de $A A^*$, tandis que les vecteurs singuliers droits de $A$ sont les vecteurs propres de $A^* A$.

   • Une fois les vecteurs singuliers gauches et droits calculés associés aux p valeurs singulières dominantes, les valeurs singulières elles-mêmes peuvent être obtenues en extrayant la diagonale de la matrice $$ \Sigma_p = U_p^* B V_p. $$ Ici $U_p$ et $V_p$ sont les matrices contenant en colonnes les vecteurs singuliers gauches et droits associés aux p valeurs singulières dominantes, respectivement.

function mySVD(B, p, niter)

    ### BEGIN SOLUTION
    n = size(B, 1)
    λ₁, U = myEigen(B*B', p, niter)
    λ₂, V = myEigen(B'B, p, niter)
    σs = U'B*V
    ### END SOLUTION
    return diag(σs), U, V
end

mySVD (generic function with 1 method)

n = 10
B = randn(n, n)
σs, U, V = mySVD(B, n, 1000)
@assert norm(U'U - I(n)) < 1e-10
@assert norm(V'V - I(n)) < 1e-10
@assert norm(U*Diagonal(σs)*V' - B) < 1e-10

4. La décomposition en valeurs singulières est très utile pour compresser des matrices. L'idée de la compression de matrices basée sur la SVD est la suivante : étant donné $p \leqslant n$, une matrice $B \in \mathbb{C}^{n\times n}$ (on considère des matrices carrées par simplicité) peut être approchée par $$ \widetilde {B} := U_p \Sigma_p V_p^*, $$ où $\Sigma_p \in \mathbb{R}^{p \times p}$ est une matrice diagonale contenant les $p$ valeurs singulières dominantes de $B$ sur sa diagonale, et où $U_p \in \mathbb{C}^{n \times p}$ et $V_p \in \mathbb{C}^{n \times p}$ sont des matrices rectangulaires contenant les vecteurs singuliers gauches et droits associés, respectivement.

   Puisqu'une image en niveaux de gris peut être représentée par une matrice contenant les valeurs d'intensité de tous les pixels, cette approche de compression de matrices peut être utilisée pour compresser des images en niveaux de gris. Utiliser cette méthode, c'est-à-dire calculer $\widetilde {B}$, pour $p \in \{5, 10, 20, 30\}$, afin de compresser l'image woman_darkhair donnée ci-dessous (d'autres images de test disponibles sont listées ici), et tracer l'image compressée pour ces valeurs de p.

   Remarques :

   • (Pour information uniquement) En pratique, au lieu de stocker la matrice complète $\widetilde {B}$, qui contient $n^2$ entrées, on peut stocker uniquement les matrices $U_p$, $\Sigma_p$ et $V_p$, qui contiennent ensemble seulement $(2n+1)p$ entrées (on ne compte que les $p$ entrées diagonales de $\Sigma_p$). Si $p \ll n$, la mémoire nécessaire pour stocker ces matrices est bien inférieure à la mémoire nécessaire pour stocker la matrice initiale $B$.

   • Une fonction pour tracer des images à partir de la matrice des valeurs d'intensité des pixels est fournie ci-dessous, et sert de test pour votre implémentation.

A = testimage("woman_darkhair")

# Convert image to matrix of Float64
M = Float64.(A)
n = size(M,1)

# Function to plot a grayscale image from the matrix
# containing the intensity values of all the pixels
function plot_matrix(B, p)
    plot(Gray.(B), ticks=false, showaxis=false, title="p=$p")
end

plots = typeof(plot())[]

for p in [5,10,20,30]
    niter = 100
    σs, U, V = mySVD(M, p, niter)
    println("p = $p, compression ratio = ", ((2n+1)*p)/(n^2))
    push!(plots,plot_matrix(U*Diagonal(σs)*V', p))
end

push!(plots, plot_matrix(M, "n"))
plot(plots...)

p = 5, compression ratio = 0.019550323486328125

p = 10, compression ratio = 

0.03910064697265625
p = 20, compression ratio = 0.0782012939453125
p = 30, compression ratio = 

0.11730194091796875

Utilisation de matrices creuses pour des problèmes de grande dimension

[Exercice 3] Algorithme PageRank

PageRank est un algorithme qui attribue un score aux sommets d'un graphe orienté. Il était autrefois utilisé par les principaux moteurs de recherche pour classer les résultats de recherche. Dans ce contexte, le graphe orienté code les liens entre les pages du World Wide Web, les sommets représentant les pages web et les arêtes représentant les connexions entre pages : il y a une arête de la page $i$ à la page $j$ si la page $i$ contient un hyperlien vers la page $j$.

Considérons un graphe orienté $G(V, E)$ avec des sommets $V = \{1, \dotsc, n\}$ et des arêtes $E$. Le graphe peut être représenté par sa matrice d'adjacence $A \in \{0, 1\}^{n \times n}$, où les entrées sont données par $$ a_{ij} = \begin{cases} 1 & \text{s'il y a une arête de $i$ à $j$,} \\ 0 & \text{sinon.} \end{cases} $$ L'idée de l'algorithme PageRank, dans sa forme la plus simple, est d'attribuer des scores $r_i$ aux sommets en résolvant le système d'équations suivant : $$ \tag{PageRank} \forall i \in V, \qquad r_i = \sum_{j \in \mathcal{N}(i)} \frac{r_j}{o_j}. $$ où $o_j$ est le degré sortant du sommet $j$, c'est-à-dire le nombre d'arêtes d'origine $j$. Ici, la somme porte sur l'ensemble des nœuds dans $\mathcal{N}(i)$, qui représente l'ensemble des voisins entrants du sommet $i$, c'est-à-dire ceux qui ont une arête pointant vers $i$.

Soit $\mathbf{r} = \begin{pmatrix} r_1 & \dots & r_n \end{pmatrix}^T$. Il est facile de montrer que résoudre le système (PageRank) est équivalent à résoudre le problème suivant : $$ \tag{PageRank-vector} \mathbf{r} = A^T \begin{pmatrix} \frac{1}{o_1} & & \\ & \ddots & \\ & & \frac{1}{o_n} \end{pmatrix} \mathbf{r} =: A^T O^{-1} \mathbf{r}. $$ Autrement dit, le problème se réduit à trouver un vecteur propre de valeur propre $1$ de la matrice $M = A^T O^{-1}$. À ce stade, on n'a prouvé ni l'existence ni l'unicité d'une solution à cette équation. La question de l'unicité d'une solution est liée à la connectivité du graphe et ne sera pas abordée ici. Cependant, nous allons démontrer qu'une solution au problème existe.

Remarque. La matrice $O^{-1} A$ est la matrice de transition d'une marche aléatoire sur le graphe orienté, où à chaque étape on se déplace vers un voisin sortant, avec probabilité égale pour chacun d'eux. Résoudre (PageRank-vector) est équivalent à trouver une distribution stationnaire de cette marche aléatoire.

1. • Remarquer que $M$ est une matrice stochastique à gauche, c'est-à-dire que la somme des éléments de chaque colonne est égale à 1.

   • Montrer que les valeurs propres de toute matrice $B \in \mathbb{R}^{n \times n}$ coïncident avec celles de $B^T$. On pourra utiliser le fait que $\det(B) = \det(B^T)$.

   • En utilisant les points précédents, montrer que $1$ est une valeur propre et que $\rho(M) = 1$. Pour la seconde partie, trouver une norme matricielle subordonnée telle que $\lVert M\rVert= 1$. Cela démontre l'existence d'une solution à (PageRank-vector), et prouve également que $1$ est la valeur propre dominante de $M$.

Nous allons appliquer PageRank pour trier les pages Wikipedia selon leur importance. Les deux cellules suivantes téléchargent et analysent les données dans des tableaux. Pour limiter le temps de calcul, seulement 5 % des articles les mieux notés ont été sélectionnés.

import Downloads
import Tar

# URL where data can be downloaded
url = "https://urbain.vaes.uk/static/wikidata.tar"

# Download the data
filename = "wikidata.tar"
isfile(filename) || Downloads.download(url, filename)

# Extract data into directory `wikidata`
directoryname = "wikidata"
isdir(directoryname) || Tar.extract(filename, directoryname)

true

import CSV
import DataFrames

# Read nodes and edges into data frames
nodes_dataframe = CSV.read("$directoryname/names.csv", DataFrames.DataFrame)
edges_dataframe = CSV.read("$directoryname/edges.csv", DataFrames.DataFrame)

# Convert data to matrices
nodes = Matrix(nodes_dataframe)
edges = Matrix(edges_dataframe)

# The data structures should be self-explanatory
edges_dataframe

10722190×2 DataFrame

10722165 rows omitted

Row	FromNode	ToNode
	Int64	Int64
1	175973	1
2	130880	2
3	145856	2
4	159190	2
5	159200	2
6	159207	2
7	159431	2
8	4	3
9	5	3
10	6887	3
11	6916	3
12	6957	3
13	11490	3
⋮	⋮	⋮
10722179	170697	157057
10722180	177550	45
10722181	30983	29807
10722182	76929	104616
10722183	86508	120348
10722184	128270	198551
10722185	152419	116918
10722186	156536	190926
10722187	170697	112376
10722188	177550	108047
10722189	181302	84697
10722190	189293	154081

Combien y a-t-il de nœuds ? Calculer la mémoire nécessaire pour stocker chaque entrée de $M$ au format Float64.

2. Implémenter une structure struct MySparseMatrix pour représenter une matrice creuse avec des entrées Float64 (au format COO), et une méthode *(M::MySparseMatrix,X::Vector{Float64}) pour calculer le produit d'une matrice creuse M par un vecteur X.

import Base.*

struct MySparseMatrix
    rows::Vector{Int}
    cols::Vector{Int}
    vals::Vector{Float64}
    m::Int
    n::Int
end

MySparseMatrix(R::Vector{Int},C::Vector{Int},V::Vector{Float64}) = MySparseMatrix(R,C,V,maximum(R),maximum(C))

@inbounds function *(M::MySparseMatrix, X::Vector{Float64})
    @assert size(X, 1) == M.n "Incompatible dimensions: M has $(M.n) columns but X has $(length(X)) rows."
    ### BEGIN SOLUTION
    Y = zeros(M.m)
    for (i,j,v) = zip(M.rows, M.cols, M.vals)
        Y[i] += v*X[j]
    end
    return Y
    ### END SOLUTION
end

* (generic function with 285 methods)

Tester votre code avec la cellule ci-dessous :

m, n = 4, 3
R = [2, 2, 2, 3, 3]
C = [1, 2, 3, 1, 3]
V = [5., 6., 7., 8., 9.]
A = MySparseMatrix(R, C, V, m, n)
b = [1.; 1.; 1.]
@assert A*b == [0.; 18.; 17.; 0.] "Multiplication does not work!"

3. Construire la matrice stochastique à gauche $M$ :

nn, ne = length(nodes), size(edges, 1)

### BEGIN SOLUTION
# Count the number of outbound edges for each node
n_outbound = zeros(Int, nn)
for e in eachrow(edges)
    n_outbound[e[1]] += 1
end

# Build matrix
R, C = edges[:, 1], edges[:, 2]
V = 1 ./ n_outbound[R]
### END SOLUTION

M = MySparseMatrix(C, R, V)

MySparseMatrix([1, 2, 2, 2, 2, 2, 2, 3, 3, 3  …  29807, 104616, 120348, 198551, 116918, 190926, 112376, 108047, 84697, 154081], [175973, 130880, 145856, 159190, 159200, 159207, 159431, 4, 5, 6887  …  30983, 76929, 86508, 128270, 152419, 156536, 170697, 177550, 181302, 189293], [0.025, 0.007936507936507936, 0.012658227848101266, 0.0625, 0.017543859649122806, 0.02564102564102564, 0.0037174721189591076, 0.07692307692307693, 0.08333333333333333, 0.0125  …  0.1, 0.5, 1.0, 1.0, 1.0, 0.1, 0.1, 0.1, 1.0, 1.0], 199903, 199903)

4. Implémenter la méthode des puissances pour calculer le vecteur propre associé à la valeur propre principale de $M$. Puisque la valeur propre $\lambda=1$ est connue, on pourra utiliser $\|Mr-r\| \leq \varepsilon\|r\|$ comme critère d'arrêt.

function power_iteration(M, x; ε=1e-12, maxiter=1000)
    ### BEGIN SOLUTION
    niter = 0
    while niter < maxiter
        niter += 1
        x = x / √(x'x)
        Mx = M*x
        λ = x'Mx
        e = norm(Mx - λ*x)
        x = Mx
        e ≤ ε && return x
    end
    return nothing
    ### END SOLUTION
    # Return only the eigenvector, not the eigenvalue
end

power_iteration (generic function with 1 method)

@assert [1, -1]'power_iteration([1. 2.; 2. 1.], [1., 0.]) |> abs < 1e-9
@assert [1, 0]'power_iteration([0. 0.; 0. 1.], [1., 1.]) |> abs < 1e-9
@assert [0, 1]'power_iteration([1. 0.; 0. .5], [1., 1.]) |> abs < 1e-9

La cellule suivante exécute PageRank :

x = ones(nn) / nn
x = @time power_iteration(M, x)

p = sortperm(x, rev=true)
sorted_nodes = view(nodes, p)

print(join(sorted_nodes[1:20], "\n"))
@assert sorted_nodes[1] == "United States"
@assert sorted_nodes[2] == "United Kingdom"
@assert sorted_nodes[3] == "World War II"
@assert sorted_nodes[4] == "Latin"
@assert sorted_nodes[5] == "France"

  6.637422 seconds (8.49 k allocations: 1.126 GiB, 2.43% gc time, 0.27% compilation time)
United States
United Kingdom
World War II
Latin
France
Germany
English language
China
Canada
India
Mathematics
Italy
Catholic Church
Australia
Greek language
Europe
England
World War I
London
Russia

5. Écrire une fonction search(keyword) pour effectuer une recherche dans la base de données. Voici un exemple de ce que cette fonction pourrait renvoyer :

   julia> search("Newton")
   47-element Vector{String}:
    "Isaac Newton"
    "Newton (unit)"
    "Newton's laws of motion"
    …
   

   Indice (cliquer pour afficher)

   • La méthode filter(condition, itr) renvoie une version filtrée de l'itérable itr ne conservant que les éléments x vérifiant condition(x) == true.
   • La méthode occursin(needle::String, haystack::String) renvoie true si et seulement si needle est une sous-chaîne de haystack.

   Bien entendu, la meilleure façon de comprendre ces méthodes est de lire la documentation .

function search(keyword)
    ### BEGIN SOLUTION
    filter(s -> occursin(keyword, s), sorted_nodes)
    ### END SOLUTION
end

search("Newton")

47-element Vector{String}:
 "Isaac Newton"
 "Newton (unit)"
 "Newton's laws of motion"
 "Newton's law of universal gravitation"
 "Newton's method"
 "Newtonian fluid"
 "Non-Newtonian fluid"
 "Newton metre"
 "Newton, Massachusetts"
 "Olivia Newton-John"
 "Modified Newtonian dynamics"
 "Newtonian telescope"
 "Newton's notation"
 ⋮
 "Newtonmore"
 "Gaussâu80u93Newton algorithm"
 "James Newton Howard"
 "Robert Newton"
 "John Gilbert Newton Brown"
 "Newton's inequalities"
 "Newton's method in optimization"
 "Thandie Newton"
 "Wayne Newton"
 "Post-Newtonian expansion"
 "Quasi-Newton method"
 "Helmut Newton"