Algoritmo de Gauss-Newton

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

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

$S(\mathit{\beta })=\sum _{i=1}^m{r}_i^2(\mathit{\beta }).$

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

${\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{(s)}-{\left({{𝐉}_{𝐫}}^{\mathrm{\top }}{𝐉}_{𝐫}\right)}^{-1}{{𝐉}_{𝐫}}^{\mathrm{\top }}𝐫({\mathit{\beta }}^{(s)})$

onde

${𝐉}_{𝐫}=\frac{\partial r_i}{\partial {\beta }_j}({\mathit{\beta }}^{(s)})$

é a Matriz Jacobiana de "r" e o símbolo ${\mathrm{\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 (x_i, y_i), as funções r_i são os resíduos

$r_i(\mathit{\beta })=y_i-f(x_i,\mathit{\beta }).$

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

${\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{(s)}-{\left({{𝐉}_{𝐟}}^{\mathrm{\top }}{𝐉}_{𝐟}\right)}^{-1}{{𝐉}_{𝐟}}^{\mathrm{\top }}𝐫({\mathit{\beta }}^{(s)}).$

Notas

A suposição m ≥ n na demonstração do algoritmo é necessária, senão a matriz J_r^TJ_r 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 r_i. Usando o Teorema de Taylor, podemos escrever em cada iteração:

$𝐫(\mathit{\beta })\approx 𝐫({\mathit{\beta }}^s)+{𝐉}_{𝐫}({\mathit{\beta }}^s)\mathrm{\Delta }$

com $\mathrm{\Delta }=\mathit{\beta }-{\mathit{\beta }}^s.$

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

$𝐦𝐢𝐧\Vert 𝐫({\mathit{\beta }}^s)+{𝐉}_{𝐫}({\mathit{\beta }}^s)\mathrm{\Delta }{\Vert }_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 J_r. 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 J_r, as iterações falharão, já que J_r^TJ_r se tornará singular.

Exemplo

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.$ Seja ${\beta }_1=V_{\text{max}}$ e ${\beta }_2=K_M.$ Encontraremos ${\beta }_1$ e ${\beta }_2$ tal que a soma dos quadrados dos resíduos

$r_i=y_i-\frac{{\beta }_1{x}_i}{{\beta }_2+x_i}$ ($i=1,\dots ,7$)

será minimizada.

A Jacobiana ${𝐉}_{𝐫}$ 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 }_1{x}_i}{{({\beta }_2+x_i)}^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. Podemos também determinar as erros: ${\hat{\beta }}_1=0,36$ ± $0,07$ e ${\hat{\beta }}_2=0,56$ ± $0,35$ com 80% de confiança. A soma dos quadrados dos resíduos 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

Pode ser mostrado 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. O algoritmo pode convergir lentamente ou nunca se a suposição inicial estiver longe do mínimo ou se a matriz $J_{𝐫}^{𝐓}{𝐉}_{𝐫}$ estiver mal condicionada. Por exemplo, considere o problema com m=2 equações e n=1 variáveis, dadas por:

$\begin{array}{cc}{r}_1(\beta )& =\beta +1\\ r_2(\beta )& =\lambda {\beta }^2+\beta -1.\end{array}$

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.

Derivação com o Método de Newton

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 β, é

${\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{(s)}-{𝐇}^{-1}𝐠$

onde g denota o vetor gradiente de S e H denota a 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(\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}).$

O método de Gauss-Newton é obtido ignorando os termos derivados de segunda ordem (o segundo termo nesta expressão). Isto é, o Hessian é aproximado 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 J_r. O gradiente e o Hessien aproximado podem ser escritos numa notação matricial como:

$𝐠=2{{𝐉}_{𝐫}}^{\mathrm{\top }}𝐫,𝐇\approx 2{{𝐉}_{𝐫}}^{\mathrm{\top }}{𝐉}_{𝐫}.$

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

${\mathit{\beta }}^{(s+1)}={\mathit{\beta }}^{(s)}+\mathrm{\Delta };\mathrm{\Delta }=-{\left({{𝐉}_{𝐫}}^{\mathrm{\top }}{𝐉}_{𝐫}\right)}^{-1}{{𝐉}_{𝐫}}^{\mathrm{\top }}𝐫.$

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.

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

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({\mathit{\beta }}^s)$ seja um ponto estacionário, mantem-se que $S({\mathit{\beta }}^s+\alpha \mathrm{\Delta })<S({\mathit{\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.

${\mathit{\beta }}^{s+1}={\mathit{\beta }}^s+\alpha \mathrm{\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 acentuada.

$\left({𝐉}^{𝐓}𝐉\mathbf{+}\mathit{\lambda }𝐃\right)\mathrm{\Delta }={𝐉}^T𝐫$,

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

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 toda vez que $\lambda$ for alterado. Uma estratégia mais eficiente é a seguinte: quando a divergência ocorre, deve-se 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âmetro 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.

Fontes

• Predefinição:Citar livro
• Predefinição:Citar livro.
• Predefinição:Citar livro

Predefinição:Isaac Newton