{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Itération simultanée pour le calcul de valeurs propres\n",
"\n",
"On veut implémenter une méthode itérative pour calculer simultanément les $p \\leq n$ valeurs propres dominantes d'une matrice $A\\in \\mathbb C^{n\\times n}$, qu'on supposera Hermitienne.\n",
"La méthode s'appuie sur la décomposition QR réduite, de la forme:\n",
"\n",
"$$ QR = M,$$\n",
"où:\n",
"- $M\\in\\mathbb C^{n\\times p}$.\n",
"- $Q\\in \\mathbb C^{n\\times p}$ et ses colonnes forment une famille orthonormée de $\\mathbb C^n$ pour le produit scalaire usuel.\n",
"- $R\\in \\mathbb C^{p\\times p}$ est une matrice triangulaire supérieure avec des entrées diagonales réelles et positives.\n",
"\n",
"#### Question 1\n",
"Donner la décomposition QR réduite de $M$ si $p=1$.\n",
"\n",
"#### Question 2\n",
"En écrivant\n",
"$$ 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}$$\n",
"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 \\R_+$ et $\\mathbf m\\in \\mathbb C^n$, donner la décomposition QR de $M$ sachant la décomposition $\\widetilde Q \\widetilde R = \\widetilde M$ de la matrice formée des $(p-1)$ premières colonnes de $M$. (Utiliser les propriétés de la décomposition QR).\n",
"\n",
"#### Question 3\n",
" Implémenter une méthode pour calculer la décomposition QR réduite.\n",
"\n",
"**Remarque.**\n",
"Par souci d'efficacité,\n",
"on peut utiliser la macro `@views`,\n",
"qui permet d'effectuer des opérations avec des sous-matrices sans en copier le contenu. (voir `?view` et `?@views` dans le REPL).\n",
"Voir [cette partie](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-views) de la documentation pour plus d'informations.\n",
"Nous reprenons ci-dessous l'exemple de la page du manuel pour illustrer ce point."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"using BenchmarkTools\n",
"x = rand(10^6)\n",
"fcopy(x) = sum(x[2:end-1])\n",
"@views fview(x) = sum(x[2:end-1])\n",
"t1 = @belapsed fcopy(x)\n",
"t2 = @belapsed fview(x)\n",
"println(t1, \" \", t2)\n",
"\n",
"# The following code may be helpful to develop your understanding\n",
"a = zeros(3)\n",
"@views b = a[1:2] # Does not copy a[1:2]\n",
"c = a[2:3] # Does copy a[2:3]\n",
"a[2] = 1.\n",
"println(b, \" \", c)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "07cea03a",
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra\n",
"\n",
"@views function myQR(M)\n",
" n,p = size(M)\n",
"\n",
" @assert p <= n \"Error: p > n\"\n",
" ## implémenter la décomposition QR ici\n",
" Q = zero(M)\n",
" R = zeros(eltype(M),p,p)\n",
"\n",
" return Q,R\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"M = randn(ComplexF64,100,50)\n",
"M = Complex{BigFloat}.(M)\n",
"Q,R = myQR(M)\n",
"\n",
"@show norm(Q'Q-I) # vérifier que Q a bien des colonnes orthonormales\n",
"@show all(isreal.(diag(R))) # R a des coefficients diagonaux réels\n",
"@show all(real.(diag(R)) .>= 0) # et positifs\n",
"@show norm(Q*R-M) # et que QR = M à une erreur d'arrondi près\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@time qr(M) #implémentation de LinearAlgebra\n",
"@time myQR(M) #implémentation maison"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Question 4\n",
"La méthode d'itération simultanée est définie, étant donné une condition initiale $X_0\\in\\mathbb C^{n\\times p}$, par la récurrence:\n",
"$$ X_{n+1} R_{n+1} = A X_n,$$\n",
"où $X_{n+1} R_{n+1}$ est la décomposition QR réduite de $AX_n \\in \\mathbb C^{n\\times p}$. Il est possible de montrer que pour $A$ Hermitienne, et sous certaines conditions initiales $X_0$ et les $p$ premières valeurs propres de $A$ de plus grand module, $X_n$ converge essentiellement colonne par colonne vers la matrice formée par les $p$ vecteurs propres associés (voir les notes de cours pour une preuve).\n",
"\n",
"Implémenter une fonction `myEigen` en Julia prenant en arguments : une matrice `M` (supposée Hermitienne), un entier `nev` (le nombre de valeurs propres à calculer), un entier `niter` (le nombre d'itérations), et une matrice initiale `X₀`. La fonction devrait renvoyer un vecteur de `nev` valeurs propres et une matrice formée d'une approximation des `nev` vecteurs propres dominants de `M` en colonnes.\n",
"\n",
"Pour des matrices dont les entrées ne sont pas des `BigFloat`, on pourra comparer les résultats ceux donnés par `LinearAlgebra.eigen` pour débugger."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n,nev = 10,3\n",
"\n",
"## seed pour le RNG\n",
"using Random\n",
"Random.seed!(2023)\n",
"##\n",
"\n",
"X = randn(n,n)/sqrt(n)\n",
"X₀ = randn(n,nev)/sqrt(n)\n",
"\n",
"\n",
"function myEigen(M,nev,niter,X₀)\n",
" ## implémenter l'itération ici\n",
" return λs,us\n",
"end\n",
"\n",
"M = BigFloat.(I + X'X)\n",
"\n",
"@time λs,us=myEigen(M,nev,500,X₀)\n",
"\n",
"display(λs)\n",
"display(us)\n",
"display(M*us - λs' .*us)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Question 5\n",
"\n",
"Étudier numériquement la vitesse de convergence des 3 premiers vecteurs propres vers leur limite en fonction du nombre d'itérations de la méthode. On pourra considérer que tous les vecteurs propres sont convergés après un grand nombre (disons 10000) itérations pour obtenir des valeurs de référence. L'écart entre deux vecteurs $u$ et $v$ sera mesuré en terme de l'angle formé par les sous-espaces engendrés:\n",
"$$ \\theta_{u,v} = \\arccos\\left(\\frac{|u^\\intercal v|}{\\|u\\|\\|v\\|}\\right).$$\n",
"Comparer l'erreur commise à la $k$-ème itération et les quantités $(\\lambda_i/\\lambda_{i+1})^{-k},\\, i=1,2,3$, sur un graphe avec une échelle log en ordonnée.\n",
"\n",
"Pour garantir numériquement que $0\\leq \\frac{|u^\\intercal v|}{\\|u\\|\\|v\\|}\\leq 1$, on pourra utiliser la fonction `clamp`."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@time λs,us=myEigen(M,n,10000,randn(n,n)) # on calcule les valeurs et vecteurs propres \"convergés\"\n",
"display(M*us - λs' .*us)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"nev=3\n",
"\n",
"iters = 1:200\n",
"errs = [Float64[] for i=1:nev]\n",
"\n",
"for niter in iters\n",
" λhats,uhats = myEigen(M,nev,niter,X₀)\n",
" ## calculer les erreurs ici\n",
"end\n",
"\n",
"using Plots,LaTeXStrings\n",
"colors = (:red,:blue,:green)\n",
"pl =plot(xlabel=\"Itération k\", ylabel=\"Erreur essentielle\",yaxis=:log,legend=:bottomleft) # on déclare un plot\n",
"\n",
"for i=1:nev\n",
" plot!(pl,iters,errs[i],label=L\"θ_{%$i}\",color=colors[i])\n",
" plot!(pl,k->(λs[i]/λs[i+1])^-k,label=L\"(λ_{%$i}/λ_{%$(i+1)})^{-k}\",linestyle=:dot,color=colors[i])\n",
"end\n",
"\n",
"pl"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "0f7dc8eb",
"metadata": {},
"source": [
"# Algorithme PageRank"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "2530e31e",
"metadata": {},
"source": [
"PageRank est un algorithme permettant d'attribuer un score aux sommets d'un graphe orienté.\n",
"Il est utilisé par de nombreux moteurs de recherche pour trier les résultats de recherche.\n",
"Dans ce contexte, le graphe orienté encode les liens entre les pages du World Wide Web :\n",
"les sommets du graphe orienté représentent des pages web,\n",
"et les arêtes représentent les connexions entre les pages:\n",
"il y a une arête allant de la page $i$ à la page $j$ si la page $i$ contient un lien hypertexte vers la page $j$.\n",
"\n",
"Considérons un graphe orienté $G(V, E)$ avec des sommets $V = \\{1, \\dotsc, n\\}$ et des arêtes $E$.\n",
"Le graphe peut être représenté par sa matrice d'adjacence $\\mathsf A \\in \\{0, 1\\}^{n \\times n}$,\n",
"dont les entrées sont données par\n",
"$$\n",
"a_{ij} =\n",
"\\begin{cases}\n",
" 1 & \\text{si il y a une arête de $i$ à $j$,} \\\\\n",
" 0 & \\text{sinon.}\n",
"\\end{cases}\n",
"$$\n",
"L'idée de l'algorithme PageRank, dans sa forme la plus simple,\n",
"consiste à attribuer des scores $r_i$ aux sommets\n",
"en résolvant le système d'équations suivant :\n",
"$$ \\tag{PageRank}\n",
" \\forall i \\in V, \\qquad\n",
" r_i\n",
" = \\sum_{j \\in \\mathcal N(i)} \\frac{r_j}{o_j}.\n",
"$$\n",
"\n",
"où $o_j$ est le degré sortant du sommet $j$,\n",
"c'est-à-dire le nombre d'arêtes ayant $j$ pour origine.\n",
"Ici, la somme s'applique à l'ensemble des nœuds dans $\\mathcal N(i)$,\n",
"qui représente l'ensemble des voisins entrants du sommet $i$,\n",
"c'est-à-dire ceux qui ont une arête pointant vers le sommet $i$.\n",
"\n",
"Soit $\\mathbf r = \\begin{pmatrix} r_1 & \\dots & r_n \\end{pmatrix}^T$.\n",
"Il est simple de montrer que résoudre le système (PageRank) revient à résoudre le problème suivant:\n",
"$$ \\tag{PageRank-vector}\n",
" \\mathbf r =\n",
" \\mathsf A^T\n",
" \\begin{pmatrix}\n",
" \\frac{1}{o_1} & & \\\\\n",
" & \\ddots & \\\\\n",
" & & \\frac{1}{o_n}\n",
" \\end{pmatrix}\n",
" \\mathbf r =: \\mathsf A^T \\mathsf O^{-1} \\mathbf r.\n",
"$$\n",
"\n",
"En d'autres termes,\n",
"le problème revient à trouver un vecteur propre de valeur propre $1$ de la matrice $\\mathsf M = \\mathsf A^T \\mathsf O^{-1}$.\n",
"Notons qu'à ce stade, nous n'avons prouvé ni l'existence,\n",
"ni l'unicité d'une solution de cette équation.\n",
"La question de l'unicité d'une solution est liée à la connectivité du graphe et ne sera pas abordée ici.\n",
"Nous allons démontrer, par contre, qu'il existe une solution au problème.\n",
"\n",
"**Remarque.** La matrice $\\mathsf O^{-1} \\mathsf A$ est la matrice de transition d'une marche aléatoire sur le graphe orienté,\n",
"où à chaque pas un déplacement se fait vers un voisin sortant,\n",
"avec une probabilité égale pour chacun d'eux.\n",
"Résoudre (PageRank-vector) revient à trouver une distribution stationaire de cette marche aléatoire.\n",
"\n",
"**Questions:**\n",
"\n",
"1. Remarquer que $\\mathsf M$ est une matrice stochastique gauche,\n",
" c'est-à-dire que la somme des éléments de chaque colonne est égale à 1.\n",
"\n",
"1. Prouver que les valeurs propres de n'importe quelle matrice $\\mathsf B \\in \\mathbb R^{n \\times n}$ coïncident avec celles de $\\mathsf B^T$.\n",
" Vous pouvez utiliser le fait que $\\det(\\mathsf B) = \\det(\\mathsf B^T)$.\n",
"\n",
"1. En utilisant les points précédents, montrer que 1 est une valeur propre et que $\\rho(\\mathsf M) = 1$.\n",
" Pour la deuxième partie, trouver une norme matricielle subordonnée telle que $\\lVert\\mathsf M\\rVert= 1$.\n",
" Ceci démontre l'existence d'une solution de (PageRank-vector),\n",
" et prouve aussi que 1 est la valeur propre dominante de $\\mathsf M$\n",
"\n",
"1. Implémenter sans utiliser le package `SparseArrays` l'algorithme PageRank afin de classer les pages de l'encyclopédie Wikipedia anglophone telle qu'elle était en 2013\n",
" (voir indication ci-dessous pour le téléchargement du jeu de données).\n",
" Après avoir attribué un score à toutes les pages,\n",
" imprimer les 10 pages les mieux classées,\n",
" ce qui devrait donner ceci:\n",
" ```\n",
" United States, United Kingdom, World War II, Latin, France, Germany, English language, China, Canada, India\n",
" ```\n",
"\n",
"1. Écrire une fonction `search(keyword)` permettant d'effectuer une recherche dans la base de données.\n",
" Voici un exemple de ce que cette fonction pourrait renvoyer :\n",
" ```\n",
" julia> search(\"Newton\")\n",
" 47-element Vector{String}:\n",
" \"Isaac Newton\"\n",
" \"Newton (unit)\"\n",
" \"Newton's laws of motion\"\n",
" …\n",
" ```"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "f870c2a9",
"metadata": {},
"source": [
"**Indications supplémentaires**:\n",
"\n",
"- Le code ci-dessous permet de télécharger le jeu de données."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d9602dc0",
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra\n",
"import Downloads\n",
"import Tar\n",
"\n",
"# URL where data can be downloaded\n",
"url = \"https://urbain.vaes.uk/static/wikidata.tar\"\n",
"\n",
"# Download the data\n",
"filename = \"wikidata.tar\"\n",
"isfile(filename) || Downloads.download(url, filename)\n",
"\n",
"# Extract data into directory `wikidata`\n",
"dirname = \"wikidata\"\n",
"isdir(dirname) || Tar.extract(filename, \"wikidata\")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "f6124b73",
"metadata": {},
"source": [
"Ce jeu de données contient un sous-ensemble des données disponibles publiquement ici.\n",
"Afin de limiter le temps de calcul, seuls 5% des articles les mieux notés ont été conservés pour cet exercice.\n",
"Une fois l'archive décompressée,\n",
"le jeu de données peut être chargé dans Julia de la manière suivante :"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b3025b8a",
"metadata": {},
"outputs": [],
"source": [
"import CSV\n",
"import DataFrames\n",
"\n",
"# Read nodes and edges into data frames\n",
"nodes_dataframe = CSV.read(\"$dirname/names.csv\", DataFrames.DataFrame)\n",
"edges_dataframe = CSV.read(\"$dirname/edges.csv\", DataFrames.DataFrame)\n",
"\n",
"# Convert data to matrices\n",
"nodes = Matrix(nodes_dataframe)\n",
"edges = Matrix(edges_dataframe)\n",
"\n",
"# The data structures should be self-explanatory\n",
"edges_dataframe"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "a5eee91d",
"metadata": {},
"source": [
"- Vous pouvez utiliser soit la version simplifiée de l'algorithme donnée en (PageRank-vector),\n",
" soit la version améliorée avec un facteur d'amortissement décrite sur Wikipedia.\n",
"\n",
"- Comme la valeur propre dominante est connue *a priori*,\n",
" le critère d'arrêt suivant pourra être utilisé:\n",
" $$\n",
" \\frac{\\lVert\\mathsf M \\widehat{\\mathbf r} - \\widehat{\\mathbf r}\\rVert_1}{\\lVert\\widehat{\\mathbf r}\\rVert_1} < 10^{-10}.\n",
" $$\n",
" où $\\widehat{\\mathbf r}$ est une approximation du vecteur propre correspondant à la valeur propre dominante.\n",
"\n",
"- Toutes les méthodes de résolution de (PageRank-vector) nécessitent de calculer des produits matrice-vecteur avec la matrice $\\mathsf M$.\n",
" Cette matrice étant très grande,\n",
" elle ne peut pas être stockée comme une matrice dense de type `Matrix`.\n",
" Pour résoudre le problème, il est donc nécessaire de définir une structure de matrice creuse.\n",
" On pourra par exemple utiliser le format de matrice creuse COO,\n",
" qui est le plus simple.\n",
" Une autre option est d'utiliser le format CSR,\n",
" qui permet de calculer les produits matrice-vecteur de manière plus efficace.\n",
" Il sera utile, pour résoudre le problème aux valeurs propres,\n",
"\n",
" * de définir la méthode de multiplication matrice-vecteur `*(A::my_sparse_matrix, b::Vector{Float64})`.\n",
"\n",
" * de définir une fonction `my_sparse` permettant de construire une matrice creuse à partir de vecteurs `rows`, `cols` et `vals` au format COO.\n",
" Évidemment, si vous choisissez d'utiliser le format COO,\n",
" cette fonction est triviale."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "062e7d05",
"metadata": {},
"outputs": [],
"source": [
"# Necessary to overload `*`\n",
"import Base: *\n",
"\n",
"# Modify if necessary\n",
"struct my_sparse_matrix\n",
" rows::Vector{Int}\n",
" cols::Vector{Int}\n",
" vals::Vector{Float64}\n",
" m::Int\n",
" n::Int\n",
"end\n",
"\n",
"function *(A::my_sparse_matrix, b::Vector{Float64})\n",
" # Fill me\n",
"end\n",
"\n",
"# Constructor from parameters in COO format (modify if necessary)\n",
"function my_sparse(rows, cols, vals, m, n)\n",
" my_sparse_matrix(rows, cols, vals, m, n)\n",
"end\n",
"\n",
"# Test the code\n",
"m, n = 4, 3\n",
"R = [2, 2, 2, 3, 3]\n",
"C = [1, 2, 3, 1, 3]\n",
"V = [5., 6., 7., 8., 9.]\n",
"A = my_sparse(R, C, V, m, n)\n",
"b = [1.; 1.; 1.]\n",
"@assert A*b == [0.; 18.; 17.; 0.] \"Multiplication does not work!\""
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "e76caa60",
"metadata": {},
"source": [
"**Remarque.** Bien qu'il soit demandé de ne pas utiliser `SparseArrays` dans votre code final,\n",
"l'utilisation de ce package pour vos tests est encouragée.\n",
"Il peut être utile, par exemple, d'implémenter PageRank d'abord avec une matrice creuse construite avec `SparseArrays`,\n",
"avant d'utiliser votre propre structure de matrice creuse.\n",
"Pour rappel, une matrice creuse peut être construite à partir des paramètres du format COO de la manière suivante avec `SparseArrays`:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bdd7304e",
"metadata": {},
"outputs": [],
"source": [
"import SparseArrays\n",
"A = SparseArrays.sparse(R, C, V, m, n)\n",
"@assert A*b == [0.; 18.; 17.; 0.] \"Multiplication does not work!\""
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "1aa0c435",
"metadata": {},
"source": [
"Nous vous proposons de découper le reste de votre code en trois cellulles:\n",
"\n",
"- Construire la matrice $\\mathsf M$:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "709cdf94",
"metadata": {},
"outputs": [],
"source": [
"# Fill me"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "dec46ca7",
"metadata": {},
"source": [
"- Résoudre le problème PageRank (la fonction `sortperm` peut être utile ici) et imprimer les noms des 10 pages les mieux classées:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5a9c8322",
"metadata": {},
"outputs": [],
"source": [
"# Fill me"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "021d3952",
"metadata": {},
"source": [
"- Définir et tester la fonction `search`\n",
"(les fonctions `filter` et `occursin` peuvent être utiles ici):"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "595b9f95",
"metadata": {},
"outputs": [],
"source": [
"function search(keyword)\n",
" # Fill me\n",
"end\n",
"search(\"Newton\")"
]
}
],
"metadata": {
"jupytext": {
"encoding": "# -*- coding: utf-8 -*-"
},
"kernelspec": {
"display_name": "Julia 1.9.0",
"language": "julia",
"name": "julia-1.9"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.9.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}