Método do gradiente conjugado

Uma comparação da convergência do método de descida do gradiente com tamanho de passo ótimo (em verde) e o método do gradiente conjugado (em vermelho) para a minimização da forma quadrática com um sistema linear dado. O gradiente conjugado, assumindo aritmética exata, converge em no máximo n passos onde n é o tamanho da matriz do sistema (no exemplo, n=2).

Em matemática, o método do gradiente conjugado é um algoritmo para a solução numérica de sistemas de equações lineares, cuja matriz é simétrica e positiva definida. O método do gradiente conjugado é um método iterativo, então ele pode ser aplicado a sistemas esparsos que são grandes demais para ser tratados por métodos diretos como a decomposição de Cholesky. Tais sistemas surgem frequentemente quando se resolve numericamente equações diferenciais parciais.

O método do gradiente conjugado também pode ser utilizado para resolver problemas de otimização irrestritos como minimização de energia. É comumente atribuído aos matemáticos Magnus Hestenes e Eduard Stiefel, que o programaram no computador Z4 e o estudaram extensivamente.

O método do gradiente biconjugado fornece uma generalização a matrizes não simétricas. Diversos métodos do gradiente conjugado não lineares procuram o mínimo de problemas de otimização não lineares.

Descrição do problema

Considera-se resolver o sistema de equações lineares

$$\mathbf {A} \mathbf {x} =\mathbf {b}$$

para o vetor ${\textstyle \mathbf {x} }$, onde $\mathbf {A} \in \mathbb {R} ^{n\times n}$ é uma matriz real simétrica (i.e., ${\textstyle \mathbf {A} ^{\intercal }=\mathbf {A} }$) positiva definida (i.e., ${\textstyle \mathbf {x} ^{\intercal }\mathbf {A} \mathbf {x} >0}$ para todo ${\textstyle \mathbf {x} \neq \mathbf {0} }$) e ${\textstyle \mathbf {b} \in \mathbb {R} ^{n}}$. Denota-se a solução única desse sistema por ${\textstyle \mathbf {x} _{\ast }}$.

Derivação do método do gradiente conjugado

O método do gradiente conjugado (CG) pode ser derivado sob diversas perspectivas. A abordagem a seguir busca minimizar uma função quadrática envolvendo a matriz ${\textstyle \mathbf {A} }$ e o vetor ${\textstyle \mathbf {b} }$.

O método da descida rápida

Descida do gradiente em 2D

O método CG pode ser obtido a partir do método da descida rápida para resolver um sistema positivo definido acrescido de modificações que melhoram a taxa de convergência. Para resolver o sistema linear ${\textstyle \mathbf {A} \mathbf {x} =\mathbf {b} }$, onde ${\textstyle \mathbf {A} }$ é matriz positiva definida, encontra-se o mínimo de uma função quadrática cujo mínimo coincide com a solução do sistema. Considere a função

$$\phi (\mathbf {x} )={\frac {1}{2}}\mathbf {x} ^{\intercal }\mathbf {A} \mathbf {x} -\mathbf {x} ^{\intercal }\mathbf {b} ,$$

onde ${\textstyle \mathbf {A} \in \mathbb {R} ^{n\times n}}$ e ${\textstyle \mathbf {b} \in \mathbb {R} ^{n}}$. Para minimizar ${\textstyle \phi (\mathbf {x} )}$ encontra-se o vetor ${\textstyle \mathbf {x} }$ tal que todas as derivadas parciais ${\textstyle {\frac {\partial \phi }{\partial x_{i}}}}$ sejam nulas. O gradiente de ${\textstyle \phi }$ é o vetor ${\textstyle \mathbf {\nabla } \phi =\left[{\frac {\partial \phi }{\partial x_{1}}}\;{\frac {\partial \phi }{\partial x_{2}}}\;\cdots \;{\frac {\partial \phi }{\partial x_{n-1}}}\;{\frac {\partial \phi }{\partial x_{n}}}\right]^{\intercal }}$, então ${\textstyle \mathbf {x} }$ é ponto de mínimo somente se ${\textstyle \mathbf {\nabla } \phi (\mathbf {x} )=0}$. Algumas manipulações algébricas mostram que

$$\mathbf {\nabla } \phi (\mathbf {x} )={\frac {1}{2}}\mathbf {A} ^{\intercal }\mathbf {x} +{\frac {1}{2}}\mathbf {A} \mathbf {x} -\mathbf {b} .$$

 

 

 

 

(1)

Como ${\textstyle \mathbf {A} }$ é simétrica,

$$\mathbf {\nabla } \phi (\mathbf {x} )=\mathbf {A} \mathbf {x} -\mathbf {b} ,$$

 

 

 

 

(2)

e o mínimo de ${\textstyle \phi (\mathbf {x} )}$ ocorre quando ${\textstyle \mathbf {A} \mathbf {x} =\mathbf {b} }$. Logo, o problema de otimização de encontrar o valor mínimo de uma função ${\textstyle \phi (\mathbf {x} )}$ é equivalente a resolver ${\textstyle \mathbf {A} \mathbf {x} =\mathbf {b} }$. Como ${\textstyle \mathbf {A} }$ é positiva definida, ${\textstyle \mathbf {A} }$ é não singular, e ${\textstyle \mathbf {x} _{\ast }=\mathbf {A} ^{-1}\mathbf {b} }$. Portanto,

$$\phi (\mathbf {x} _{\ast })={\frac {1}{2}}\mathbf {x} _{\ast }^{\intercal }\mathbf {A} \mathbf {A} ^{-1}\mathbf {b} -\mathbf {x} _{\ast }^{\intercal }\mathbf {b} =-{\frac {1}{2}}\mathbf {x} _{\ast }^{\intercal }\mathbf {b} .$$

 

 

 

 

(3)

Gráfico de um paraboloide elíptico

Gráfico do campo vetorial do gradiente e das linhas de contorno de um paraboloide elíptico

A Equação 1 revela o tipo de superfície gerado pela função ${\textstyle \phi (x)}$. Como ${\textstyle \mathbf {A} }$ é uma matriz real e simétrica, o Teorema Espectral garante que existe uma matriz ortogonal ${\textstyle \mathbf {P} }$ tal que ${\textstyle \mathbf {A} =\mathbf {P} \mathbf {D} \mathbf {P} ^{\intercal }}$, onde ${\textstyle \mathbf {D} }$ é uma matriz diagonal contendo os autovalores de ${\textstyle \mathbf {A} }$. Fazendo a mudança de coordenadas ${\textstyle \mathbf {u} =\mathbf {P} ^{\intercal }(\mathbf {x} -\mathbf {x} _{\ast })}$, i.e., ${\textstyle \mathbf {x} =\mathbf {P} \mathbf {u} +\mathbf {x} _{\ast }}$, e substituindo na Equação 1, obtém-se

$${\begin{aligned}{\tilde {\phi }}(u)&={\frac {1}{2}}(\mathbf {P} \mathbf {u} +\mathbf {x} _{\ast })^{\intercal }\mathbf {A} (\mathbf {P} \mathbf {u} +\mathbf {x} _{\ast })-(\mathbf {P} \mathbf {u} +\mathbf {x} _{\ast })^{\intercal }\mathbf {b} \\&={\frac {1}{2}}\mathbf {u} ^{\intercal }\mathbf {D} \mathbf {u} -{\frac {1}{2}}\mathbf {x} _{\ast }^{\intercal }\mathbf {b} \\&={\frac {1}{2}}\sum _{i=1}^{n}\lambda _{i}\mathbf {u} _{i}^{2}-{\frac {1}{2}}\mathbf {x} _{\ast }^{\intercal }\mathbf {b} ,\end{aligned}}$$

notando que ${\textstyle \mathbf {D} =\mathbf {P} ^{\intercal }\mathbf {A} \mathbf {P} }$ e que ${\textstyle \mathbf {A} \mathbf {x} _{\ast }=\mathbf {b} }$. Tem-se que ${\textstyle \lambda _{i}>0,1\leq i\leq n}$, porque ${\textstyle \mathbf {A} }$ é positiva definida. Para ${\textstyle n=2}$, ${\textstyle {\tilde {\phi }}(u)=\lambda _{1}\mathbf {u} _{1}^{2}+\lambda _{2}\mathbf {u} _{2}^{2}-{\frac {1}{2}}\mathbf {x} _{\ast }^{\intercal }\mathbf {b} }$ representa um paraboloide elíptico. Em geral, a transformação ortogonal ${\textstyle \mathbf {u} =\mathbf {P} ^{\intercal }(\mathbf {x} -\mathbf {x} _{\ast })}$ preserva as distâncias e os ângulos entre os vetores, então ${\textstyle \phi (x)}$ é um paraboloide no ${\textstyle \mathbb {R} ^{n}}$. Isso garante que o ponto crítico onde o gradiente se anula (${\textstyle \mathbf {\nabla } \phi (x)=\mathbf {0} }$) é o mínimo global da função.

Na descida rápida, o valor mínimo ${\textstyle \phi (\mathbf {x} _{\ast })}$ é encontrado iterativamente, tomando um ponto inicial e decaindo em direção a ${\textstyle \mathbf {x} _{\ast }}$. Como a função quadrática ${\textstyle \phi }$ decresce mais rapidamente na direção negativa do gradiente, a Equação 2 mostra que ${\textstyle -\mathbf {\nabla } \phi (\mathbf {x} _{i})=\mathbf {b} -\mathbf {A} \mathbf {x} _{i}}$. Seja

$$\mathbf {r} _{i}=\mathbf {b} -\mathbf {A} \mathbf {x} _{i}$$

o resíduo no ponto ${\textstyle \mathbf {x} _{i}}$. Então o próprio resíduo aponta na direção de maior decrescimento da função partindo de ${\textstyle \mathbf {x} _{i}}$ e é chamado de vetor diretor. Se o resíduo é zero, o problema está resolvido, caso contrário, move-se ao longo da reta determinada pelo vetor diretor. Seja ${\textstyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+\alpha _{i}\mathbf {r} _{i}}$, onde ${\textstyle \alpha _{i}}$ é um escalar que determina o tamanho do passo na direção do resíduo ${\textstyle \mathbf {r} _{i}}$. No método da descida rápida, ${\textstyle \alpha _{i}}$ minimiza ${\textstyle \phi }$ para cada passo na reta correspondente. Para isso, considere a função

$$f(\alpha _{i})=\phi (\mathbf {x} _{i}+\alpha _{i}\mathbf {r} _{i}).$$

Usando a simetria de ${\textstyle \mathbf {A} }$ e a comutatividade do produto interno, tem-se

$${\begin{aligned}f(\alpha _{i})&={\frac {1}{2}}(\mathbf {x} _{i}+\alpha _{i}\mathbf {r} _{i})^{\intercal }\mathbf {A} (\mathbf {x} _{i}+\alpha _{i}\mathbf {r} _{i})-(\mathbf {x} _{i}+\alpha _{i}\mathbf {r} _{i})^{\intercal }\mathbf {b} \\&={\frac {1}{2}}\mathbf {x} _{i}^{\intercal }\mathbf {A} \mathbf {x} _{i}+{\frac {\alpha _{i}^{2}}{2}}\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}+\alpha _{i}\mathbf {x} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}-\mathbf {x} _{i}^{\intercal }\mathbf {b} -\alpha _{i}\mathbf {r} _{i}^{\intercal }\mathbf {b} \end{aligned}}$$

O mínimo de ${\textstyle \phi }$ ocorre quando ${\textstyle f'(\alpha _{i})=0}$, e

$${\begin{aligned}f'(\alpha _{i})&=\alpha _{i}\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}+\mathbf {x} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}-\mathbf {r} _{i}^{\intercal }\mathbf {b} \\&=\alpha _{i}\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}+\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {x} _{i}-\mathbf {r} _{i}^{\intercal }\mathbf {b} \\&=\alpha _{i}\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}+\mathbf {r} _{i}^{\intercal }(\mathbf {A} \mathbf {x} _{i}-\mathbf {b} )\\&=\alpha _{i}\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}-\mathbf {r} _{i}^{\intercal }\mathbf {r} _{i}\end{aligned}}$$

Então, o mínimo ocorre quando ${\textstyle \alpha _{i}={\frac {\mathbf {r} _{i}^{\intercal }\mathbf {r} _{i}}{\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}}}}$. Move-se de ${\textstyle \mathbf {x} _{i}}$ a ${\textstyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+\alpha _{i}\mathbf {r} _{i}}$ e depois usa ${\textstyle \mathbf {x} _{i+1}}$ e ${\textstyle \mathbf {r} _{i+1}}$ para computar ${\textstyle \mathbf {x} _{i+2}=\mathbf {x} _{i+1}+\alpha _{i+1}\mathbf {r} _{i+1}}$, o próximo ponto em direção ao mínimo. Esse desenvolvimento motiva o algoritmo a seguir.

Entradas

matriz ${\textstyle \mathbf {A} \in \mathbb {R} ^{n\times n}}$ simétrica positiva definida, lado direito ${\textstyle \mathbf {b} \in \mathbb {R} ^{n}}$, aproximação inicial ${\textstyle \mathbf {x} _{0}}$, tolerância tol e número máximo de iterações maxiter.

Saídas

solução aproximada ${\textstyle \mathbf {x} }$, norma-2 do resíduo relativo ${\textstyle r}$ e número de iterações ${\textstyle i}$

Inicialização

$\mathbf {r} _{0}:=\mathbf {b} -\mathbf {A} \mathbf {x} _{0}$

$i:=0$

Algoritmo

Enquanto $\mathbf {r} _{i}$ for maior que a tolerância desejada e o número de iterações for menor que o número máximo de iterações, faça:

$\alpha _{i}:={\frac {\mathbf {r} _{i}^{\intercal }\mathbf {r} _{i}}{\mathbf {r} _{i}^{\intercal }A\mathbf {r} _{i}}}$

$\mathbf {x} _{i+1}:=\mathbf {x} _{i}+\alpha _{i}\mathbf {r} _{i}$

$\mathbf {r} _{i+1}:=\mathbf {b} -\mathbf {A} \mathbf {x} _{i+1}$

${\textstyle i=i+1}$

Se $\mathbf {r} _{i+1}$ for suficientemente pequeno, então pare o ciclo.

Fim Enquanto

Se o número de iterações for maior que o número máximo de iterações, faça ${\textstyle i=-1}$

Retorna ${\textstyle \mathbf {x} =\mathbf {x} _{i+1}}$, ${\textstyle r={\frac {\Vert \mathbf {b} -\mathbf {A} \mathbf {x} \Vert _{2}}{\Vert \mathbf {b} \Vert _{2}}}}$, ${\textstyle i}$

O método da descida rápida converge globalmente. Se $\lambda _{\max }(\mathbf {A} )$ e $\lambda _{\min }(\mathbf {A} )$ são, respectivamente, o maior e o menor autovalor de uma matriz simétrica $\mathbf {A}$, o número de condicionamento da matriz $\mathbf {A}$ é dado por

$\kappa (\mathbf {A} )={\frac {\lambda _{\max }(\mathbf {A} )}{\lambda _{\min }(\mathbf {A} )}}.$

Pode ser mostrado que

$$\left(\phi (\mathbf {x} _{i})+{\frac {1}{2}}\mathbf {b} ^{\intercal }\mathbf {A} ^{-1}\mathbf {b} \right)\leq \left(1-{\frac {1}{\kappa (\mathbf {A} )}}\right)\left(\phi (\mathbf {x} _{i-1})+{\frac {1}{2}}\mathbf {b} ^{\intercal }\mathbf {A} ^{-1}b\right).$$

 

 

 

 

(3)

Definindo a função erro ${\textstyle E(\mathbf {x} _{i})=\phi (\mathbf {x} _{i})-\phi (\mathbf {x} _{\ast })=\phi (\mathbf {x} _{i})+{\frac {1}{2}}\mathbf {b} ^{\intercal }\mathbf {A} ^{-1}\mathbf {b} }$, a Equação 3 pode ser reescrita como

$$E(\mathbf {x} _{i})\leq \left(1-{\frac {1}{\kappa (\mathbf {A} )}}\right)E(\mathbf {x} _{i-1}),$$

o que indica que o algoritmo da descida rápida converge independentemente do chute inicial ${\textstyle \mathbf {x} _{0}}$, pois ${\textstyle \kappa (\mathbf {A} )\geq 1}$. Se ${\textstyle \kappa (\mathbf {A} )}$ é grande, então ${\textstyle 1-1/\kappa (\mathbf {A} )}$ é próximo de 1 e a convergência é lenta. Isso pode ser visto geometricamente. Para qualquer matriz ${\textstyle n\times n}$, as direções de buscas sucessivas ${\textstyle \mathbf {r} _{i}}$ são ortogonais. Sabendo que ${\textstyle \mathbf {A} \mathbf {x} _{i}=\mathbf {b} -\mathbf {r} _{i}}$,

$${\begin{aligned}\mathbf {r} _{i}^{\intercal }\mathbf {r} _{i+1}&=\mathbf {r} _{i}^{\intercal }(\mathbf {b} -\mathbf {A} \mathbf {x} _{i+1})\\&=\mathbf {r} _{i}^{\intercal }[\mathbf {b} -\mathbf {A} (\mathbf {x} _{i}+\alpha _{i}\mathbf {r} _{i})]\\&=\mathbf {r} _{i}^{\intercal }\mathbf {b} -\mathbf {r} _{i}^{\intercal }(\mathbf {A} \mathbf {x} _{i}+\alpha _{i}\mathbf {A} \mathbf {r} _{i})\\&=\mathbf {r} _{i}^{\intercal }\mathbf {b} -\mathbf {r} _{i}^{\intercal }[(\mathbf {b} -\mathbf {r} _{i})+\alpha _{i}\mathbf {A} \mathbf {r} _{i}]\\&=\mathbf {r} _{i}^{\intercal }\mathbf {b} -\mathbf {r} _{i}^{\intercal }\mathbf {b} +\mathbf {r} _{i}^{\intercal }\mathbf {r} _{i}-\alpha _{i}\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}\\&=\mathbf {r} _{i}^{\intercal }\mathbf {r} _{i}-\alpha _{i}\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}\end{aligned}}$$

Substituindo pelo valor de ${\textstyle \alpha _{i}}$, vem

$$\mathbf {r} _{i}^{\intercal }\mathbf {r} _{i}-\alpha _{i}\mathbf {r} _{i}^{\intercal }\mathbf {A} \mathbf {r} _{i}=0,$$

então ${\textstyle \mathbf {r} _{i}^{\intercal }\mathbf {r} _{i+1}=0}$. Se ${\textstyle \mathbf {x} _{i}}$ é o ponto atual na descida ao mínimo, o vetor de direção de busca ${\textstyle \mathbf {r} _{i}}$ é ortogonal ao conjunto de nível ${\textstyle \phi (\mathbf {x} _{i})=k_{i}}$. A próxima direção de busca ${\textstyle \mathbf {r} _{i+1}}$ é ortogonal à direção de busca anterior ${\textstyle \mathbf {r} _{i}}$ e ao conjunto de nível ${\textstyle \phi (\mathbf {x} _{i+1})=k_{i+1}}$. Para ${\textstyle n=2}$, esse caminho segue um padrão de zigue-zague que pode não ser o mais eficiente até o mínimo. Como os autovalores de uma matriz simétrica ${\textstyle \mathbf {A} }$ são iguais aos seus valores singulares, a forma das elipses que configuram as linhas de contorno da função ${\textstyle \phi (\mathbf {x} )}$ dependem do número de condicionamento da matriz. Se ${\textstyle \kappa (\mathbf {A} )}$ é pequeno, as elipses são próximas de círculos, e a descida rápida converge rapidamente ao mínimo do paraboloide. Se ${\textstyle \kappa (\mathbf {A} )}$ é grande, essas elipses se tornam mais alongadas e estreitas; o paraboloide tem um cânion perto do mínimo que faz os valores de ${\textstyle \mathbf {x} _{i}}$ se moverem de um lado para o outro nas paredes íngremes do cânion, desacelerando a convergência.

Exemplo numérico

Comparação das linhas de contorno das formas quadráticas de um sistema bem condicionado A (à esquerda) e um mal condicionado B (à direita).

Aplicando o método da descida rápida às matrizes ${\textstyle \mathbf {A} ={\begin{bmatrix}5&1\\1&4\end{bmatrix}}}$, ${\textstyle \kappa (\mathbf {A} )=1.6612}$ e ${\textstyle \mathbf {B} ={\begin{bmatrix}7&6.99\\6.99&7\end{bmatrix}}}$, ${\textstyle \kappa (\mathbf {B} )=1399}$, com ${\textstyle \mathbf {b} ={\begin{bmatrix}46.463\\17.499\end{bmatrix}}}$, ${\textstyle \mathbf {x} _{0}={\begin{bmatrix}5\\5\end{bmatrix}}}$ e tolerância de ${\textstyle 1.0\times 10^{-12}}$ em ambos os casos, obtém-se a aproximação desejada em 21 iterações para a matriz ${\textstyle \mathbf {A} }$ e em 8462 iterações para a matriz ${\textstyle \mathbf {B} }$.

Uma analogia para entender a descida rápida

Neblina nas montanhas

A intuição por trás da descida rápida pode ser ilustrada pelo seguinte cenário hipotético. Um grupo de pessoas encontra-se perdido nas montanhas e tenta descer de volta à entrada da trilha (i.e., tenta achar o mínimo global). Há uma densa neblina que torna a visibilidade extremamente baixa, então o caminho até a base da montanha não é visível e elas precisam usar informação local para se guiarem. Elas podem usar o método da descida rápida, que envolve calcular a declividade da ladeira na sua posição atual e prosseguir na direção com a maior declividade. Caso elas quisessem ir ao topo da montanha (i.e., o máximo global), bastava seguir a direção com maior aclive. Usando esse método, elas eventualmente encontrariam o caminho da descida, ou então ficariam presas em algum buraco (i.e., ponto de mínimo local ou ponto de sela), como um lago da montanha. No entanto, suponha que encontrar a direção de maior declive não é imediatamente óbvio por mera observação, mas que requer algum tipo de instrumento sofisticado para calcular que uma das pessoas coincidiu de trazer consigo para a trilha. Calcular a declividade em cada posição toma certo tempo, então elas devem minimizar o uso do instrumento caso elas queiram descer da montanha antes do pôr do sol. A dificuldade é então escolher a frequência com a qual calcular a declividade na posição atual, de modo a não perderem o rumo.

Nessa analogia, as pessoas representam o algoritmo, e o caminho a ser tomado para descer a montanha representa a sequência de configurações que o algoritmo vai explorar. A declividade da ladeira representa a inclinação da função naquele ponto. O instrumento usado para medir a declividade é a diferenciação. A direção que o grupo escolhe para descer a montanha é a direção do gradiente negativo da função naquele ponto. O tempo que elas caminham antes de fazer outra medição é o tamanho do passo.

Da descida rápida ao gradiente conjugado

O método da descida rápida faz uma busca baseada no gradiente conforme ${\textstyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+\alpha _{i}\mathbf {r} _{i}}$, onde os ${\textstyle \mathbf {r} _{i}}$ são os vetores do resíduo. Esse algoritmo é um algoritmo guloso; ele resolve o problema de minimização de ${\textstyle \phi }$ apenas localmente. Uma vez que as direções de busca ${\textstyle \mathbf {r} _{i}}$ são sucessivamente ortogonais, o caminho em zigue-zague pode desfazer o progresso passado e ainda repetir direções na descida ao mínimo da função ${\textstyle \phi }$. Uma alternativa é então escolher direções de busca de forma a não perder a minimização atingida em qualquer direção de busca anterior.

O método do gradiente conjugado computa ${\textstyle \mathbf {x} _{i+1}=\mathbf {x} _{i}+\alpha _{i}\mathbf {p} _{i}}$, onde os vetores diretores ${\textstyle \mathbf {p} _{i}}$ são escolhidos tal que as aproximações ${\textstyle \mathbf {x} _{i}}$ decrescem muito mais rapidamente ao mínimo de ${\textstyle \phi }$. Isso pode ser feito exigindo que as direções de busca ${\textstyle \mathbf {p} _{i}}$ sejam conjugadas em relação a uma matriz ${\textstyle \mathbf {A} }$, i.e., elas devem ser ortogonais após serem esticadas e distorcidas por ${\textstyle \mathbf {A} }$. Desse modo, o método utiliza informações dos passos anteriores e evita o padrão de zigue-zague característico da descida rápida.

O gradiente conjugado como um método direto

Embora seja usado como um método iterativo, o método do gradiente conjugado foi originalmente desenvolvido por Hestenes e Stiefel como um método direto para resolver sistemas lineares ${\textstyle n\times n}$ positivo definidos que encontra a solução exata em ${\textstyle n}$ passos, assim como a eliminação gaussiana com pivotamento parcial. Seja qual for a abordagem, toda derivação do método do gradiente conjugado se apoia em duas propriedades fundamentais:

1. Os vetores do resíduo ${\textstyle \{\mathbf {r} _{i}\}}$ são ortogonais: ${\textstyle \mathbf {r} _{i}^{\intercal }\mathbf {r} _{j}=0}$ (${\textstyle i>j}$);
2. Os vetores diretores ${\textstyle \{\mathbf {p} _{i}\}}$ são ${\textstyle \mathbf {A} }$-conjugados: ${\textstyle \mathbf {p} _{i}^{\intercal }\mathbf {p} _{j}=0}$ (${\textstyle i>j}$).

Para saber como escolher as direções de busca ${\textstyle \mathbf {p} _{i}}$, é necessário definir a norma-${\textstyle \mathbf {A} }$ (ou norma de energia). Dizemos que dois vetores não nulos ${\textstyle \mathbf {u} }$ e ${\textstyle \mathbf {v} }$ são conjugados em relação a uma matriz ${\textstyle \mathbf {A} }$ (ou ${\textstyle \mathbf {A} }$-conjugados) se

$$\mathbf {u} ^{\intercal }\mathbf {A} \mathbf {v} =0.$$

Sendo ${\textstyle \mathbf {A} }$ simétrica e positiva definida, o lado esquerdo da equação define um produto interno

$$\mathbf {u} ^{\intercal }\mathbf {A} \mathbf {v} =\langle \mathbf {u} ,\mathbf {v} \rangle _{\mathbf {A} }:=\langle \mathbf {A} \mathbf {u} ,\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} ^{\intercal }\mathbf {v} \rangle =\langle \mathbf {u} ,\mathbf {A} \mathbf {v} \rangle .$$

Então, dois vetores são ${\textstyle \mathbf {A} }$-conjugados se, e somente se, eles são ortogonais com respeito a esse produto interno. Essa é uma relação simétrica: se ${\textstyle \mathbf {u} }$ é ${\textstyle \mathbf {A} }$-conjugado a ${\textstyle \mathbf {v} }$, então ${\textstyle \mathbf {v} }$ é ${\textstyle \mathbf {A} }$-conjugado a ${\textstyle \mathbf {u} }$. Suponha que

$${\mathcal {P}}=\{\mathbf {p} _{1},...,\mathbf {p} _{n}\}$$

é um conjunto de ${\textstyle n}$ vetores ${\textstyle \mathbf {A} }$-conjugados dois a dois, i.e., ${\textstyle \langle \mathbf {p} _{i},\mathbf {p} _{j}\rangle _{\mathbf {A} }=0}$ para todo ${\textstyle i\neq j}$. Então ${\textstyle {\mathcal {P}}}$ forma uma base para ${\textstyle \mathbb {R} ^{n}}$, e podemos escrever a solução ${\textstyle \mathbf {x} _{\ast }}$ de ${\textstyle \mathbf {A} \mathbf {x} =\mathbf {b} }$ nessa base:

$$\mathbf {x} _{\ast }=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{i}\Rightarrow \mathbf {A} \mathbf {x} _{\ast }=\sum _{i=1}^{n}\alpha _{i}\mathbf {A} \mathbf {p} _{i}$$

Multiplicando pela esquerda a equação ${\textstyle \mathbf {A} \mathbf {x} =\mathbf {b} }$ por ${\textstyle \mathbf {p} _{j}^{\intercal }}$, fica

$$\mathbf {p} _{j}^{\intercal }\mathbf {b} =\mathbf {p} _{j}^{\intercal }\mathbf {A} \mathbf {x} _{\ast }=\sum _{i=1}^{n}\alpha _{i}\mathbf {p} _{j}^{\intercal }\mathbf {A} \mathbf {p} _{i}=\sum _{i=1}^{n}\alpha _{i}\langle \mathbf {p} _{j},\mathbf {p} _{i}\rangle _{\mathbf {A} }=\alpha _{j}\langle \mathbf {p} _{j},\mathbf {p} _{j}\rangle _{\mathbf {A} }$$

Isolando ${\textstyle \alpha _{j}}$, chegamos em

$$\alpha _{j}={\frac {\mathbf {p} _{j}^{\intercal }\mathbf {b} }{\langle \mathbf {p} _{j},\mathbf {p} _{j}\rangle _{\mathbf {A} }}}={\frac {\langle \mathbf {p} _{j},\mathbf {b} \rangle }{\langle \mathbf {p} _{j},\mathbf {p} _{j}\rangle _{\mathbf {A} }}}.$$

Isso resulta no seguinte método para resolver a equação matricial ${\textstyle \mathbf {A} \mathbf {x} =\mathbf {b} }$: encontre uma sequência de ${\textstyle n}$ direções conjugadas, e então calcule os coeficientes ${\textstyle \alpha _{j}}$.

O gradiente conjugado como um método iterativo

Se os vetores conjugados ${\textstyle \mathbf {p} _{i}}$ forem escolhidos cuidadosamente, é possível que nem todos eles sejam necessários para obter uma boa aproximação para a solução ${\textstyle \mathbf {x} _{\ast }}$. Assim, o método do gradiente conjugado pode ser considerado como um método iterativo. Isso possibilita encontrar soluções aproximadas de sistemas de ordem ${\textstyle n}$ suficientemente grande tal que o método direto levaria muito tempo para resolver.

Simlarmente ao método da descida rápida, o algoritmo começa com uma aproximação ${\textstyle \mathbf {x} _{0}}$. A primeira direção de busca é tomada como ${\textstyle \mathbf {p} _{0}=\mathbf {b} -\mathbf {A} \mathbf {x} _{0}}$. Os outros vetores da base serão conjugadas ao gradiente, então o nome método do gradiente conjugado. Observa que ${\textstyle \mathbf {p} _{0}}$ é também o resíduo proveniente do passo inicial do algoritmo.

Seja ${\textstyle \mathbf {r} _{i}=\mathbf {b} -\mathbf {A} \mathbf {x} _{i}}$ o resíduo no i-ésimo passo. Como notado acima, ${\textstyle \mathbf {r} _{i}}$ é o negativo do gradiente de ${\textstyle \phi }$ no ponto ${\textstyle \mathbf {x} _{i}}$, então a descida rápida exigiria um passo na direção ${\textstyle \mathbf {r} _{i}}$. Dessa vez, no entanto, faz-se com que as direções ${\textstyle \mathbf {p} _{i}}$ sejam ${\textstyle \mathbf {A} }$-conjugadas. Uma forma de fazer isso é exigindo que a próxima direção de busca seja construída do resíduo e de todas as direções de busca anteriores. Essa restrição de conjugação é um tipo de restrição de ortogonalização, logo o algoritmo pode ser visto como um exemplo do processo de Gram-Schmidt. Isso leva à expressão:

$$\mathbf {p} _{i}=\mathbf {r} _{i}-\sum _{j<i}c_{ij}\mathbf {p} _{j},\quad c_{ji}={\frac {\langle \mathbf {r} _{i},\mathbf {p} _{j}\rangle _{\mathbf {A} }}{\langle \mathbf {p} _{j},\mathbf {p} _{j}\rangle _{\mathbf {A} }}}.$$

 

 

 

 

(4)

Seguindo essa direção, o próximo ponto é dado por

$$\mathbf {x} _{i+1}=\mathbf {x} _{i}+\alpha _{i}\mathbf {p} _{i},$$

com

$$\alpha _{i}={\frac {\mathbf {p} _{i}^{\intercal }(\mathbf {b} -\mathbf {Ax} _{i})}{\mathbf {p} _{i}^{\intercal }\mathbf {A} \mathbf {p} _{i}}}={\frac {\mathbf {p} _{i}^{\intercal }\mathbf {r} _{i}}{\mathbf {p} _{i}^{\intercal }\mathbf {A} \mathbf {p} _{i}}}.$$

A expressão para ${\textstyle \alpha _{i}}$ pode ser deduzida da mesma forma que na descida rápida, substituindo ${\textstyle \mathbf {r} _{i}}$ por ${\textstyle \mathbf {p} _{i}}$ e minimizando a função ${\textstyle f(\alpha _{i})=\phi (\mathbf {x} _{i+1})=\phi (\mathbf {x} _{i}+\alpha _{i}\mathbf {p} _{i})}$. Com algumas manipulações algébricas, é possível simplificar essa expressão para

$$\alpha _{i}={\frac {\Vert \mathbf {r} _{i}\Vert _{2}^{2}}{\langle \mathbf {A} \mathbf {p} _{i},\mathbf {p} _{i}\rangle }}.$$

Da Equação 4, pode parecer que o algoritmo descrito requer o armazenamento de todos os vetores diretores e resíduos, como também de várias multiplicações matriz-vetor, o que apresenta um custo computacional alto. Entretanto, da ortogonalidade dos resíduos pode-se mostrar que

$$c_{ji}={\frac {\langle \mathbf {r} _{i},\mathbf {p} _{j}\rangle _{\mathbf {A} }}{\langle \mathbf {p} _{j},\mathbf {p} _{j}\rangle _{\mathbf {A} }}}={\begin{cases}-{\frac {\Vert \mathbf {r} _{i}\Vert _{2}^{2}}{\Vert \mathbf {r} _{i-1}\Vert _{2}^{2}}},\quad j=i-1\\0,\quad j<i-1,\end{cases}}$$

ou seja, é suficiente computar apenas o coeficiente ${\textstyle c_{i-1,i}}$ da última direção de busca ${\textstyle \mathbf {p} _{i-1}}$ e ainda se ter a propriedade de conjugação das direções de busca. Definindo

$$\beta _{i}={\frac {\Vert \mathbf {r} _{i}\Vert _{2}^{2}}{\Vert \mathbf {r} _{i-1}\Vert _{2}^{2}}},$$

obtém-se a fórmula para atualizar a direção de busca

$$\mathbf {p} _{i}=\mathbf {r} _{i}+\beta _{i}\mathbf {p} _{i-1},$$

que leva ao seguinte algoritmo.

Entradas

matriz ${\textstyle \mathbf {A} \in \mathbb {R} ^{n\times n}}$ simétrica positiva definida, lado direito ${\textstyle \mathbf {b} \in \mathbb {R} ^{n}}$, aproximação inicial ${\textstyle \mathbf {x} _{0}}$, tolerância tol e número máximo de iterações maxiter

Saídas

solução aproximada ${\textstyle \mathbf {x} }$, norma-2 do resíduo relativo ${\textstyle r}$ e número de iterações ${\textstyle i}$

Inicialização

$$\mathbf {r} _{0}:=\mathbf {b} -\mathbf {A} \mathbf {x} _{0}$$ $$\mathbf {p} _{0}:=\mathbf {r} _{0}$$

$i:=0$

Algoritmo

Enquanto ${\textstyle \mathbf {r} _{i}}$ for maior que a tolerância desejada e o número de iterações for menor que o número máximo de iterações, faça:

$${\begin{aligned}&\alpha _{i}:={\frac {\mathbf {r} _{i}^{\intercal }\mathbf {r} _{i}}{\mathbf {p} _{i}^{\intercal }\mathbf {A} \mathbf {p} _{i}}}\\&\mathbf {x} _{i+1}:=\mathbf {x} _{i}+\alpha _{i}\mathbf {p} _{i}\\&\mathbf {r} _{i+1}:=\mathbf {r} _{i}-\alpha _{i}A\mathbf {p} _{i}\\&\beta _{i}:={\frac {\mathbf {r} _{i+1}^{\intercal }\mathbf {r} _{i+1}}{\mathbf {r} _{i}^{\intercal }\mathbf {r} _{i}}}\\&\mathbf {p} _{i+1}:=\mathbf {r} _{i+1}+\beta _{i}\mathbf {p} _{i}\\&i:=i+1\end{aligned}}$$

Se ${\textstyle \mathbf {r} _{i+1}}$ for suficientemente pequeno, então pare o ciclo.

Fim Enquanto

Se o número de iterações for maior que o número máximo de iterações, faça ${\textstyle i=-1}$

Retorna ${\textstyle \mathbf {x} =\mathbf {x} _{i+1}}$, ${\textstyle r={\frac {\Vert \mathbf {b} -\mathbf {A} \mathbf {x} \Vert _{2}}{\Vert \mathbf {b} \Vert _{2}}}}$, ${\textstyle i}$

Reinicializações

Por conta de erros de arredondamento inerentes à aritmética de ponto flutuante, os vetores ${\textstyle \mathbf {p} _{i}}$ podem perder a propriedade de serem ${\textstyle \mathbf {A} }$-conjugados. Uma reinicialização consiste em zerar o coeficiente ${\textstyle \beta _{i}}$ após ${\textstyle n}$ iterações (ou quando o resíduo estagnar), isto é, descartar todas as direções de busca anteriores ${\textstyle \mathbf {p} _{i-1},\mathbf {p} _{i-2},\ldots ,\mathbf {p} _{0}}$ e recomeçar o algoritmo. Isso faz com que a próxima direção de busca ${\textstyle \mathbf {p} _{i+1}=\mathbf {r} _{i+1}}$ tenha a direção do resíduo, que é exatamente um passo da descida rápida. O uso de reinicializações é uma estratégia simples e barata para limpar os erros de arredondamento e melhorar a estabilidade do algoritmo, embora possa desacelerar a convergência.

Exemplo numérico

Considere o sistema linear ${\textstyle \mathbf {A} \mathbf {x} =\mathbf {b} }$ dado por

$\mathbf {A} \mathbf {x} ={\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\end{bmatrix}}={\begin{bmatrix}1\\2\end{bmatrix}}=\mathbf {b} .$

Faremos duas iterações do método do gradiente conjugado, começando com a aproximação inicial ${\textstyle \mathbf {x} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}}$ para encontrar uma solução aproximada do sistema.

Solução

Para fins comparativos, a solução exata é

$\mathbf {x} ={\begin{bmatrix}{\frac {1}{11}}\\{\frac {7}{11}}\end{bmatrix}}\approx {\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.$

O primeiro passo é computar o vetor resíduo ${\textstyle \mathbf {\mathbf {r} _{0}} }$ associado a ${\textstyle \mathbf {\mathbf {x} _{0}} }$. Usando a fórmula ${\textstyle \mathbf {\mathbf {r} _{0}} =\mathbf {b} -\mathbf {A} \mathbf {\mathbf {x} _{0}} }$, obtém-se

$\mathbf {r} _{0}={\begin{bmatrix}1\\2\end{bmatrix}}-{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}2\\1\end{bmatrix}}={\begin{bmatrix}-8\\-3\end{bmatrix}}=\mathbf {p} _{0}.$

Como essa é a primeira iteração, usamos o resíduo ${\textstyle \mathbf {\mathbf {r} _{0}} }$ como a direção de busca inicial ${\textstyle \mathbf {\mathbf {p} _{0}} }$; o método para calcular ${\textstyle \mathbf {\mathbf {p} _{i}} }$ será outro nas iterações seguintes.

Agora computamos o escalar ${\textstyle \alpha _{0}}$ (tamanho do passo na direção de ${\textstyle \mathbf {\mathbf {p} _{0}} }$) através da relação

$\alpha _{0}={\frac {\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}{\mathbf {p} _{0}^{\mathsf {T}}\mathbf {Ap} _{0}}}={\frac {{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}={\frac {73}{331}}\approx 0.2205.$

Calculamos agora ${\textstyle \mathbf {x} _{1}}$ com a fórmula

$\mathbf {x} _{1}=\mathbf {x} _{0}+\alpha _{0}\mathbf {p} _{0}={\begin{bmatrix}2\\1\end{bmatrix}}+{\frac {73}{331}}{\begin{bmatrix}-8\\-3\end{bmatrix}}\approx {\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}.$

Isso completa a primeira iteração, cujo resultado é uma melhor aproximação ${\textstyle \mathbf {x} _{1}}$ para a solução do sistema. Prosseguindo com o método, computamos o próximo resíduo ${\textstyle \mathbf {\mathbf {r} _{1}} }$ como

$\mathbf {r} _{1}=\mathbf {r} _{0}-\alpha _{0}\mathbf {A} \mathbf {p} _{0}={\begin{bmatrix}-8\\-3\end{bmatrix}}-{\frac {73}{331}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}\approx {\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}.$

O próximo passo no processo é calcular o escalar ${\textstyle \beta _{0}}$ que será eventualmente usado para determinar a direção de busca ${\textstyle \mathbf {\mathbf {p} _{1}} }$

$\beta _{0}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {r} _{0}^{\mathsf {T}}\mathbf {r} _{0}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-8&-3\end{bmatrix}}{\begin{bmatrix}-8\\-3\end{bmatrix}}}}=0.0088.$

Agora, usando ${\textstyle \beta _{0}}$, podemos computar ${\textstyle \mathbf {\mathbf {p} _{1}} }$ usando a relação

$\mathbf {p} _{1}=\mathbf {r} _{1}+\beta _{0}\mathbf {p} _{0}\approx {\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}+0.0088{\begin{bmatrix}-8\\-3\end{bmatrix}}={\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}.$

Computamos então o escalar ${\textstyle \alpha _{1}}$ usando a direção de busca ${\textstyle \mathbf {\mathbf {p} _{1}} }$

$\alpha _{1}={\frac {\mathbf {r} _{1}^{\mathsf {T}}\mathbf {r} _{1}}{\mathbf {p} _{1}^{\mathsf {T}}\mathbf {Ap} _{1}}}\approx {\frac {{\begin{bmatrix}-0.2810&0.7492\end{bmatrix}}{\begin{bmatrix}-0.2810\\0.7492\end{bmatrix}}}{{\begin{bmatrix}-0.3511&0.7229\end{bmatrix}}{\begin{bmatrix}4&1\\1&3\end{bmatrix}}{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}}}=0.4122.$

Finalmente, encontramos ${\textstyle \mathbf {x} _{2}}$ usando o mesmo método que aquele usado para calcular ${\textstyle \mathbf {x} _{1}}$.

$\mathbf {x} _{2}=\mathbf {x} _{1}+\alpha _{1}\mathbf {p} _{1}\approx {\begin{bmatrix}0.2356\\0.3384\end{bmatrix}}+0.4122{\begin{bmatrix}-0.3511\\0.7229\end{bmatrix}}={\begin{bmatrix}0.0909\\0.6364\end{bmatrix}}.$

O resultado, ${\textstyle \mathbf {x} _{2}}$, é uma melhor aproximação da solução do sistema que ${\textstyle \mathbf {x} _{1}}$ e ${\textstyle \mathbf {x} _{0}}$. Em aritmética exata, a solução exata do sistema é encontrada após ${\textstyle n=2}$ iterações (${\textstyle n}$ sendo a ordem do sistema).

Propriedade da terminação finita

Em aritmética exata, o número de iterações necessário não é maior que a ordem da matriz. Esse comportamento é conhecido como a propriedade de terminação finita do método do gradiente conjugado. Ela se refere à capacidade do método de chegar à solução exata de um sistema linear em um número finito de passos—no máximo igual à dimensão do sistema—quando utilizada aritmética exata. Essa propriedade decorre do fato de que, em cada iteração, o método gera um vetor resíduo que é ortogonal a todos os resíduos anteriores. Esses vetores formam um conjunto ortogonal.

Em um espaço ${\textstyle n}$-dimensional, é impossível construir mais que ${\textstyle n}$ vetores linearmente independentes e ortogonais dois a dois a menos que um deles seja o vetor nulo. Então, assim que o resíduo nulo aparecer, o método terá chegado à solução e deve terminar. Isso garante que o método do gradiente conjugado converge em no máximo ${\textstyle n}$ passos.

Aplicação a sistemas esparsos

A propriedade de terminação finita também tem implicações práticas na resolução de sistemas lineares grandes e esparsos, que surgem frequentemente em aplicações da ciência e das engenharias. Por exemplo, a discretização da equação de Laplace bidimensional ${\textstyle \mathbf {\nabla } ^{2}\mathbf {u} =0}$ usando diferenças finitas em uma malha uniforme leva ao sistema linear esparso ${\textstyle \mathbf {A} \mathbf {x} =\mathbf {b} }$, onde ${\textstyle \mathbf {A} }$ é simétrica e positiva definida.

O uso de uma malha interior de tamanho ${\textstyle 5\times 5}$ produz um sistema ${\textstyle 25\times 25}$, e a matriz dos coeficientes ${\textstyle \mathbf {A} }$ tem um padrão de estêncil de cinco pontos. Cada linha de ${\textstyle \mathbf {A} }$ contém no máximo cinco entradas não nulas correspondentes ao ponto central e a seus vizinhos imediatos. Por exemplo, a matriz gerada de tal malha pode ter a forma

$\mathbf {A} ={\begin{bmatrix}4&-1&0&\cdots &-1&0&\cdots \\-1&4&-1&\cdots &0&0&\cdots \\0&-1&4&-1&0&0&\cdots \\\vdots &\vdots &\ddots &\ddots &\ddots &\vdots \\-1&0&\cdots &-1&4&-1&\cdots \\0&0&\cdots &0&-1&4&\cdots \\\vdots &\vdots &\cdots &\cdots &\cdots &\ddots \end{bmatrix}}$

Embora o sistema tenha dimensão 25, o método do gradiente conjugado deve convergir em no máximo 25 iterações em aritmética exata. Na prática, a convergência frequentemente ocorre em menos iterações devido a propriedades espectrais da matriz. Essa eficiência faz do método do gradiente conjugado um método iterativo atraente para resolver sistemas lineares de larga escala oriundos de equações diferenciais parciais, como aqueles encontrados em condução do calor, dinâmica dos fluidos e eletrostática.

Taxa de convergência

Embora o método do gradiente conjugado possa ser interpretado como um método direto em aritmética exata, sua utilização prática é normalmente analisada como a de um método iterativo. Em aritmética exata, o método produz a solução exata após um número finito de iterações, não maior que a dimensão da matriz. Em computadores reais, entretanto, erros de arredondamento impedem que essa propriedade seja observada de forma exata.

Como método iterativo, o gradiente conjugado melhora a aproximação da solução na norma de energia. A velocidade dessa convergência depende das propriedades espectrais da matriz do sistema. Um parâmetro importante é o número de condicionamento da matriz $\mathbf {A}$. Para uma matriz simétrica positiva definida, ele pode ser escrito como

$\kappa (\mathbf {A} )={\frac {\lambda _{\max }(\mathbf {A} )}{\lambda _{\min }(\mathbf {A} )}},$

onde $\lambda _{\max }(\mathbf {A} )$ e $\lambda _{\min }(\mathbf {A} )$ são, respectivamente, o maior e o menor autovalor de $\mathbf {A}$. Em geral, quanto maior é $\kappa (\mathbf {A} )$, mais lenta tende a ser a convergência do método.

O número de condicionamento, porém, não descreve sozinho todo o comportamento do método. A distribuição dos autovalores também influencia a convergência. Matrizes com autovalores agrupados podem apresentar convergência mais rápida do que a sugerida apenas pelo valor de $\kappa (\mathbf {A} )$, enquanto certas distribuições espectrais podem causar fases de estagnação no erro iterativo.

Teorema de convergência

Define-se o subconjunto de polinômios

$\Pi _{k}^{\ast }:=\{p\in \Pi _{k}:p(0)=1\},$

onde $\Pi _{k}$ é o conjunto dos polinômios de grau máximo $k$.

Sejam $(\mathbf {x} _{k})_{k}$ as aproximações iterativas da solução exata $x_{\ast }$ do sistema $\mathbf {A} \mathbf {x} =\mathbf {b}$, e defina-se o erro por

$\mathbf {e} _{k}:=\mathbf {x} _{k}-\mathbf {x} _{\ast }.$

A taxa de convergência do método do gradiente conjugado pode ser limitada por

$\|\mathbf {e} _{k}\|_{\mathbf {A} }=\min _{p\in \Pi _{k}^{\ast }}\|p(\mathbf {A} )\mathbf {e} _{0}\|_{\mathbf {A} }\leq \min _{p\in \Pi _{k}^{\ast }}\max _{\lambda \in \sigma (\mathbf {A} )}|p(\lambda )|\|\mathbf {e} _{0}\|_{\mathbf {A} }$

$\leq 2\left({\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\right)^{k}\|\mathbf {e} _{0}\|_{\mathbf {A} }\leq 2\exp \left({\frac {-2k}{\sqrt {\kappa (\mathbf {A} )}}}\right)\|\mathbf {e} _{0}\|_{\mathbf {A} },$

onde $\sigma (\mathbf {A} )$ denota o espectro de $\mathbf {A}$, $\kappa (\mathbf {A} )$ denota o número de condicionamento da matriz, e

$\|v\|_{\mathbf {A} }={\sqrt {v^{\intercal }\mathbf {A} v}}$

é a norma de energia associada a $\mathbf {A}$.

Essa estimativa mostra que

$k={\frac {1}{2}}{\sqrt {\kappa (\mathbf {A} )}}\log \left(\|\mathbf {e} _{0}\|_{\mathbf {A} }\varepsilon ^{-1}\right)$

iterações são suficientes para reduzir o erro a $2\varepsilon$, para qualquer $\varepsilon >0$.

No limite em que $\kappa (\mathbf {A} )$ tende ao infinito, tem-se

${\frac {{\sqrt {\kappa (\mathbf {A} )}}-1}{{\sqrt {\kappa (\mathbf {A} )}}+1}}\approx 1-{\frac {2}{\sqrt {\kappa (\mathbf {A} )}}},\qquad \kappa (\mathbf {A} )\gg 1.$

Essa relação indica uma taxa de convergência mais rápida do que a de métodos iterativos clássicos como Jacobi ou Gauss-Seidel, cujas estimativas tipicamente escalam como

$1-{\frac {2}{\kappa (\mathbf {A} )}}.$

O teorema de convergência não assume erros de arredondamento. Ainda assim, essa cota costuma descrever bem o comportamento prático do método em muitas situações.

Convergência prática

Quando o método é inicializado de maneira aleatória, a primeira fase de iterações costuma ser a mais rápida, pois o erro é reduzido dentro do subespaço de Krylov inicialmente gerado, que pode refletir um número de condicionamento efetivo menor. A segunda fase da convergência costuma ser bem descrita pela cota teórica envolvendo ${\sqrt {\kappa (\mathbf {A} )}}$, mas pode apresentar comportamento superlinear, dependendo da distribuição dos autovalores da matriz $\mathbf {A}$ e da distribuição espectral do erro.

Na fase final, a menor precisão atingível pela aritmética de ponto flutuante é alcançada, e a convergência pode estagnar ou até começar a divergir. Em aplicações típicas de computação científica com aritmética de dupla precisão, o método do gradiente conjugado usa um critério de parada baseado em tolerância, de modo que as iterações geralmente são encerradas durante a primeira ou a segunda fase da convergência.

Precondicionamento para o método do gradiente conjugado

Na maioria dos casos, o precondicionamento é usado para acelerar a convergência do método do gradiente conjugado. A ideia é substituir o sistema linear original por um sistema equivalente, mas com propriedades espectrais mais favoráveis. Para isso, escolhe-se uma matriz $\mathbf {M}$, chamada precondicionador, que deve aproximar $\mathbf {A}$, mas ser mais barata de aplicar do que resolver diretamente o sistema original.

Se $\mathbf {M} ^{-1}$ é simétrica positiva definida e $\mathbf {M} ^{-1}\mathbf {A}$ possui número de condicionamento melhor que $\mathbf {A}$, pode-se usar o método do gradiente conjugado precondicionado. Uma forma comum do algoritmo é:

$\mathbf {r} _{0}:=\mathbf {b} -A\mathbf {x} _{0}$

resolva $\mathbf {M} \mathbf {z} _{0}:=\mathbf {r} _{0}$

$\mathbf {p} _{0}:=\mathbf {z} _{0}$

$k:=0$

Enquanto o critério de parada não for satisfeito, faça:

$\alpha _{k}:={\frac {\mathbf {r} _{k}^{\intercal }\mathbf {z} _{k}}{\mathbf {p} _{k}^{\intercal }\mathbf {A} \mathbf {p} _{k}}}$

$\mathbf {x} _{k+1}:=\mathbf {x} _{k}+\alpha _{k}\mathbf {p} _{k}$

$\mathbf {r} _{k+1}:=\mathbf {r} _{k}-\alpha _{k}A\mathbf {p} _{k}$

Se $\mathbf {r} _{k+1}$ for suficientemente pequeno, então pare o ciclo.

resolva $\mathbf {M} \mathbf {z} _{k+1}:=\mathbf {r} _{k+1}$

$\beta _{k}:={\frac {\mathbf {r} _{k+1}^{\intercal }\mathbf {z} _{k+1}}{\mathbf {r} _{k}^{\intercal }\mathbf {z} _{k}}}$

$\mathbf {p} _{k+1}:=\mathbf {z} _{k+1}+\beta _{k}\mathbf {p} _{k}$

$k:=k+1$

Fim Enquanto.

O resultado é $\mathbf {x} _{k+1}$.

A formulação acima é equivalente à aplicação do método do gradiente conjugado usual ao sistema precondicionado

$\mathbf {E} ^{-1}\mathbf {A} (\mathbf {E} ^{-1})^{\intercal }\mathbf {\hat {x}} =\mathbf {E} ^{-1}\mathbf {b} ,$

onde

$\mathbf {E} \mathbf {E} ^{\intercal }=\mathbf {M}$

e

$\mathbf {\hat {x}} =\mathbf {E} ^{\intercal }\mathbf {x} .$

A decomposição de Cholesky do precondicionador preserva a simetria e a positividade definida do sistema transformado. Entretanto, essa decomposição não precisa ser calculada explicitamente, basta saber aplicar o efeito de $\mathbf {M} ^{-1}$, isto é, resolver sistemas lineares auxiliares envolvendo $\mathbf {M}$.

A matriz precondicionadora $\mathbf {M}$ deve ser simétrica positiva definida e fixa, isto é, não deve mudar de uma iteração para outra. Se essas hipóteses forem violadas, o comportamento do método do gradiente conjugado precondicionado pode se tornar imprevisível.

Decomposição Cholesky incompleta

Um exemplo comum de precondicionador para matrizes simétricas positivas definidas é a decomposição de Cholesky incompleta. Na decomposição de Cholesky usual, uma matriz simétrica positiva definida é fatorada como

$\mathbf {A} =\mathbf {L} \mathbf {L} ^{\intercal },$

onde $\mathbf {L}$ é uma matriz triangular inferior. Para matrizes grandes e esparsas, entretanto, a fatoração completa pode introduzir muitos elementos não nulos adicionais, efeito denominado de fill-in (preenchimento), aumentando o custo de armazenamento na memória e de cálculo. A decomposição de Cholesky incompleta busca uma aproximação

$\mathbf {A} \approx \mathbf {\tilde {L}} \mathbf {\tilde {L}} ^{\intercal },$

mantendo parte da estrutura esparsa da matriz original. O precondicionador correspondente é

$\mathbf {M} =\mathbf {\tilde {L}} \mathbf {\tilde {L}} ^{\intercal }.$

Na aplicação prática do precondicionador, não se calcula explicitamente $\mathbf {M} ^{-1}$. Em vez disso, para resolver

$\mathbf {M} \mathbf {z} =\mathbf {r} ,$

usa-se a fatoração incompleta. Como

$\mathbf {M} =\mathbf {\tilde {L}} \mathbf {\tilde {L}} ^{\intercal },$

tem-se

$\mathbf {\tilde {L}} \mathbf {\tilde {L}} ^{\intercal }\mathbf {z} =\mathbf {r} .$

Tomando um vetor intermediário $\mathbf {a}$, resolve-se primeiro

$\mathbf {\tilde {L}} \mathbf {a} =\mathbf {r}$

por substituição para frente (forward substitution). Em seguida, resolve-se

$\mathbf {\tilde {L}} ^{\intercal }\mathbf {z} =\mathbf {a}$

por substituição para trás (backward substitution). Desse modo, obtém-se $\mathbf {z} =\mathbf {M} ^{-1}\mathbf {r}$ sem calcular explicitamente $\mathbf {M} ^{-1}$.

A qualidade da decomposição incompleta depende da estratégia usada para descartar elementos durante a fatoração (drop tolerance). Manter mais elementos tende a produzir um precondicionador melhor, mas aumenta o custo de armazenamento e de aplicação. Descartar mais elementos reduz o custo por iteração, mas pode diminuir o efeito do precondicionamento.

Precondicionador SSOR

Outro precondicionador usado em sistemas simétricos positivos definidos é o SSOR (symmetric successive over-relaxation). Ele é derivado do método SOR e usa a decomposição da matriz $\mathbf {A}$ em partes diagonal, triangular inferior e triangular superior.

Se

$\mathbf {A} =\mathbf {D} +\mathbf {L} +\mathbf {L} ^{\intercal },$

onde $\mathbf {D}$ é a parte diagonal de $\mathbf {A}$ e $\mathbf {L}$ é a parte triangular inferior estrita, uma forma comum do precondicionador SSOR é

$\mathbf {M} _{\mathrm {SSOR} }={\frac {1}{\omega (2-\omega )}}(\mathbf {D} +\omega \mathbf {L} )\mathbf {D} ^{-1}(\mathbf {D} +\omega \mathbf {L} ^{\intercal }),$

com ${\textstyle 0<\omega <2.}$

O parâmetro $\omega$ é o parâmetro de relaxação. Sua escolha influencia a eficiência do precondicionador, onde valores inadequados podem reduzir o benefício do método, enquanto uma escolha apropriada pode acelerar a convergência.

Em comparação com a decomposição de Cholesky incompleta, o SSOR é conceitualmente simples e pode ser construído diretamente a partir das partes da matriz. No entanto, sua eficiência depende do problema considerado e da escolha do parâmetro $\omega$.

Comparação dos métodos iterativos

Para ilustrar a performance de diferentes métodos iterativos, considere o sistema linear $\mathbf {A} \mathbf {x} =\mathbf {b}$, com

$\mathbf {A} ={\begin{bmatrix}0.2&0.1&1&1&0\\0.1&4&-1&1&-1\\1&-1&60&0&-2\\1&1&0&8&4\\0&-1&-2&4&700\end{bmatrix}}\quad {\text{e}}\quad \mathbf {b} ={\begin{bmatrix}1\\2\\3\\4\\5\end{bmatrix}}$

cuja solução é

$\mathbf {x} _{\ast }=(7.859713071,0.4229264082,-0.07359223906,-0.5406430164,0.01062616286)^{\intercal }.$

A tabela abaixo lista os resultados obtidos ao utilizar os métodos iterativos de Jacobi, Gauss-Seidel e SOR (com $\omega =1.25$) aplicados ao sistema com $\mathbf {A}$ com uma tolerância de 0.01, bem como aqueles de quando o CG é aplicado tanto em sua forma não precondicionada quanto utilizando uma matriz de precondicionamento. O método do gradiente conjugado precondicionado não apenas fornece as aproximações mais precisas, mas também utiliza o menor número de iterações.

Método	Número de Iterações	$\mathbf {x} _{i}$	$\|\mathbf {x} _{\ast }-\mathbf {x} _{i}\|_{\infty }$
Jacobi	49	$(7.86277141,0.42320802,-0.07348669,-0.53975964,0.01062847)^{\intercal }$	0.00305834
Gauss-Seidel	15	$(7.83525748,0.42257868,-0.07319124,-0.53753055,0.01060903)^{\intercal }$	0.02445559
SOR ($\omega =1.25$)	7	$(7.85152706,0.42277371,-0.07348303,-0.53978369,0.01062286)^{\intercal }$	0.00818607
Gradiente Conjugado	5	$(7.85341523,0.42298677,-0.07347963,-0.53987920,0.008628916)^{\intercal }$	0.00629785
Gradiente Conjugado (Precondicionado)	4	$(7.85968827,0.42288329,-0.07359878,-0.54063200,0.01064344)^{\intercal }$	0.00009312

Gradiente conjugado aplicado às equações normais

O método do gradiente conjugado pode ser aplicado a uma matriz retangular arbitrária ${\textstyle \mathbf {A} _{m\times n}}$ através das equações normais ${\textstyle \mathbf {A} ^{\intercal }\mathbf {A} \mathbf {x} =\mathbf {A} ^{\intercal }\mathbf {b} }$, pois ${\textstyle \mathbf {A} ^{\intercal }\mathbf {A} }$ é positiva semidefinida (i.e., ${\textstyle \mathbf {x} ^{\intercal }\mathbf {A} ^{\intercal }\mathbf {A} \mathbf {x} >0}$ para todo ${\textstyle \mathbf {x} \neq 0}$) para qualquer matriz ${\textstyle \mathbf {A} }$. Isso resulta no método CGNR (do inglês, conjugate gradient normal equation residual).

Como um método iterativo, não é necessário formar as matrizes ${\textstyle \mathbf {A} ^{\intercal }\mathbf {A} }$ explicitamente, mas apenas computar multiplicações matriz-vetor e matriz transposta-vetor. Assim, o CGNR é particularmente útil quando ${\textstyle \mathbf {A} }$ é uma matriz esparsa, pois essas operações são geralmente realizadas de forma extremamente eficiente. A desvantagem de resolver as equações normais ${\textstyle \mathbf {A} ^{\intercal }\mathbf {A} }$ é que o número de condicionamento ${\textstyle \kappa (\mathbf {A} ^{\intercal }\mathbf {A} )}$ é igual a ${\textstyle \kappa ^{2}(\mathbf {A} )}$ e então a taxa de convergência do CGNR pode ser lenta e a acurácia da solução pode ser mais sensível a erros de arredondamento. Encontrar um precondicionador adequado torna-se fundamental para utilizar o método CGNR.

Leitura adicional

O método do gradiente conjugado foi originalmente proposto em

• Hestenes, Magnus R.; Stiefel, Eduard (dezembro de 1952). «Methods of Conjugate Gradients for Solving Linear Systems» (PDF). Journal of Research of the National Bureau of Standards. 49 (6). Cópia arquivada (PDF) em 29 de maio de 2026 

Descrições do método podem ser encontradas nos seguintes livros texto:

• Kendell A. Atkinson (1988), An introduction to numerical analysis (2nd ed.), Section 8.9, John Wiley and Sons. ISBN 0-471-50023-2.
• Mordecai Avriel (2003). Nonlinear Programming: Analysis and Methods. Dover Publishing. ISBN 0-486-43227-0.
• Gene H. Golub and Charles F. Van Loan, Matrix computations (3rd ed.), Chapter 10, Johns Hopkins University Press. ISBN 0-8018-5414-8.

Ligações externas

• «Método do gradiente conjugado». por Nadir Soualem. 
• «Método do gradiente conjugado com precondicionamento». por Nadir Soualem. 
• «An Introduction to the Conjugate Gradient Method Without the Agonizing Pain» (PDF). por Jonathan Richard Shewchuk. 
• «Métodos iterativos para sistemas lineares esparsos». por Yousef Saad 
• «LSQR: Equações esparsas e mínimos quadrados». por Christopher Paige e Michael Saunders.