L'objectif est de minimiser la fonction
où A est une matrice carrée symétrique définie positive de taille n.
Le calcul montre qu'une solution du problème est la solution du système
: en effet, on a
.
Intuitivement, la fonction f peut donc être vue comme une primitive (littéralement un potentiel scalaire) du résidu
. En annulant le gradient de f, on obtient le vecteur x qui minimise l'erreur.
La méthode du gradient conjugué vue comme une méthode directe
On rappelle que deux vecteurs non nuls u et v sont conjugués par rapport à A si

Sachant que A est symétrique définie positive, on en déduit un produit scalaire

Ainsi, deux vecteurs sont conjugués s'ils sont orthogonaux pour ce produit scalaire.
La conjugaison est une relation symétrique : si u est conjugué à v pour A, alors v est conjugué à u.
Supposons que {pk} est une suite de n directions conjuguées deux à deux. Alors les {pk} forment une base de Rn, ainsi la solution x* de Ax= b dans cette base :

Les coefficients sont donnés par

(car
sont conjugués deux à deux)

On a ainsi l'idée directrice de la méthode pour résoudre le système Ax=b : trouver une suite de n directions conjuguées, et calculer les coefficients αk.
La méthode du gradient conjugué vue comme une méthode itérative
En choisissant correctement les directions conjuguées pk, il n'est pas nécessaire de toutes les déterminer pour obtenir une bonne approximation de la solution x*. Il est ainsi possible de considérer la méthode du gradient conjugué comme une méthode itérative. Ce choix permet ainsi de considérer la résolution de systèmes de très grande taille, où le calcul de l'ensemble des directions aurait été très long.
On considère ainsi un premier vecteur x0, qu'on pourra supposer nul (sinon, il faut considérer le système Az=b − Ax0). L'algorithme va consister, partant de x0, à se « rapprocher » de la solution x* inconnue, ce qui suppose la définition d'une métrique. Cette métrique vient du fait que la solution x* est l'unique minimiseur de la forme quadratique :

Ainsi, si f(x) diminue après une itération, alors on s'approche de x*.
Ceci suggère donc de prendre la première direction p1 comme l'opposé du gradient de f à x=x0. Le gradient vaut Ax0 – b= – b, d'après notre première hypothèse. Les vecteurs suivants de la base seront ainsi conjugués au gradient, d'où le nom « méthode du gradient conjugué ».
Soit rk le résidu à la ke itération :

Notons que rk est l'opposé du gradient de f en x=xk, ainsi, l'algorithme du gradient indique d'évoluer dans la direction rk. On rappelle que les directions pk sont conjuguées deux à deux. On veut aussi que la direction suivante soit construite à partir du résidu courant et des directions précédemment construites, ce qui est une hypothèse raisonnable en pratique.
La contrainte de conjugaison est une contrainte d'orthonormalité, aussi le problème partage des similitudes avec le procédé de Gram-Schmidt.
On a ainsi

Suivant cette direction, le point suivant est donné par

- où le pas
est déterminé de manière à minimiser
;
- le minimum de g est atteint pour
et comme A est définie positive,
,
- on a donc :

- Une analyse plus détaillée de cette algorithme (un raisonnement par récurrence) montre que
est orthogonal à
, i.e.
pour
(voir ci-après) et que
est
-orthogonal à
, i.e.
, pour
.
Algorithme
Pour amorcer la récurrence, il faut partir d’une estimation initiale x0 du vecteur x recherché ; et le nombre d'itérations N nécessaire pour que
(où ε est un nombre positif arbitrairement proche de zéro) dépend du x0 choisi. Malheureusement, les méthodes de « préconditionnement » à la fois sûres et générales (c'est-à-dire efficaces pour toutes sortes de matrices symétriques positives) pour former un x0 correct sont aussi elles-mêmes coûteuses en temps de calcul. En pratique, l'intuition physique, guidée par la nature physique du problème à résoudre, suggère parfois une initialisation efficace : ces idées ont donné lieu depuis plus de trente ans à une littérature spécialisée abondante[2].
Algorithme itératif en pseudo-code
L'algorithme ci-dessous résout Ax = b, où A est une matrice réelle, symétrique, et définie positive.
Le vecteur d'entrée x0 peut être une approximation de la solution initiale ou 0. Cette algorithme est issu de la méthode itérative exacte présentée dans le paragraphe précédent, les valeurs des coefficients
semblent différentes mais en utilisant les relations de l'algorithme ci-dessous, en particulier :
, et le fait que les résidus
et
soient orthogonaux , on peut montrer par récurrence que l'on a bien :
. Ci-dessous le coefficient
est choisi pour que les résidus
et
soient orthogonaux et non pas minimiser g comme dans le paragraphe précédent, mais en fait les deux approches reviennent à la même formule pour
et le coefficient
est choisi pour que
soit A-conjugué de
.

- Cet algorithme du gradient conjugué est celui qui est très souvent utilisé. Il est intéressant de constater que le coefficient
de cet algorithme a la même expression que celui qui est utilisé dans la méthode 'Fletcher-Reeves' du gradient conjugué non-linéaire.
- On peut également remarquer que
est déduit de
en utilisant la méthode du gradient et que prendre
, revient à appliquer la méthode du gradient et peut donc être utilisé pour réinitialiser un calcul de gradient conjugué en cours. Réinitialiser un calcul en cours peut ralentir la convergence, mais peu également en augmenter la stabilité en particulier en réduisant les erreurs dus à l'accumulation d'imprécisions numériques (round-off errors).
- Les formules
et
, qui sont exactes, impliquent que les deux formules
et
sont mathématiquement équivalentes. La première formule
est utilisé dans l'algorithme pour éviter une multiplication supplémentaire par
car le vector
est déjà calculé pour évaluer
. La deuxième formule
peut par contre être plus précise car elle réduit l'accumulation des imprécisions numérique et elle est donc parfois recommandée[3].
Exemple numérique
Considérons le système linéaire Ax = b suivant :

Deux itérations de l'algorithme du gradient conjugué vont être réalisées pas à pas à en partant du vecteur initial

- Il est rappelé, que s'il n'y avait pas d'imprécision numérique dans les calculs le présent exemple serait résolu en seulement n = 2 itérations, puisque ici la matrice est de dimension 2x2 et donc que x2 devrait, aux erreurs numériques près, retrouver la solution exacte.
Solution
La solution exacte de ce système linéaire est obtenu en inversant simplement la matrice A :
et 
Les valeurs initiales de l'algorithme itératifs sont :


Lors de la première itération la valeur x1 est une meilleure approximation de la solution exacte:

et le résidu associé r1 est égal à :

Ce résidu est encore grand et donc il est nécessaire de poursuivre la deuxième itération et donc de calculer les paramètres qui vont être utile pour calculer x2: (i) d'abord β0 puis (ii) la nouvelle direction de recherche p1 et ensuite (iii) α1 :




x2 est une bien meilleure approximation de la solution que x1 et x0 et cela est numériquement déterminé en évaluant le résidu r2.

Il est également intéressant de vérifier que :
:
et

Convergence
On peut montrer le résultat suivant sur la convergence de l'algorithme :

où
désigne le condisertionnement de la matrice et 
La méthode du gradient conjugué a donc une convergence superlinéaire, qui peut être mise à mal par un mauvais conditionnement de la matrice. Elle reste toutefois meilleure que les algorithmes à direction de plus forte pente.