Discussion :

[Seltymar]

Bonjour,
je dois résoudre une équation linéaire utilisant le gradient conjugué. Je tiens à préciser que les math ne sont pas mon point fort.

voilà mon équation :
a = argmin a^T*L*a + (a-â)^T*D*(a-â) + (a-â)^T*F*(a-â)

Je dois donc trouver a sachant que je connais les matrices L, D et F de tailles N*N ainsi que le vecteur colonne â de taille N. L est définie symétrique positive et D et F sont toutes deux des matrices diagonales.
a est un vecteur colonne de taille N et connu partiellement.


Sur les cours que j'ai trouvé sur internet, ils parlent de l'équation A*x=b mais je n'arrive pas à faire le lien avec mon équation.

Toute aide est la bienvenue.
Merci.

PS: Pour resituer le contexte, j'essaye de déterminer la composante alpha d'une image en implémentant l'algorithme "Shared Sampling for Real-Time Alpha Matting". L'équation ci-dessus est décrite dans la partie 4.2 de la publication (j'ai supprimé les constantes pour ne pas surcharger).

[Aleph69]

Bonjour,

ta fonctionnelle J atteint son minimum en un point stationnaire, c'est-à-dire en un point a annulant le gradient de J.
En calculant ce gradient, tu vas tomber sur un système d'équations linéaires du type A*x=b.

[Seltymar]

Merci pour ta réponse.

Je suis pas sur d'avoir tout compris alors j'ai essayé de le faire sur un exemple simple. Corrige moi si je dis des âneries.
Prenons J(a)=a^T*L*a et a = (a1;a2) et
L= (l11,l12) une matrice 2*2
(l21,l22)

Donc je dois calculer le gradient de J et c'est possible car J(a) est un scalaire.
Donc je calcule d'abord J(a).

J(a) = a1^2*l11+a1*a2*l21+a1*a2*l12+a2^2*l22.

Je calcule le gradient :
grad(J)= (a1*2*l11+a2*l21+a2*l12)
(a2*2*l22+a1*l12+a1*l21)

ce qui fait
(a1*(2*l11) + a2*(l21+l12))
(a1*(l21+l12) + a2*(2*l22))

je peux transformer ça en L'*a avec L' =
(2*l11, l12+l21)
(l12+l21, 2*l22)

Comme grad(J)=0, je me retrouve avec L'*a=b avec b = 0 vecteur colonne de taille 2.

Est-ce que j'ai bon ?

Si je reprends l'équation de ma question, ça me donne :
L'*a+D'*(a-â)+ F'*(a-â)=0

Je dois maintenant appliquer l'algorithme du gradient conjugué.
Maintenant, comme mes matrices et vecteurs ont une taille au alentour de 200 000 voir plus, calculer L' peut s'avérer long non ? Y'a-t-il une méthode pour le faire ou faut le faire de façon brut ?

[Aleph69]

Bonjour,

il suffit de revenir à la définition de la différentielle d'une fonction pour faire les calculs. Ce n'est pas utile de calculer le gradient composante par composante. Pour une matrice symétrique S et un vecteur x, le gradient de la fonctionnelle J(x)=(x,S*x) appliqué en x est 2*S*x. De même, pour un vecteur b fixé, le gradient de la fonctionnelle K(x)=(b,x) appliqué en x est b. De ceci résulte, que la solution unique d'un système linéaire symétrique inversible S*x=b est également l'unique minimum de la fonctionnelle convexe L(x)=0.5*(x,S*x)-(b,x).

Tu te trouves exactement dans cette situation au coefficient 0.5 près.

Pour ta question sur la taille de ton problème, tu n'as besoin de calculer qu'un peu plus de la moitié de ta matrice par symétrie. Cela dit, si elle est dense je doute que tu puisses la stocker en mémoire vive, d'une part, et il ne faut pas utiliser la méthode des gradients conjugués mais plutôt la factorisation de Cholesky d'autre part.

[Seltymar]

Bonjour,

Effectivement, j'avais zappé que la matrice était symétrique.

Donc pour mon équation de base j'obtiens :
2*L*a+2*D*(a-â)+ 2*F*(a-â)=0
Avec la méthode des gradients conjugués, je comprends qu'il me suffit de donner un a0 initial et j'obtiens le résiduel.

Mais avec la factorisation de Cholesky je vois pas trop comment m'y prendre.
Soit L=L'L'^T avec L' la matrice triangulaire. D=D'D'^T avec D' la matrice triangulaire de D et F=F'F'^T avec F' la matrice triangulaire de F.

Si j'ai bien compris je dois résoudre 2 équations :
2*L'*x + 2*D'*y + 2*F'*z=0
L'^T*a + D'^T *(a-â) + F'^T*(a-â) = x+y+z

Mais dans la première équation j'ai 3 inconnus non ?

[Aleph69]

Bonjour,

l'équation 2*L*a+2*D*(a-â)+ 2*F*(a-â)=0 s'écrit également sous la forme A*x=b avec A=L+D+F, x=a et b=(D+F)*â. Il reste à factoriser la matrice A.

Pour information, la factorisation de Cholesky d'une matrice diagonale D est le produit de deux matrices diagonales identiques contenant les racines carrées des coefficients de D. Je le précise par rapport à ce que tu as écrit dans ton dernier message (factorisations de D et F) mais ça n'est d'aucune utilité pratique ici.

[Seltymar]

Bonjour,

Merci beaucoup pour ton aide, cela m'a bien aidé. Je me rends compte que j'ai vraiment de grosses lacunes en math particulièrement avec les matrices.

Précédemment, tu as dit qu'il fallait mieux utiliser la factorisation de Cholesky plutôt que la méthode des gradients conjugués. Pourquoi une méthode directe est-elle préférable à une méthode itérative ?

[Aleph69]

Bonjour,

ce n'est jamais évident de savoir quand utiliser un solveur direct ou un solveur itératif.

Les solveurs directs sont plus gourmands en mémoire vive. La factorisation de Cholesky est numériquement stable. Ce n'est pas le cas de la factorisation de Gauss en théorie, mais en pratique on obtient généralement de bonnes approximations des solutions des systèmes d'équations linéaires. Dans le cas de la résolution de systèmes linéaires denses, il intéressant de choisir une méthode directe car on peut stocker les facteurs dans la matrice elle-même.

Les solveurs itératifs (méthodes de Krylov) sont plus intéressants quand les matrices sont creuses avec une largeur de bande importante et quand les matrices sont singulières.

Ce ne sont que des petites règles empiriques. Tu pourras toujours trouver un système pour lequel telle ou telle méthode est meilleure (théorèmes "no free lunch").

[Seltymar]

Bonjour,

Comme ma matrice peut aller jusqu'à plus de 2 000 000 x 2 000 0000 (image HD 1080) et qu'elle est très parsemée, je vais opter pour le conjugué de gradient. De toute façon, si je veux utiliser la méthode pré-conditionnée je vais devoir passer par la décomposition de Cholesky donc je changerai peut-être de méthode à ce moment là. Mais dans tous les cas, je crains un temps de calcul important si je ne parallélise pas tous mes calculs sur GPU.

Merci pour ton aide. Je mets le sujet en résolu.