Quadrados mínimos não lineares

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 é $\phi:\R\times\R^n\to\R$, $(x;c)\mapsto \phi(x;c),$ onde $x\in\R$ é a variável livre e $c\in\R^n$ são parâmetros (coeficientes) da função $\phi$, a serem determinados, de modo que o resíduo seja o menor possível. Por exemplo, queremos determinar $c = (c_1,c_2)^T$ tais que $$\phi(x_k;c) \equiv c_1\exp(c_2x_k) \approx y_k.$$ Mais precisamente, o vetor $c$ ótimo é definido como solução do problema \begin{equation}\label{qm09:qm1} \min_c \sum_{k=1}^m [\phi(x_k;c) - y_k]^2. \end{equation} Para melhorar a notação, defina $F:\R^n\to\R^m$ como \begin{equation} F(c) \equiv \left[\begin{array}{c} \phi(x_1;c) - y_1\\\ \phi(x_2;c) - y_2\\\ \vdots \\\ \phi(x_m;c) - y_m \end{array}\right]. \end{equation} Desta forma, o problema \eqref{qm09:qm1} torna-se determinar $c$, solução de \begin{equation}\label{qm09:qm3} \min_c\|F(c)\|_2^2. \end{equation} Este é um problema de ajuste por quadrado míninos não linear.

No caso em que $F:\R^n\to\R^m$, com $n \lt m$, não tenha vindo especificamente de problema de ajuste, dizemos que o vetor $c\in\R^n$, solução de \eqref{qm09:qm3}, é a solução de quadrados mínimos do sistema não linear $F(c) = 0.$

Linearizações sucessivas

Como já sabemos determinar a solução de quadrados mínimos para um sistema linear, vamos tentar recair nesse problema, aproximando $F(c)$ por um modelo linear, ao redor de um vetor $c$ de referência. Para isso, seja $J(c)\in\R^{m\times n}$ a matriz Jacobiana de $F$ em $c$, dada por $$ [J(c)]_{ij} = \frac{\partial f_i(c)}{\partial c_j}, $$ onde $f_i$ é a $i$-ésima componente de $F.$ Se $J$ for diferenciável , então é possível mostrar que $$ F(c+s) = F(c) + J(c)s + {\cal O}(\|s\|^2) \equiv M(c) + {\cal O}(\|s\|^2). $$ Dizemos que $M(c) \equiv F(c) + J(c)s$ é a aproximação linear de $F(c+s)$. Buscar a solução de quadrados mínimos para o modelo linear é uma forma de aproximar a solução de quadrados mínimos para o problema não linear. Ao passar do vetor $c$ para o vetor $c+s$, devemos estar nos aproximando da solução de quadrados mínimos desejada. Caso ainda não estejamos satisfeitos, o processo pode ser repetido. Este processo de sucessivas linearizações dá origem ao método de Gauss-Newton, sumarizado no algoritmo a seguir.

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.

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 é $\phi(x;c) = c_1+\cos(c_2x),$ defina os vetores com os pontos amostrados e a curva, com

# 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$ é um vetor que computa a diferença $\phi(x_i;c)-y_i$ para $i=1,\ldots,5$. Fazendo uso da notação vetorial do Octave, defina $F$ com:

F = @(c) phi(x,c) - y;

A seguir, defina a matriz Jacobina, que tem como primeira coluna as derivadas de $F$ com respeito a $c_1$ e como segunda colunas as derivadas de $F$ com respeito a $c_2$.

# 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 $\phi(x;c) = c_1\exp(c_2x)$. Suponha que os pontos amostrados são $(x_1,y_1),\ldots,(x_m,y_m).$ Defina $F(c)$. Exiba a matriz Jacobiana $J(c)$.

2. Determine $c_1$ e $c_2$ tais que $\phi(x) = e^{c_1x} - 2 e^{c_2x}$ se ajuste aos dados da tabela.