‎

1. Détails numériques des expérimentations

1 Détails numériques des expérimentations

On se propose ici de décrire comment les expérimentations numériques furent réalisées. Le répertoire writings comporte les '.org' des expérimentations. Le code python est 'tangle' depuis ces fichiers vers le répertoire DD. Une première section listera les différents fichiers, en précisant leurs dépendances. Dans une seconde section, l'objectif/contenu des fichiers sera donné. Enfin une troisième section décrira numériquement les expérimentations. On rappelle avant que l'objectif de ces expérimentations est l'étude numérique du gradient conjugué déflaté, pour les résolutions successives de seconds membres multiples. Le système résoudre est de la forme: Ax^s=b^s avec s=1,2,…,ν et tel que b^s+1 soit fonction de la solution x^s calculée à l'étape s. On note n la dimension de A. Le gradient conjugué déflaté noté CG-DEF cherche la solution x^s à l'étape s dans l'espace K^s_km(A,W^s-1,r₀^s=vect(W^s-1,V^s_m), où m est le nombre d'itérations du CG-DEF, et r₀^s=b-Ax₀^s le résidu à l'itération initiale de l'étape s. On note aussi V_m^s=[v₀^s|v₁^s|…|v_m-1^s] la matrice de taille n× n dont les vecteurs colonnes forment une base pour le sous-espace de Krylov K_m^s(A,r₀^s). Et W^s=[w₁|w₂|…|w_k] la matrice de déflation (de taille \(n\times k\)) de déflation dont les vecteurs colonnes sont k approximations des vecteurs propres de A.

1.1 Listes des fichiers

Dans cette section, on donne la liste des fichiers '.py' situés dans le répertoire DD. On précise aussi leurs dépendances.

1.1.1 Fichier pour les solveurs

Le fichier DDMPY.py contient les solveurs pour les problèmes de la forme Ax=b et Aw=λ w. On donne ici la liste librairies python utilisées:

__future__
- import print_function
collections
- import OrderedDict
traceback
time
os
- utilisée pour lancer les script '.R' pour les sorties graphiques
scipy
- import io
scipy.sparse = import [coo|csr|csc]_matrix
scipy.linalg
scipy.sparse.linalg
- import inv et eigsh
matplotlib.pyplot
csv
- import * utilise la fonction csv_cg() pour sauvegarder les données du cg()

La fonction cg() de DDMPY.py utilise la fonction csv_cg(). Cette fonction est tangle depuis writings/MATRIX4EXP.org vers DD/csv.py.

1.1.2 Fichiers Pour les expériementations numériques

Ces fichiers sont tangles depuis le répertoire writings vers le répertoires DD. Ils ont pour objectifs de réaliser les expérimentations numériques du rapport ../report/PFE_LEDERER.pdf. Les dépendances avec les librairies sont:

DDMPY.py pour utiliser le gradient conjugué déflaté via la classe ConjGrad() et la fonction cg() et la fonction setmatrix4exp() (tangle depuis writings/MATRIX4EXP.org vers DD/csv.py) pour

la construction des matrices cas tests et la fonction csv_eigen() pour sauvegarder dans un fichier la qualité de l'information spectrale (tangle depuis writings/MATRIX4EXP.org vers DD/csv.py)

os/sys pour executer les script '.R' en system command lines

La liste des fichiers d'expérimentations est la suivante:

EXP1_1.py Pour étudier la convergence du CG-DEF lorsque les vecteurs colonnes de la matrice de déflation W sont les vecteurs propres exactes de \(A\)
- EXP1_1.R et REXP1_1.R pour la gestion des données sauvegarder et les sorties graphiques via ggplot2
EXP1_2.py Pour étudier la convergence du CG-DEF lorsque les vecteurs colonnes de la matrice de déflation W sont les vecteurs propres perturbés de \(A\)
- EXP1_2.R pour la gestion des données sauvegarder et les sorties graphiques via ggplot2
EXP2_1.py Pour étudier la convergence la qualité de l'information spectral généré par la méthode de Lanczos Rayleigh-Ritz (RR) et Harmonic-Ritz (HR)
- utilise la fonction lanczos_RR() tangle depuis writings/MATRIX4EXP.org vers DD/csv.py pour exécuter la procédure de Lanczos RR ou HR
- REXP2_1.R pour la gestion des données sauvegarder et les sorties graphiques via ggplot2
EXP3_1.py pour étudier la convergence du CG-DEF lorsque ν=2 et avec W^(s=2) calculé via l'approche basée RR ou HR
- EXP3_1.R, REXP3_1.R et REXP3_1_L.R pour la gestion des données sauvegarder et les sorties graphiques via ggplot2
EXP3_2.py pour étudier la convergence du CG-DEF lors des résolutions successive à seconds membres multiples (ν ≥ 2) et avec W^(s) calculé via l'approche basée RR ou HR
- EXP3_2.R, REXP3_2.R et REXP3_2_L.R pour la gestion des données sauvegarder et les sorties graphiques via ggplot2

1.1.3 Fichiers pour les sorties graphiques

Le fichier csv.py comporte les fonctions csv_cg(), csv_eigen(), lanczos_RR() et setmatrix4exp() qui sont tangles depuis writings/MATRIX4EXP.org vers DD/csv.py. On rappelle que csv.py est une dépendance de DDMPY.py=et que les dépendances de =csv.py sont numpy, scipy.linalg et os.

csv_cg() est utilisée dans la fonction cg() pour sauvegarder la convergence du CG-DEF au cours des itérations
csv_eigen() est utilisée dans les fichiers EXP pour sauvegarder la qualité de l'information spectrale
lanczos_RR() pour la méthode de Lanczos RR et HR
setmatrix4exp() pour construire les matrices cas test

Les fichiers '.R' permettent de gérer les données issues de csv_cg() sauvegardées dans les fichiers et de générer les sorties graphiques (.pdf et .png) dans les répertoires nommés en fonctions des paramètres de l'expérimentation (n,k, RR, HR, EXP[1_1|1_2|3_1|3_2]) par exemple ../DD/HR_n_100_k_4_EXP3_2 pour les sorties du CG-DEF. Dans ces répertoires, l'utilisateur va trouver deux sous-répertoires EIGEN et ITE. Le répertoire EIGEN pour sauvegarder la qualité de l'information spectrale et le répertoire ITE pour sauvegarder la convergence du CG-DEF. Ces répertoires contiennent les sous-répertoires nommés en fonctions de la matrice utilisée A2 comme exemple. Les sorties graphiques sont réalisées à ce niveau. Les fichers '.pdf' et '.png' sont nommés en fonctions des paramètres du cas test:

BA pour backward error, DE pour direct error et IBA pour iterated backward error
[A0|A1|A2|A3|A4] suivant la matrice cas test utilisée
L_l1_l2 avec [l1,l2] l'intervalle sur lequel le CG-DEF sauvegarde les coefficients du CG-DEF pour calculer W (exemple BA_A2_L₀₅₁)
ZOOM une sortie graphique sur les premières itérations pour mieux observer l'effet de l'approximation de rang faible
[LD|MD|LMD] des sorties graphiques via des facet_grid plus poussées suivant que la déflation soit à partir de LD, MD ou LMD (exemple LD1BA)

Les sorties liées à l'expérimentation EXP2_1.py (étude de la procédure de Lanczos RR et HR) sont dans le répertoire LANCZOS.

1.2 Objectifs des fonctions

Cette sections présentent les fonctions utilisées pour le CG-DEF et les expérimentations.

1.2.1 Dans `DDMPY.py`

On présente ici les deux fonctions permettant de gérer la déflation, ainsi que les nouvelles entrées/sorties que la déflation apporte a la fonction cg() et à la classe ConjGrad().

Gestion des paramètres de déflation

Pour faire varier les paramètres ℓ et k de la déflation utilise plusieurs cas tests:

SI

['S0', 'S1', 'S2', 'S3'] liste de strings permettant de gérer les paramètres l et k de la déflation

'S0' pour utiliser le CG (pas de déflation) et 'S1', 'S2' ou 'S3' pour utiliser lz CG-DEF
'S1' pour paramétrer l et k au nombre d'itérations réalisé lors de la résolution précédente
'S2' pour paramétrer l like 'S1' avec k dans [1, 2, 3, 4, 10, 50,…,n] (par exemple, (custom list))
'S3' pour fixer explicitement l et k

En pratique pour ℓ on utilise deux entiers l1 et l2 tel que si le CG-DEF converge en n_iter itérations, on sauvegarde les coefficients du CG-DEF calculés entre les itérations l1 et l2 dans les vecteurs/matrices globaux _cg_gamma, _cg_ommega, _cg_d, _cg_p, _cg_z, _cg_mu . La fonction apply_strat(strat, W, n_iter) (dans DDMPY.py) s'assure que =SI soit bien valide en fonction de n_iter et reshape en conséquence les vecteurs /matrices. Les entrées/sorties de la fonction sont les suivantes:

strat : list object
      [Wgiven, k, l1, l2, 'S0/S1/S2/S3/S4', 'LD/MD,LMD', s, 'RR/HR', sigma]
      - Wgiven : bolean, set to True  if the user uses its own deflation subspace.
      - k : integer, number of eigenvector to compute.
      - [l1, l2] : integer, used to reduce information from the cg iterations.
      - str['S0|S1|S2|S3|S4'] : used to indicate how to reduce :
	    'S0': no deflation, 'S1': l1=0|l2=k=n_iter, 'S2': l1=0|l2=n_iter
	    'S3': user parameters default
      - str['LD|MD|LMD'] : used to indicate which part of the spectrum need to be
	    computed. 'LD' Least Dominant, 'MD' Most Dominant, 'LMD' From both sides.
      - s : integer, current system to solve.
      - str['RR|HR'] : if 'RR' use Lanczos ritz procedure to compute eigenpair
	    if 'HR' use Harmonic ritz.
      - sigma : real, Find eigenvalues near sigma.
W : Matrix-like object, the deflation subspace.
n_iter : integer, number of iteration computed by cg to solve the s-1 system.

Explicit Returns : strat.
Implicit Returns : _cg_omega, _cg_gamma, _cg_d, _cg_z, _cg_p, _cg_mu
		  global variables.

Construction et résolution du problème aux valeurs propres RR/HR

On réutilise les coefficients calculés lors de la résolution précédente pour calculer les approximations W des vecteurs propres de A. Pour ce faire on utilise la fonction lanczos_hr(strat, n_iter, A, M, W) avec les entrées/sorties suivantes:

strat : list object
      [Wgiven, k, l1, l2, 'S0/S1/S2/S3/S4', 'LD/MD,LMD', s, 'RR/HR', sigma]
      - Wgiven : bolean, set to True  if the user uses its own deflation subspace.
      - k : integer, number of eigenvector to compute.
      - [l1, l2] : integer, used to reduce information from the cg iterations.
      - str['S0|S1|S2|S3|S4'] : used to indicate how to reduce :
	    'S0': no deflation, 'S1': l1=0|l2=k=n_iter, 'S2': l1=0|l2=n_iter
	    'S3': user parameters default
      - str['LD|MD|LMD'] : used to indicate which part of the spectrum need to be
	    computed. 'LD' Least Dominant, 'MD' Most Dominant, 'LMD' From both sides.
      - s : integer, current system to solve.
      - str['RR|HR'] : if 'RR' use Lanczos ritz procedure to compute eigenpair
	    if 'HR' use Harmonic ritz.
      - sigma : real, Find eigenvalues near sigma.
W : Matrix-like object, the deflation subspace.
A : Matrix-like object from Ax=b
M : Matrix-like object, a preconditionner 
n_iter : integer, number of iteration computed by cg to solve the s-1 system.

Explicit Returns : W the updated deflation matrix

Pour calculer W on utilise l'approche basée RR ou HR suivant la valeur de strat[7]. La liste strat permet de gérer les paramètres de déflation et les cas test.

Nouvelles entrées/sorties de la fonction cg() et de la classe ConjGrad()

Pour le cg() on a rajouté les entrées suivantes:

- W : matrix-like object, the deflation matrix. Default is None.
- strat : list object
      [Wgiven, k, l1, l2, 'S0/S1/S2/S3/S4', 'LD/MD,LMD', s, 'RR/HR', sigma]
      - Wgiven : bolean, set to True  if the user uses its own deflation subspace.
      - k : integer, number of eigenvector to compute.
      - [l1, l2] : integer, used to reduce information from the cg iterations.
      - str['S0|S1|S2|S3|S4'] : used to indicate how to reduce :
	    'S0': no deflation, 'S1': l1=0|l2=k=n_iter, 'S2': l1=0|l2=n_iter
	    'S3': user parameters default
      - str['LD|MD|LMD'] : used to indicate which part of the spectrum need to be
	    computed. 'LD' Least Dominant, 'MD' Most Dominant, 'LMD' From both sides.
      - s : integer, current system to solve.
      - str['RR|HR'] : if 'RR' use Lanczos ritz procedure to compute eigenpair
	    if 'HR' use Harmonic ritz.
      - sigma : real, Find eigenvalues near sigma.
-tags : list object ['A0/A1/A2/A3/A4', str(Irun), Stheta, rep, irr, str(enum_c), ISI, Ipart]
      - tags[0] in str['A0/A1/A2/A3/A4'] used to indicates wich matrix is used
      - tags[1] = str(Irun) where Irun is the run's number
      - tags[2] = Stheta a string in ['1E2/1E4/1E6'] used to indicates the eigenvalue orderr of the LD, MD or LMD
      - tags[3] = rep is the repository name for outputs
      - tags[4] = irr = "1" or "2" used to write header only one time
      - tags[5] = str(enum_c) where enum_c is an integer used to count cg() execution
      - tags[6] = ISI, where ISI is a string in str['S0|S1|S2|S3|S4']
      - tags[7] = Ipart, where Ipart is a string in str['LD|MD|LMD']

La liste tags est utilisée pour les sorties graphiques, notamment pour nommer les labels, titres, répertoires, ect des sorties. La fonction cg() retourne initialement la solution calculée x et le nombre d'itération réalisés noté i dans le code. Lorsque l'on utilise la déflation, la fonction cg() retourne en plus la matrice W. Concernant la classe ConjGrad(), elle prend aussi en entrée W, strat et tags. En sortie, si strat[4]=='S0', elle retourne uniquement la solution x (pas de déflation CG), sinon elle retourne en plus W (déflation CG-DEF).

1.3 Description des expériences numériques.

Pour étudier la convergence du CG-DEF suivant les nombreux paramètres, on utilise l'instruction for suivant ces paramètres: n et nk sont utilisés respectivement pour fixer la dimension de la matrice A et le nombre de valeurs propres clusterisées. Les valeurs par défaut sont n=500 et nk =1. Elles sont paramètrables à l'exécution :

python3 EXP.py n nk

Les autres paramètres sont: Other parameters are:

mat

["A0", "A1", "A2", "A3", "A4"] Une liste de string pour paramétrer par nom l'opérateur linéaire de Ax=b

le constructeur est setmatrix4exp(which=Imat, n=n, k=nk, theta=Itheta)
Imat l'itérand pour la liste mat

theta

[1e2, 1e4, 1e6] une liste de réel pour paramétrer LD (1/theta) ou MD (theta) eigenvalue order.

Itheta l'itérand de theta

maxrun

un entier, fixé à 3 par défaut, donnant le nombre de run. Doit aussi être modifé dans le '.R'

Irun l'itérand de maxrun

part

["LD", "MD", "LMD"] une liste de string indiquant quelle partie du spectre de A doit être déflatée

"LD" déflation depuis la lowest dominant eigenvalue, "MD" depuis la most dominant eigenvalue et "LMD" depuis "LD" et "MD"
Ipart l'itérand de part

SI

['S0', 'S1', 'S2', 'S3'] liste de string pour gérer les paramètres de déflation (expliqué précédemment)

ISI l'itérand sur SI

Ikliste

liste pour le nombre de vecteurs à déflater

Ik et k les itérands sur Ikliste

(no term)

i l'itérand sur le nombre de seconds membres

Le schéma général pour les expériences est le suivant:

for Imat in mat: # Matrix loop
    for Itheta in theta: # Matrix parameter loop
        A,DA,QTA = setmatrix4exp(which=Imat, n=n, k=nk, theta=Itheta) # Matrix constructor
        for ISI in SI: # Deflation parameters loop
            for Ipart in part: # Deflation parameter loop
                for Ik in Ikliste:
                    #Set k; l1; l2; RHR
                    for Irun in maxrun:
                        xexa = np.random.rand(A.shape[0],nrhs)
                        b    = A@xexa
                        # ser Deflation parameters
                        for i in nrhs:
                            s =  ConjGrad() # set solver
                            x, W = s.dot(b[...,i:i+1]) # solve
                            csv_eigen(A, W, DA, QTA, strat, tags=Tags) # Eigenvalues/vectors data
# R script execution using system command lines

Table of Contents