Quadrados mínimos não lineares
Como buscar a solução de quadrados mínimos de problemas não lineares
Vamos lá?
Já vimos como determinar a solução de quadrados mínimos para um sistema linear. Mas e se o sistema for não linear? E de onde podem surgir sistemas não lineares?
- 1
- 2
O problema de quadrados mínimos e o problema de ajuste de curva são o mesmo problema?
Vimos que o problema de ajuste de curvas recai em um sistema linear sobredeterminado, quando a curva buscada é combinação linear de funções conhecidas. Como dificilmente um sistema sobredeterminado tem solução, nos contentamos em buscar a solução de quadrados mínimos, quer seja pela resolução do sistema normal, quer seja através da decomposição QR.
Considere o caso mais geral, em que que a função de ajuste é
No caso em que
Linearizações sucessivas
Como já sabemos determinar a solução de quadrados mínimos para um sistema linear, vamos tentar recair nesse problema, aproximando
Sejam dados
- Determine
que minimiza - Se
, retorne e FIM. - Se
, retorne ao primeiro passo. Caso contrário, encerre sem convergência.
Exemplo no Octave
No vídeo, vimos um exemplo de aplicação do método de Gauss-Newton a um problema de ajuste de curva. Vejamos como reproduzir esse exemplo com o Octave.
Inicie implementando o algoritmo acima em uma função para ser chamada no Octave.
gaussnewton.m
function [c, k] = gaussnewton(F,J,c0,epsilon,K) ## [c, k] = gaussnewton(F,J,c0,epsilon,K) ## ## Implementação simples do método de Gauss-Newton. ## ## F: R^n -> R^m ## J: R^n -> R^{m x n}, Jacobiana de F ## c0, aproximação inicial ## epsilon, tolerância para o critério de parada ## K, número máximo de iterações k = 0; c = c0; do ## Solução de quadrados mínimos para F(c)+J(c)s=0 s = - J(c)\F(c); if (k==0) norms0 = norm(s,inf); norms = norms0; else norms = norm(s,inf); endif c = c + s; if (norms < epsilon * norms0) return endif k = k +1; until (k == K) endfunction
Para usar esta função precisamos definir os dados particulares do problema. Lembrando que curva desejada é
# Dados amostrados x = [0.0 0.5 1.00 1.5 2.0]'; y = [2.5 2.0 0.99 0.5 1.0]'; # Curva phi phi = @(x,c) c(1) + cos(c(2)*x);
A função
F = @(c) phi(x,c) - y;
A seguir, defina a matriz Jacobina, que tem como primeira coluna as derivadas de
# Derivadas parciais dFdc1 = @(c) ones(size(x)); # vetor de 1's, do tamanho do vetor x dFdc2 = @(c) -x.*sin(c(2)*x); # Matriz Jacobina J = @(c) [dFdc1(c) dFdc2(c)];
Já temos os ingredientes que definem o problema. Resta agora escolher a aproximação inicial, a tolerância para o critério de parada, o limite máximo para a quantidade de iterações, e chamar a função gaussnewton
.
c0 = [1;1]; epsilon = 1e-3; K = 20; [c,k] = gaussnewton(F,J,c0,epsilon,K) c = 1.4978 2.0968 k = 6 # Gráfico dessa primeira aproximação plot(x,y,'or',X,phi(X,c))
Referência
Åke Björck. Numerical Methods for Least Squares Problems. SIAM, 1996.
J. E. Dennis e R. B. Schnabel. Numerical Methods for Unconstrained Optimization and Nonlinear Equations. SIAM, 1996.
1. Seja
2. Determine