Algoritmo de Gauss-Newton

Origem: Wikipédia, a enciclopédia livre.
Ir para: navegação, pesquisa

O algoritmo de Gauss-Newton é um método usado para resolver problemas de mínimos quadrados não lineares. Ele pode ser visto como uma modificação do Método de Newton para achar o mínimo de uma função. Diferentemente do Método de Newton, o Algoritmo de Gauss-Newton apenas pode ser usado para minimizar uma soma dos valores quadrados da função, mas tem a vantagem de que as derivadas segundas, que podem ser difíceis de calcular, não são necessárias.

Problemas de mínimos quadrados não lineares surgem, por exemplo, em regressão não linear, onde os parâmetros de um modelo são procurados de forma que o modelo esteja em concordância com as observações disponíveis.

O método foi nomeado a partir dos matemáticos Carl Friedrich Gauss e Isaac Newton.

Descrição[editar | editar código-fonte]

Dada "m" funções r = (r1, …, rm) de n variáveis 'β = (β1, …, βn), com m  ≥ n', o Algoritmo de Gauss-Newton iterativamente encontra o mínimo das somas dos quadrados[1]

 S(\boldsymbol \beta)= \sum_{i=1}^m r_i^2(\boldsymbol \beta).

Começando com uma estimativa inicial \boldsymbol \beta^{(0)} para o mínimo, o método prossegue com as iterações

 \boldsymbol \beta^{(s+1)} = \boldsymbol \beta^{(s)} - \left(\mathbf{J_r}^\top \mathbf{J_r} \right)^{-1} \mathbf{ J_r} ^\top \mathbf{r}(\boldsymbol \beta^{(s)})

onde

 \mathbf{J_r} = \frac{\partial r_i}{\partial \beta_j}(\boldsymbol \beta^{(s)})

é a Matriz Jacobiana de "r" e o símbolo ^\top denota a matriz transposta.

Na montagem de dados, onde o objetivo é encontrar os parâmetros β tais que uma dada função modelo y = f(x, β) ajuste melhor alguns pontos de dados (xi, yi), as funções ri são os resíduos

r_i(\boldsymbol \beta)= y_i - f(x_i, \boldsymbol \beta).

Então, o método de Gauss-Newton pode ser expresso em termos da Jacobiana da função "f" como

 \boldsymbol \beta^{(s+1)} = \boldsymbol \beta^{(s)} + \left(\mathbf{J_f}^\top \mathbf{J_f} \right)^{-1} \mathbf{ J_f} ^\top \mathbf{r}(\boldsymbol \beta^{(s)}).

Notas[editar | editar código-fonte]

A suposição m ≥ n na demostração do algoritmo é necessária, senão a matriz JrTJr não será invertível e as equações normais não poderão ser resolvidas (pelo menos exclusivamente).

O Algoritmo de Gauss-Newton pode ser obtido por aproximação linear do vetor das funções ri. Usando o Teorema de Taylor, podemos escrever em cada iteração:

\mathbf{r}(\boldsymbol \beta)\approx \mathbf{r}(\boldsymbol \beta^s)+\mathbf{J_r}(\boldsymbol \beta^s)\Delta

com \Delta=\boldsymbol \beta-\boldsymbol \beta^s.

A tarefa de encontrar Δ minimizando a soma dos quadrados do lado direito, por exemplo:

\mathbf{min}\|\mathbf{r}(\boldsymbol \beta^s)+\mathbf{J_r}(\boldsymbol \beta^s)\Delta\|_2^2

é um problema linear de mínimos quadrados, que pode ser resolvido explicitamente, obtendo-se as equações normais no algoritmo.

As equações normais são m equações normais simultâneas no desconhecido incremento Δ. Elas podem ser solucionadas em um passo, usando decomposição de Cholesky, ou, melhor, a fatoração QR de Jr. Para sistemas de grandes dimensões, um método iterativo, tal como o método do gradiente conjugado, pode ser mais eficiente. Se existe uma dependência linear entre as colunas de Jr, as iterações falharão, já que JrTJr se tornará singular.

Exemplo[editar | editar código-fonte]

Curva calculada obtida com \hat\beta_1=0.362 e \hat\beta_2=0.556 (em azul) e os dados observados (em vermelho).

Neste exemplo, o Algoritmo de Gauss-Newton será usado para ajustar um modelo a alguns dados, minimizando os erros das somas dos quadrados entre os dados e as previsões do modelo.

Numa experiência de biologia, estudando a relação entre a concentração de substrato ["S"] e a taxa de reação de uma reação mediada por enzimas, foram obtidos os dados da tabela a seguir:

i 1 2 3 4 5 6 7
[S] 0.038 0.194 0.425 0.626 1.253 2.500 3.740
taxa 0.050 0.127 0.094 0.2122 0.2729 0.2665 0.3317

É desejado encontrar uma curva (função modelo) com a fórmula

\text{taxa}=\frac{V_\text{max}[S]}{K_M+[S]}

que melhor se adapte aos dados, no sentido dos mínimos quadrados, com os parâmetros V_\text{max} e K_M a serem determinados.

Denotado por x_i e y_i o valor de [S] e a taxa da tabela i=1, \dots, 7. Permite que \beta_1=V_\text{max} e \beta_2=K_M. Acharemos \beta_1 e \beta_2 tal que a soma dos quadrados dos resíduos

r_i = y_i - \frac{\beta_1x_i}{\beta_2+x_i}   (i=1,\dots, 7)

será minimizada.

A Jacobiana \mathbf{J_r} do vetor dos resíduos r_i em relação às incógnitas \beta_j é uma matriz 7\times 2 com i-th linhas de entrada

\frac{\partial r_i}{\partial \beta_1}= -\frac{x_i}{\beta_2+x_i},\  \frac{\partial r_i}{\partial \beta_2}= \frac{\beta_1x_i}{\left(\beta_2+x_i\right)^2}.

Começando com as estimativas iniciais de \beta_1=0.9 e \beta_2=0.2, depois de cinco iterações do Algoritmo de Gauss-Newton os valores otimizados \hat\beta_1=0.362 e \hat\beta_2=0.556 são obtidos. A soma dos quadrados dos residuos diminuiu desde o valor inicial de 1.445 para 0.00784 depois da quinta iteração. O gráfico à direita mostra a curva determinada pelo modelo dos parâmetros otimizados e os dados observados.

Propriedades de Convergência[editar | editar código-fonte]

Pode ser mostrado[2] que o incremento Δ é um sentido descendente para "S", e, se o algoritmo converge, então o limite é um ponto estacionário de "S". Contudo, a convergência não é garantida, nem mesmo a convergência local como no Método de Newton.

A taxa de convergência do Algoritmo de Gauss-Newton pode se aproximar do quadrático.[3] O algoritmo pode convergir lentamente ou nunca se a suposição inicial estiver longe do mínimo ou se a matriz \mathbf{J_r^T  J_r} estiver mal condicionada. Por exemplo, considere o problema com m=2 equações e n=1 variáveis, dadas por:

 \begin{align}
r_1(\beta) &= \beta + 1 \\
r_2(\beta) &= \lambda \beta^2 + \beta - 1.
\end{align}

A otimização é em \beta = 0. Se |λ| = 0, então o problema é de fato linear e o método converge em uma iteração. Se |λ| < 1, então o método converge linearmente e o erro decresce assintóticamente com um fator |λ| em cada iteração. No entanto, se |λ| > 1, então o método não converge nem mesmo localmente.[4]

Derivação com o Método de Newton[editar | editar código-fonte]

No que segue, o Algoritmo de Gauss-Newton será derivado com o Método de Newton por otimização da função através de uma aproximação. Como consequência, a taxa de convergência do Algoritmo de Gauss-Newton será no máximo quadrática.

A relação de recorrência do Método de Newton para minimizar uma função "S"de parâmetros β, é

 \boldsymbol\beta^{(s+1)} = \boldsymbol\beta^{(s)} - \mathbf H^{-1} \mathbf g \,

onde g denota o gradient|vetor gradiente de S e H denota a Hessian matrix|Matriz de Hessian de S. Uma vez que  S = \sum_{i=1}^m r_i^2, o gradiente é dado por:

g_j=2\sum_{i=1}^m r_i\frac{\partial r_i}{\partial \beta_j}.

Elementos de Hessien são calculados através da diferenciação dos elementos do gradiente, g_j, com respeito a \beta_k

H_{jk}=2\sum_{i=1}^m \left(\frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}+r_i\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k} \right).

O método de Gauss-Newton é obtido ignorando os termos derivados de segunda ordem (o segundo termo nesta expressão). Isto é, o Hessian é aproxiamado por:

H_{jk}\approx 2\sum_{i=1}^m J_{ij}J_{ik}

onde J_{ij}=\frac{\partial r_i}{\partial \beta_j} são entradas da Jacobiana Jr. O gradiente e o Hessien aproximado podem ser escritos numa notação matricial como:

\mathbf g=2\mathbf{J_r}^\top \mathbf{r}, \quad \mathbf{H} \approx 2 \mathbf{J_r}^\top \mathbf{J_r}.\,

Estas expressões são substituídas na relação de recorrência acima pra obter as equações operacionais.

 \boldsymbol{\beta}^{(s+1)} = \boldsymbol\beta^{(s)}+\Delta;\quad \Delta = -\left( \mathbf{J_r}^\top \mathbf{J_r} \right)^{-1} \mathbf{J_r}^\top \mathbf{r}.

A convergência do método de Gauss-Newton não é garantida em todas as instâncias. A aproximação

\left|r_i\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k}\right| \ll \left|\frac{\partial r_i}{\partial \beta_j}\frac{\partial r_i}{\partial \beta_k}\right|

que precisa ser capaz de ignorar os termos derivados de segunda ordem, pode ser válida em dois casos, nos quais a convergência é de se esperar.[5]

  1. Os valores da função r_i são pequenos em magnitude, ao menos em torno do mínimo.
  2. As funções são apenas "ligeiramente"não lineares, de modo que \frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k} seja relativamente pequena em magnitude.

Versões Aperfeiçoadas[editar | editar código-fonte]

Com o método de Gauss-Newton a soma dos quadrados "S"pode não decrescer em cada iteração. Contudo, uma vez que Δ é um sentido descendente, a menos que S(\boldsymbol \beta^s) seja um ponto estacionário, afirma que S(\boldsymbol \beta^s+\alpha\Delta) < S(\boldsymbol \beta^s) para todos suficientemente pequenos \alpha>0. Assim, se ocorre divergência, uma solução é a de empregar uma fração \alpha, do vetor incremento Δ na atual fórmula.

 \boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\alpha\  \Delta.

Em outras palavras, o vetor incremento é muito longo, mas ele aponta "para baixo", então somente uma parte irá decrescer o objetivo da função "S". Um valor ideal para \alpha pode ser encontrado usando um algoritmo de linha de pesquisa, ou seja, a magnitude de \alpha é determinada encontrando o valor que minimiza "S", geralmente usando um método de pesquisa direto no intervalo 0<\alpha<1.

Nos casos em que a direção do vetor de deslocamento é tal que a fração otimizada,  \alpha , é próxima de zero, um método alternativo para o tratamento de divergência é a utilização do algoritmo de Levenberg-Marquardt, também conhecido como o "método da região de segurança". As equações normais são modificadas de tal forma que o vetor de incremento seja virado na direção de descida mais ecentuada.

\left(\mathbf{J^TJ+\lambda D}\right)\Delta=\mathbf{J}^T \mathbf{r},

onde "D"é uma matriz diagonal positiva. Note que quando "D"é a matriz identidade e \lambda\to+\infty, então  \Delta/\lambda\to \mathbf{J}^T \mathbf{r}, portanto a direção de Δ se aproxima da direção do gradiente  \mathbf{J}^T \mathbf{r}.

O chamado parâmetro de Marquardt, \lambda, também pode ser otimizado por uma linha de pesquisa, mas é ineficiente já que o vetor de deslocamento deve ser recalculado tpda vez que \lambda for aletreado. Uma estratégia mais eficiente é a seguinte: quando a divergência ocorre, se deve aumentar o parâmetro de Marquartdt até que haja uma diminuição em "S". Então deve-se manter o valor de uma iteração para a outra, mas, se possível, diminuí-lo até que um valor de corte seja atingido quando o parâmentro de Marquardt pode ser definido em zero; a minimização de "S" se torna então um padrão de minimização de Gauss-Newton.

Notas[editar | editar código-fonte]

  1. Björck (1996)
  2. Björck (1996) p260
  3. Björck (1996) p341, 342
  4. Fletcher (1987) p.113
  5. Nocedal (1997)[falta página]

Referências[editar | editar código-fonte]