Polinômio dos quadrados mínimos condicionado

The constrained least squares polynomial

NOTAS E DISCUSSÃO

Polinômio dos quadrados mínimos condicionado

The constrained least squares polynomial

Oclide J. Dotto¹; Adalberto A. Dornelles Filho

Departamento de Matemática e Estatística, Universidade de Caxias do Sul, Caxias do Sul, RS, Brasil

RESUMO

O ajustamento polinomial a dados experimentais pelo método de quadrados mínimos é muito usado, e aqui estudamos a introdução de uma simples modificação, que, em certos casos, é mais adequada, isto é, condicionamos o polinômio a conter um dado ponto. Incluímos uma implementação do correspondente algoritmo em MATLAB.

Palavras-chave: ajustamento, quadrados mínimos, condicionamento.

ABSTRACT

Since the least-squares polynomial fit is used very often, we study here a simple modification that is more suitable for some cases, that is, we constrain the polynomial to pass through a given point, besides satisfying the least-squares property. We include an implementation of the corresponding algorithm in MATLAB.

Keywords: data fitting, least squares, constraint.

1. Introdução

É popular o ajustamento polinomial discreto dos quadrados mínimos (q. m.). Para evocar o processo desse ajustamento, consideremos os dados da Tabela 1, onde a primeira linha forma uma seqüência finita crescente.

Suponhamos que queiramos ajustar a esses dados o polinômio

de ordem k < n dos q. m. A melhor maneira de obter a solução é o uso da álgebra linear como segue. Começamos gerando a matriz de planejamento (design matrix [1])

de ordem (n+1)×(k+1), determinada pelo vetor x = [ x₀x₁ ... x_n ]^t. A solução única (porque X é não-singular) da equação normal

é o vetor a dos coeficientes do polinômio dos q. m. (1).

Do ponto de vista dos métodos de resolução numérica de (3), o uso da matriz de planejamento pode levar a imprecisões uma vez que X se torna mal-condicionada quando sua ordem é alta. No entanto, isso não é relevante na maioria dos casos práticos, onde dificilmente a ordem do polinômio ajustado é maior que 3.

O desenvolvimento mais completo da teoria e solução numérica dos quadrados mínimos foi feita por Åke Björck em [2], mas pode ser encontrado em muitos bons livros de álgebra linear, como [3, 4, 5]. Uma introdução rápida sobre a solução dos q. m. encontra-se em [6].

Façamos aqui uma pequena digressão útil. Na estimativa da solução Xa = y mediante o método dos q. m., é conveniente às vezes ponderar essa estimativa, isto é, usar e equação normal ponderada

em vez da equação normal ordinária (3), escolhendo, para matriz de ponderação W, uma conveniente matriz simétrica definida positiva, que minimize o vetor residual r : = Xa - y no sentido dos q. m. associados ao produto interno definido por W, o que significa simplesmente minimizar r^tWr. É conhecido [7] que a escolha ótima é a inversa da matriz V = [ v_ij ] das variâncias-covariâncias dos vetores aleatórios das observações, y_i = [ y_1iy_2i ... y_mi]^t, i = 0, 1, ..., n, isto é, W = V^-1 (cf. também [2, 8]). Ainda, caso os desvios possam ser supostos não-correlacionados, todas as covariâncias são nulas e V se reduz a uma matriz diagonal, cujos elementos diagonais são as variâncias, que medem a dispersão dos desvios em torno da média, e, nesse caso, para inverter V, basta inverter esses elementos diagonais (que são não-nulos). O efeito do uso de V^-1 para ponderação é terem menor peso os desvios maiores.

Voltemos agora ao assunto central do artigo. A curva que representa geometricamente o polinômio dos q. m. não passa pelos pontos definidos pelos dados na Tabela 1, a não ser que coincida com a do polinômio de interpolação, o que ocorre quando k = n. No entanto, em certos casos, pela natureza do problema, sabemos de antemão que deve passar por determinado ponto, mesmo que esse ponto não figure entre os dados.

A um polinômio dos q. m., forçado a passar por um ponto P₀ escolhido previamente, chamamos de polinômio dos q. m. condicionado, ou polinômio dos q. m. ancorado no ponto P₀.

2. Reta dos quadrados mínimos ancorada na origem

A lei de Hooke estabelece que, quando uma força F é aplicada a uma mola homogênea, o alongamento x da mola varia linearmente com essa força, dentro de limites, isto é,

onde a é dita a constante de Hooke da mola. Para determinar a, executamos experimentos que fornecem uma tabela de duas linhas com os valores de F dados e os valores de x resultantes. Em seguida ajustamos a reta dos q. m. aos dados e a declividade dessa reta é uma estimativa da constante de Hooke da mola. Essa reta quase certamente não passa pela origem, mesmo que a origem se inclua nos dados. No entanto, idealmente, deveria passar pela origem, pois o alongamento é nulo para uma força aplicada nula.

Para forçar a reta dos q. m. a passar pela origem, basta achar a solução dos q. m. do sistema linear

Uma relação de proporcionalidade direta entre duas grandezas, como em (5), ocorre com freqüência. Outros exemplos são: a intensidade i e a tensão v de uma corrente elétrica aplicada a um resistor, a força F e aceleração resultante a aplicada a um móvel, etc.

Exemplo 1

Vejamos um caso específico. A Tabela 2 mostra os valores do alongamento x (em centímetros) de uma mola vs. força de tração F (em Newtons) aplicada, ambos em módulo.

Aqui a reta dos q. m. tem por equação

e a reta dos q. m. ancorada na origem

A Fig. 1 mostra ambas as retas. A primeira não contém a origem. Vemos que as declividades, em unidades coerentes, não coincidem, embora sejam próximas, e é razoável dizer que o valor a = 0,76 é uma estimativa melhor que a = 0,78 para a constante da mola.

Para a reta dos q. m. ancorada na origem, a equação matricial (3) pode ser reescrita como

De (9) deduzimos que a declividade a dessa reta é dada por

Em laboratórios de ensino, em vez da fórmula (10), é usada comumente a fórmula

que expressa a média aritmética a das declividades y_i/x_i das retas, cada uma das quais contendo a origem e um ponto (x_i,y_i). A expressão (11) é bastante intuitiva e produz um resultado levemente diferente (a = 0,73, no caso) e, por outro lado, parece não ter nenhuma relação com (10). De fato, elas não são equivalentes. Mas é interessante observar o que segue.

Primeiro, adotando a equação normal ponderada (4) e a inversa V^-1 da matriz das variâncias-covariâncias para matriz de ponderação, a fórmula (9) passa a escrever-se

Segundo, suponhamos que, para cada i Î {0, 1, 2, ...,n }, tenhamos feito um conjunto de m observações y_1i, y_2i, ..., y_mi; então podemos falar da matriz das variâncias-covariâncias dos vetores aleatórios dessas observações, y_i = [ y_1i y_2i ... y_mi]^t, i = 0, 1, 2, ..., n. Usando a fórmula normal (12), ponderada pela inversa da matriz V = [ v_ij ] das variâncias-covariâncias, ou ponderando as declividades y_i/x_i em (11) pelos inversos das correspondentes variâncias dessas declividades (ao invés de tomar a média aritmética), obtemos por ambos os caminhos a fórmula

onde pusemos v_i : = v_ii (variâncias) e admitimos v_ij = 0, se i ¹ j (covariâncias nulas). Com relação a ponderação em (13), lembramos que, por uma propriedade da variância, temos:

var ,

e o denominador principal em (13) é a soma dos inversos dessas variâncias.

Notamos que a fórmula (10) é recuperada de (13) supondo v₀ = v₁ = ... = v_n, e, com uma simplificação maior, que a fórmula (11) segue de (13), usando as correspondentes variâncias do vetor das declividades

todas iguais a 1/(n+1), isto é, para i = 0, 1, 2 ..., n. Isso explica porque as fórmulas (10) e (11) não são equivalentes.

3. Polinômio dos quadrados mínimos condicionado

Em vez de ancorar na origem a reta dos q. m., podemos ancorá-la num ponto qualquer. Mas tratemos do caso geral: mostremos como ancorar o polinômio de ordem k dos q. m. num ponto escolhido (x₀, y₀). O procedimento é feito em três etapas:

1. operamos uma translação de eixos de maneira que a nova origem seja o ponto (x₀, y₀); para isso basta subtrair x₀ das abscissas e y₀ das ordenadas dos dados;

2. ajustamos aos dados transformados o polinômio q dos q. m. ancorado na origem, como no Exemplo 1;

3. fazemos, por último, a transformação inversa da executada na etapa 1, obtendo o polinômio final dos q. m., p(x) = q(x - x₀), ancorado em (x₀, y₀).

Suponhamos que os dados na Tabela 1 já tenham sofrido a alteração descrita na etapa 1. Para operacionalizar a etapa 2, suprimimos a última coluna da matriz de planejamento (2), obtendo uma matriz que indicaremos com U, e achamos a solução a dos q. m. do sistema linear

O vetor-solução a será formado pelos coeficientes do polinômio q ancorado na origem, descrito na etapa 2. Para o caso de uma reta, a equação matricial (15) torna-se uma equação escalar, da forma (9), onde x e y são os vetores dos dados.

No ^apêndice implementamos esse algoritmo em MATLAB, com o nome polqmcond. O programa produz também a correspondente visualização geométrica.

Exemplo 2

Consideremos a integral elíptica completa

Mediante integração numérica, obtemos os valores de E(x) (com x dado em graus) e construímos a Tabela 3.

Se olharmos para a distribuição desses dados, veremos que um ajustamento parabólico é adequado. Além disso, é óbvio que E(0) = p/2. Então é indicado ancorar a parábola dos q. m. no ponto (0, p/2). Utilizando o programa polqmcond, obtemos essa parábola

A parábola ''livre'' dos q. m. é dada por

A Fig. 2 exibe ambas as parábolas. Visualmente elas se superpõem em quase todo o trajeto.

Para compararmos os dois ajustamentos, os testamos nos pontos x = 2, 12, 17, 27. Obtivemos a Tabela 4 na qual a segunda linha contém os valores de E(x) obtidos por integração numérica aqui considerados exatos, a terceira linha os valores de P(x) do ajustamento dos q. m. ancorado, e a quarta linha os valores de Q(x) do ajustamento livre.

Vemos que, nesse exemplo, o ajustamento ancorado é melhor para os pontos próximos ao ponto ancorador, como era esperado, e longe dele os ajustamentos tendem a ser equivalentes. Claramente, para um ajustamento linear, por exemplo, a situação longe da âncora é outra, como mostra a Fig. 1.

4. Comentários finais

Nos laboratórios de ensino, é bastante comum realizar experimentos para obter o coeficiente de proporcionalidade entre duas grandezas diretamente proporcionais. Se for usado o recurso do ajuste linear dos q. m., disponível em algumas calculadoras científicas, dificilmente obteremos uma reta que passe exatamente pela origem. Uma elegante solução para essa dificuldade é usar a fórmula (10) que decorre do ajustamento dos q. m. ancorado na origem.

Mais geral, mostramos que é possível e relativamente fácil, ao ajustarmos um polinômio dos quadrados mínimos, condicioná-lo a conter um ponto de interesse. Poder-se-ia tentar desenvolver um algoritmo que o condicione a conter mais de um ponto, contanto que ainda retenha um número de graus de liberdade suficientemente elevado. Por exemplo, um caso extremo oposto é condicionar uma reta dos q. m. a conter dois pontos, o que corresponde a ela ficar com zero graus de liberdade e equivale ao problema trivial de achar a equação da reta que passa por esses dois pontos. Não temos certeza se há interesse prático nessa generalização de ancoragem. De qualquer maneira, parece-nos que o caso interessante, exeqüível e útil, é a ancoragem em um ponto apenas.

Agradecimento

Os autores agradecem ao parecerista a gentil apreciação do artigo, assim como as sugestões.

Recebido em 20/2/2006; Revisado em 15/5/2006; Aceito em 29/5/2006

• [1] Richard A. Johnson and Dean W. Wichern, Applied Multivariate Analysis (Prentice Hall, Nova York, 2002)
• [2] Ĺke Björck, Numerical Methods for Least Squares Problems (SIAM, Philadelphia, 1966).
• [3] James W. Demmel, Applied Numerical Linear Algebra (SIAM, Philadelphia, 1997).
• [4] Gene H. Golub and Charles F. Van Loan, Matrix Computations (The John Hopkins University Press, London, 1996).
• [5] Lloyd N. Trefethen and David Bau III, Numerical Linear Algebra (SIAM, Philadelphia, 1997).
• [6] Mário Barone Júnior, Álgebra Linear (IME-USP, Săo Paulo, 2005), 3Ş ed.
• [7] Gilbert Strang and Kay Borre, Linear Algebra, Geodesy and GPS (Wellesley-Cambridge Press, Wellesley, 1997).
• [8] José Henrique Vuolo, Fundamentos da Teoria de Erros (Ed. Edgard Blücher Ltda., Săo Paulo, 2005), 2Ş ed.

apêndice

Datas de Publicação

• Publicação nesta coleção
  06 Nov 2006
• Data do Fascículo
  2006

Histórico

• Aceito
  29 Maio 2006
• Revisado
  15 Maio 2006
• Recebido
  20 Fev 2006