Resumo: Método do Gradiente Conjugado

Prof. Doherty Andrade

Texto baseado em Burden & Faires.

O método do gradiente conjugado de Hestenes e Stiefel [HS] foi originalmente desenvolvido como um método direto criado para resolver um sistema linear definido positivo \( n \times n \). Como um método direto, é geralmente inferior à eliminação de Gauss com pivotamento, já que ambos os métodos exigem \( n \) passos para determinar uma solução, e os passos do método do gradiente conjugado são computacionalmente mais caros que os da eliminação de Gauss.

Todavia, o método do gradiente conjugado é muito útil quando empregado como um método de aproximação iterativa para resolver grandes sistemas esparsos com elementos não nulos ocorrendo em padrões previsíveis. Esses problemas freqüentemente surgem na solução de problemas de contorno. Quando a matriz tiver sido pré-condicionada para tornar os cálculos mais eficientes, bons resultados são obtidos em cerca de apenas \( \sqrt{n} \) passos. Empregado dessa maneira, o método é preferido ao da eliminação de Gauss e aos métodos iterativos anteriormente discutidos.

Magnus Hestenes (1906-1991) e Edward Stiefel publicaram o artigo original sobre o método do gradiente conjugado em 1952, enquanto trabalhavam no Instituto de Análise Numérica no campus da UCLA.

Nesta seção, supomos que a matriz \( A \) seja definida positiva. Utilizaremos a notação de produto interno.

\[ \langle x, y \rangle = x' y, \tag{7.25} \]

em que \( x \) e \( y \) são vetores \( n \) dimensionais. Também precisaremos de resultados adicionais da álgebra linear.

O próximo resultado segue facilmente a partir das propriedades da transposta. (Veja o Exercício 12.)

Quando \( A \) for definida positiva, \( \langle x, Ax \rangle = x'Ax > 0 \) a menos que \( x = 0 \). Também, como \( A \) é simétrica, temos \( x'Ay = x'A'y = (Ax)'y \), e assim, além dos resultados no Teorema 7.30, temos para cada \( x \) e \( y \),

\[ \langle x, Ay \rangle = \langle Ax, y \rangle. \tag{7.26} \]

O resultado seguinte é uma ferramenta básica no desenvolvimento do método do gradiente conjugado.

Demonstração Sejam \( x \) e \( v \neq 0 \) vetores fixos e \( t \) uma variável real. Temos

\[ \begin{aligned} g(x + tv) &= \langle x + tv, Ax + tAv \rangle - 2\langle x + tv, b \rangle \\ &= \langle x, Ax \rangle + t\langle v, Ax \rangle + t\langle x, Av \rangle + t^2\langle v, Av \rangle - 2\langle x, b \rangle - 2t\langle v, b \rangle \\ &= \langle x, Ax \rangle - 2\langle x, b \rangle + 2t\langle v, Ax \rangle - 2t\langle v, b \rangle + t^2\langle v, Av \rangle, \end{aligned} \]

e assim,

\[ g(x + tv) = g(x) + 2t\langle v, Ax - b \rangle + t^2\langle v, Av \rangle. \tag{7.27} \]

Como \( x \) e \( v \) são fixos, podemos definir a função quadrática \( h \) em \( t \) por meio de

\[ h(t) = g(x + tv). \]

Então, \( h \) assume um valor mínimo quando \( h'(t) = 0 \), pois o seu coeficiente de \( t^2, \langle v, Av \rangle, \delta \) possui. Já que

\[ h'(t) = 2\langle v, Ax - b \rangle + 2t\langle v, Av \rangle, \]

o mínimo ocorre quando

\[ \hat{t} = -\frac{\langle v, Ax - b \rangle}{\langle v, Av \rangle} = \frac{\langle v, b - Ax \rangle}{\langle v, Av \rangle}, \]

e, da Equação (7.27),

\[ \begin{aligned} h(\hat{t}) &= g(x) - 2\frac{\langle v, b - Ax \rangle}{\langle v, Av \rangle}\langle v, b - Ax \rangle + \left( \frac{\langle v, b - Ax \rangle}{\langle v, Av \rangle} \right)^2\langle v, Av \rangle \\ &= g(x) - \frac{\langle v, b - Ax \rangle^2}{\langle v, Av \rangle}. \end{aligned} \]

Assim, para qualquer vetor \( v \neq 0 \), temos \( g(x + \hat{t}v) < g(x) \) a menos que \( \langle v, b - Ax \rangle = 0 \), caso no qual \( g(x) = g(x + \hat{t}v) \). Esse é o resultado básico de que precisamos para demonstrar o Teorema 7.31.

Suponha que \( x^* \) satisfaça a \( Ax^* = b \). Então, \( \langle v, b - Ax^* \rangle = 0 \) para qualquer vetor \( v \), e \( g(x) \) não pode se tornar menor que \( g(x^*) \). Assim, \( x^* \) minimiza \( g \).

Por outro lado, suponha que \( x^* \) seja um vetor que minimiza \( g \). Então, para qualquer vetor \( v \) temos \( g(x^* + \hat{t}v) \geq g(x^*) \). Assim, \( \langle v, b - Ax^* \rangle = 0 \). Isso implica que \( b - Ax^* = 0 \) e, conseqüentemente, que \( Ax^* = b \).

Para começar o método do gradiente conjugado, escolhemos \( x \), uma solução aproximada para \( Ax^* = b \), e \( v \neq 0 \), o que fornece uma direção de busca para se afastar de \( x \) para melhorar a aproximação. Seja \( r = b - Ax \) o vetor resíduo associado a \( x \) e

\[ t = \frac{\langle v, b - Ax \rangle}{\langle v, Av \rangle} = \frac{\langle v, r \rangle}{\langle v, Av \rangle}. \]

Se \( r \neq 0 \) e se \( v \) e \( r \) não forem ortogonais, então \( x + tv \) fornece um valor menor para \( g \) que \( g(x) \) está presumivelmente mais próximo de \( x^* \) do que \( x \). Isso sugere o método seguinte.

Seja \( x^{(0)} \) uma aproximação inicial de \( x^* \), e seja \( v^{(1)} \neq 0 \), um sentido inicial de busca. Para \( k = 1, 2, 3, \ldots \), calculamos

\[ \begin{aligned} t_k &= \frac{\langle v^{(k)}, b - A x^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle}, \\ x^{(k)} &= x^{(k-1)} + t_k v^{(k)}. \end{aligned} \]

e escolhemos uma nova direção de busca \( v^{(k+1)} \). O objetivo é fazer esta escolha de modo que a seqüência de aproximações \( \{x^{(k)}\} \) convirja rapidamente para \( x^* \).

Para escolher a direção de busca, consideramos g como uma função das componentes de \( x = (x_1, x_2, \ldots, x_n) \). Assim,

\[ g(x_1, x_2, \ldots, x_n) = \langle x, Ax \rangle - 2\langle x, b \rangle = \sum_{i=1}^n \sum_{j=1}^n a_{ij} x_i x_j - 2 \sum_{i=1}^n x_i b_i. \]

Tornando a derivada parcial em relação às variáveis componentes \( x_k \), temos

\[ \frac{\partial g}{\partial x_k} (x) = 2 \sum_{i=1}^n a_{ij} x_i - 2b_k. \]

Assim, o gradiente de g é

\[ \nabla g(x) = \left( \frac{\partial g}{\partial x_1} (x), \frac{\partial g}{\partial x_2} (x), \ldots, \frac{\partial g}{\partial x_n} (x) \right)' = 2(A x - b) = -2r, \]

em que o vetor r é o vetor resíduo para x.

Do cálculo de várias variáveis, sabemos que a direção de maior diminuição no valor de \( g(x) \) é a direção dada por \( -\nabla g(x) \); isto é, na direção do resíduo \(r\). O método que escolhe

\[ v^{(k+1)} = r^{(k)} = b - A x^{(k)} \]

é denominado o método do declive máximo. Embora venhamos a ver na Seção 10.4 que este método tem méritos para sistemas não-lineares e problemas de otimização, ele não é utilizado para sistemas lineares por causa de convergência lenta.

Uma abordagem alternativa usa um conjunto de vetores direção não nulos \( \{v^{(1)}, \ldots, v^{(n)}\} \) que satisfaçam a

\[ \langle v^{(i)}, A v^{(j)} \rangle = 0 \quad \text{se} \quad i \neq j. \]

Isso é chamado uma condição de A-ortogonalidade e diz-se que o conjunto de vetores \( \{v^{(1)}, \ldots, v^{(n)}\} \) é A-ortogonal. Não é difícil de mostrar que um conjunto de vetores A-ortogonais associados à matriz definida positiva A é linearmente independente. (Consultar o Exercício 13(a)). Esse conjunto de direções de busca dá

\[ \begin{aligned} t_k &= \frac{\langle v^{(k)}, b - A x^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle} = \frac{\langle v^{(k)}, r^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle} \\ \mbox{ e } x^{(k)} &= x^{(k-1)} + t_k v^{(k)}. \end{aligned} \]

O teorema seguinte mostra que esta escolha de direções de busca fornece convergência em, no máximo, \(n\) passos, de modo que, como um método direto, produz a solução exata, supondo que a aritmética seja exata.

Demonstração Como, para cada \( k = 1, 2, \ldots, n \),

\[ x^{(k)} = x^{(k-1)} + t_k v^{(k)}, \]

temos

\[ \begin{aligned} A x^{(n)} &= A x^{(n-1)} + t_n A v^{(n)} \\ &= (A x^{(n-2)} + t_{n-1} A v^{(n-1)}) + t_n A v^{(n)} \\ &= A x^{(0)} + t_1 A v^{(1)} + t_2 A v^{(2)} + \cdots + t_n A v^{(n)} \end{aligned} \]

e, subtraindo \(b\) deste resultado, vem

\[ A x^{(n)} - b = A x^{(0)} - b + t_1 A v^{(1)} + t_2 A v^{(2)} + \cdots + t_n A v^{(n)} \]

Agora, tomamos o produto interno de ambos os lados com o vetor \( v^{(k)} \) e usamos as propriedades de produto interno e o fato de \( A \) ser simétrica para obter

\[ \begin{aligned} \langle A x^{(n)} - b, v^{(k)} \rangle &= \langle A x^{(0)} - b, v^{(k)} \rangle + t_1 \langle A v^{(1)}, v^{(k)} \rangle + \cdots + t_n \langle A v^{(n)}, v^{(k)} \rangle \\ &= \langle A x^{(0)} - b, v^{(k)} \rangle + t_1 \langle v^{(1)}, A v^{(k)} \rangle + \cdots + t_n \langle v^{(n)}, A v^{(k)} \rangle \end{aligned} \]

A propriedade de \( A \)-ortogonalidade dá, para cada \( k \),

\[ \langle A x^{(n)} - b, v^{(k)} \rangle = \langle A x^{(0)} - b, v^{(k)} \rangle + t_k \langle v^{(k)}, A v^{(k)} \rangle \tag{7.28} \]

Todavia,

\[ t_k = \frac{\langle v^{(k)}, b - A x^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle} \]

de modo que,

\[ \begin{aligned} t_k \langle v^{(k)}, A v^{(k)} \rangle &= \langle v^{(k)}, b - A x^{(k-1)} \rangle \\ &= \langle v^{(k)}, b - A x^{(0)} + A x^{(0)} - A x^{(1)} + \cdots - A x^{(k-2)} + A x^{(k-2)} - A x^{(k-1)} \rangle \\ &= \langle v^{(k)}, b - A x^{(0)} \rangle + \langle v^{(k)}, A x^{(0)} - A x^{(1)} \rangle + \cdots + \langle v^{(k)}, A x^{(k-2)} - A x^{(k-1)} \rangle \end{aligned} \]

Mas, para qualquer \( i \),

\[ x^{(i)} = x^{(i-1)} + t_i v^{(i)} \quad \text{e} \quad A x^{(i)} = A x^{(i-1)} + t_i A v^{(i)} \]

e então,

\[ A x^{(i-1)} - A x^{(i)} = - t_i A v^{(i)} \]

Assim,

\[ t_k \langle v^{(k)}, A v^{(k)} \rangle = \langle v^{(k)}, b - A x^{(0)} \rangle - t_1 \langle v^{(k)}, A v^{(1)} \rangle - \cdots - t_{k-1} \langle v^{(k)}, A v^{(k-1)} \rangle \]

Devido à \( A \)-ortogonalidade, \( \langle v^{(k)}, A v^{(i)} \rangle = 0 \), para \( i \neq k \), de modo que

\[ t_k \langle v^{(k)}, A v^{(k)} \rangle = \langle v^{(k)}, b - A x^{(0)} \rangle \]

Da Equação (7.28),

\[ \begin{aligned} \langle A x^{(n)} - b, v^{(k)} \rangle &= \langle A x^{(0)} - b, v^{(k)} \rangle + \langle v^{(k)}, b - A x^{(0)} \rangle \\ &= \langle A x^{(0)} - b, v^{(k)} \rangle + \langle v^{(k)}, b - A x^{(0)} \rangle \\ &= \langle A x^{(0)} - b, v^{(k)} \rangle - \langle A x^{(0)} - b, v^{(k)} \rangle \\ &= 0 \end{aligned} \]

O vetor \( A x^{(n)} - b \) é ortogonal ao conjunto de vetores \( A \)-ortogonais \( \{v^{(1)}, \cdots, v^{(n)}\} \). Disso, segue que \( A x^{(n)} - b = 0 \). (Consulte o Exercício 13(b).)

Antes de discutir como determinar o conjunto A-ortogonal, continuaremos o desenvolvimento. A utilização de um conjunto A-ortogonal \( \{v^{(1)}, \ldots , v^{(n)}\} \) de vetores direção fornece o que é denominado um método de direção conjugada. O teorema seguinte mostra a ortogonalidade dos vetores residuais \( r^{(k)} \) e dos vetores direção \( v^{(j)} \). Uma demonstração desse resultado, usando indução matemática, é considerada no Exercício 14.

O método do gradiente conjugado de Hestenes e Stiefel escolhe as direções de busca \( \{v^{(k)}\} \) durante o processo iterativo, de modo que os vetores residuais \( \{r^{(k)}\} \) sejam mutuamente ortogonais. Para construir os vetores direção \( \{v^{(1)}, v^{(2)}, \ldots \} \) e as aproximações \( \{x^{(1)}, x^{(2)}, \ldots \} \), começamos com uma aproximação inicial \( x^{(0)} \) e usamos a direção de máximo declive \( r^{(0)} = b - Ax^{(0)} \) como a primeira direção de busca \( v^{(1)} \).

Suponha que as direções conjugadas \( v^{(1)}, \ldots , v^{(k-1)} \) e as aproximações \( x^{(1)}, \ldots , x^{(k-1)} \) foram calculadas com

\[ x^{(i)} = x^{(i-1)} + t_i v^{(i)}, \]

em que

\[ \langle v^{(i)}, A v^{(j)} \rangle = 0 \quad \text{e} \quad \langle r^{(i)}, r^{(j)} \rangle = 0 \quad \text{para} \quad i \neq j. \]

Se \( x^{(k-1)} \) é a solução para \( Ax = b \), terminamos. Caso contrário, \( r^{(k-1)} = b - A x^{(k-1)} \neq 0 \) e o Teorema 7.33 implica que \( \langle r^{(k-1)}, v^{(i)} \rangle = 0 \) para \( i = 1, 2, \ldots , k-1 \). Então usamos \( r^{(k-1)} \) para gerar \( v^{(k)} \), fazendo

\[ v^{(k)} = r^{(k-1)} + s_{k-1} v^{(k-1)}. \]

Escolheremos \( s_{k-1} \) de modo que

\[ \langle v^{(k-1)}, A v^{(k)} \rangle = 0. \]

Como

\[ A v^{(k)} = A r^{(k-1)} + s_{k-1} A v^{(k-1)} \]

temos \( \langle v^{(k-1)}, A v^{(k)} \rangle = 0 \) quando

\[ s_{k-1} = -\frac{\langle v^{(k-1)}, A r^{(k-1)} \rangle}{\langle v^{(k-1)}, A v^{(k-1)} \rangle}. \]

Também pode ser mostrado que com esta escolha de \( s_{k-1} \) temos \( \langle v^{(i)}, A v^{(k)} \rangle = 0 \), para cada \( i = 1, \ldots, k-2 \) (consultar [Lu, pág. 245]). Assim, \( \{v^{(1)}, \ldots , v^{(k)}\} \) é um conjunto A-ortogonal.

\[ \begin{aligned} t_k &= \frac{\langle v^{(k)}, r^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle} = \frac{\langle r^{(k-1)} + s_{k-1} v^{(k-1)}, r^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle} \\ &= \frac{\langle r^{(k-1)}, r^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle} + s_{k-1} \frac{\langle v^{(k-1)}, r^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle}. \end{aligned} \]

Pelo Teorema 7.33, temos \( \langle v^{(k-1)}, r^{(k-1)} \rangle = 0 \), e então

\[ t_k = \frac{\langle r^{(k-1)}, r^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle}. \tag{7.29} \]

Assim,

\[ x^{(k)} = x^{(k-1)} + t_k v^{(k)}. \]

Para calcular \( r^{(k)} \), multiplicamos por \( A \) e subtraímos \( b \) para obter

\[ A x^{(k)} - b = A x^{(k-1)} - b + t_k A v^{(k)} \quad \text{ou} \quad r^{(k)} = r^{(k-1)} - t_k A v^{(k)}. \]

Assim,

\[ \langle r^{(k)}, r^{(k)} \rangle = \langle r^{(k-1)}, r^{(k)} \rangle - t_k \langle A v^{(k)}, r^{(k)} \rangle = -t_k \langle r^{(k)}, A v^{(k)} \rangle. \]

E mais, da Equação (7.29),

\[ \langle r^{(k-1)}, r^{(k-1)} \rangle = t_k \langle v^{(k)}, A v^{(k)} \rangle, \]

de modo que

\[ \begin{aligned} s_k &= -\frac{\langle v^{(k)}, A r^{(k)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle} = -\frac{\langle r^{(k)}, A v^{(k)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle} \\ &= \frac{(1/t_k)\langle r^{(k)}, r^{(k)} \rangle}{(1/t_k)\langle r^{(k-1)}, r^{(k-1)} \rangle} = \frac{\langle r^{(k)}, r^{(k)} \rangle}{\langle r^{(k-1)}, r^{(k-1)} \rangle}. \end{aligned} \]

Em suma, temos as fórmulas:

Em vez de apresentar um algoritmo para o método do gradiente conjugado usando essas fórmulas, estendemos o método para incluir pré-condicionamento.

Método do Gradiente Conjugado Précondicionado

Para aplicar o método a um sistema mais bem condicionado, queremos escolher uma matriz condicionante não-singular \( C \) tal que

\[ \hat{A} = C^{-1} A (C^{-1})^t \]

seja mais bem condicionada. Para simplificar a notação, utilizaremos a matriz \( C^{-t} \) para nos referir a \( (C^{-1})^t \).

Considere o sistema linear

\[ \hat{A} \hat{x} = \hat{b}, \]

em que \( \hat{x} = C^t x \) e \( \hat{b} = C^{-1} b \). Então

\[ \hat{A} \hat{x} = (C^{-1} A C^{-t}) (C^t x) = C^{-1} A x. \]

Assim, poderíamos resolver \( \hat{A} \hat{x} = \hat{b} \) e então obter \( x \) multiplicando por \( C^{-t} \). Todavia, em vez de reescrever a Equação (7.30) utilizando \( \hat{r}^{(k)} \), \( \hat{v}^{(k)} \), \( \hat{t}_k \), \( \hat{x}^{(k)} \) e \( \hat{s}_k \), incorporamos implicitamente o pré-condicionamento.

Como

\[ \hat{x}^{(k)} = C^t x^{(k)}, \]

temos

\[ \begin{aligned} \hat{r}^{(k)} &= \hat{b} - \hat{A} \hat{x}^{(k)} = C^{-1} b - (C^{-1} A C^{-t}) C^t x^{(k)} \\ &= C^{-1} (b - A x^{(k)}) = C^{-1} r^{(k)}. \end{aligned} \]

Sejam \( \hat{v}^{(k)} = C^t v^{(k)} \) e \( w^{(k)} = C^{-1} r^{(k)} \). Então

\[ \hat{s}_k = \frac{\langle \hat{r}^{(k)}, \hat{r}^{(k)} \rangle}{\langle \hat{r}^{(k-1)}, \hat{r}^{(k-1)} \rangle} = \frac{\langle C^{-1} r^{(k)}, C^{-1} r^{(k)} \rangle}{\langle C^{-1} r^{(k-1)}, C^{-1} r^{(k-1)} \rangle}, \]

de modo que

\[ s_k = \frac{\langle w^{(k)}, w^{(k)} \rangle}{\langle w^{(k-1)}, w^{(k-1)} \rangle}. \tag{7.31} \]

Assim,

\[ \begin{aligned} \hat{t}_k &= \frac{\langle \hat{r}^{(k-1)}, \hat{r}^{(k-1)} \rangle}{\langle \hat{v}^{(k)}, \hat{A} \hat{v}^{(k)} \rangle} = \frac{\langle C^{-1} r^{(k-1)}, C^{-1} r^{(k-1)} \rangle}{\langle C^t v^{(k)}, C^{-1} A C^{-t} C^t v^{(k)} \rangle} \\ &= \frac{\langle w^{(k-1)}, w^{(k-1)} \rangle}{\langle C^t v^{(k)}, C^{-1} A v^{(k)} \rangle}. \end{aligned} \]

e, já que

\[ \begin{aligned} \langle C^t v^{(k)}, C^{-1} A v^{(k)} \rangle &= [C^t v^{(k)}]^t C^{-1} A v^{(k)} \\ &= [v^{(k)}]^t C C^{-1} A v^{(k)} = [v^{(k)}]^t A v^{(k)} = \langle v^{(k)}, A v^{(k)} \rangle, \end{aligned} \]

temos

\[ t_k = \frac{\langle w^{(k-1)}, w^{(k-1)} \rangle}{\langle v^{(k)}, A v^{(k)} \rangle}. \tag{7.32} \]

Prosseguindo,

\[ \hat{x}^{(k)} = \hat{x}^{(k-1)} + \hat{t}_k \hat{v}^{(k)}, \quad \text{e assim,} \quad C^t x^{(k)} = C^t x^{(k-1)} + t_k C^t v^{(k)} \]

e

\[ x^{(k)} = x^{(k-1)} + t_k v^{(k)}. \tag{7.33} \]

Continuando,

\[ \hat{r}^{(k)} = \hat{r}^{(k-1)} - \hat{t}_k \hat{A} \hat{v}^{(k)}, \]

de modo que

\[ \begin{aligned} C^{-1} r^{(k)} &= C^{-1} r^{(k-1)} - t_k C^{-1} A C^{-t} \hat{v}^{(k)}, \\ r^{(k)} &= r^{(k-1)} - t_k A C^{-t} C^t v^{(k)}, \end{aligned} \]

e

\[ r^{(k)} = r^{(k-1)} - t_k A v^{(k)}. \tag{7.34} \]

Finalmente,

\[ \hat{v}^{(k+1)} = \hat{r}^{(k)} + \hat{s}_k \hat{v}^{(k)}, \quad \text{e} \quad C^t v^{(k+1)} = C^{-1} r^{(k)} + s_k C^t v^{(k)}, \]

e então,

\[ v^{(k+1)} = C^{-t} C^{-1} r^{(k)} + s_k v^{(k)} = C^{-t} w^{(k)} + s_k v^{(k)}. \tag{7.35} \]

O método do gradiente conjugado pré-condicionado baseia-se na utilização das Equações (7.31) a (7.35) na ordem (7.32), (7.33), (7.34), (7.31), (7.35). O Algoritmo 7.5 implementa esse procedimento.

O pré-condicionamento substitui um sistema dado por outro tendo as mesmas soluções, porém com melhores características de convergência.

O próximo exemplo ilustra os cálculos em um problema simples.