Quadrados_Mínimos

Prévia do material em texto

Quadrados Mínimos
Paulo J. S. Silva
12 de setembro de 2017
1 Introdução
Um dos temas principais de problemas numéricos é a obtenção de funções, ou leis, que ajustem dados expe-
rimentais. Por exemplo, imagine que temos dados experimentais de um fenômeno que pode ser explicado
por um polinômio do segundo grau. Mesmo que o fenômeno seja perfeitamente explicável por esse tipo de
função os erros de medidas nunca vão nos permitir observar exatamente o gráfico da verdadeira função.
Vamos simular isso no código abaixo.
In [ ]: using PyPlot
# Uma função do segundo grau.
lei(x) = *x.^2 + 2.0.*x + 1.73
# Vamos simular medidas reais de valores da função introduzindo um erro aleatório
# em [-5, 5].
n_pontos = 500
x = linspace(-4, 4, n_pontos)
medida = lei(x) + 10*(rand(n_pontos) - 0.5)
scatter(x, medida)
Ao observar a figura acima podemos ter a intruição que o valor medido vem de fato de um polinômio
do segundo grau. Mas como descobrir qual é esse polinômio. Ou seja, como recuperar os coeficientes
a2 = π, a1 = 2, a0 = 1.73? Este é o assunto dessa parte do curso.
2 Modelando o problema de ajuste polinômios de grau 2
Vamos tentar achar uma formulação matemática razoável para o problema e depois resolvê-lo. O que temos
são medidas com erros de função p cujo o formato geral conhecemos
p(x) = a2x2 + a1x + a0,
para coeficientes a2, a1, a0 que não conhecemos. Desejamos descobrir os valores dos seus coeficientes a2, a1
e a0. A informação que temos disponível são medidas (xi, yi), i = 1, . . . , m e que yi = p(xi) mais algum
tipo de ruído aleatório. Note que se não houvesse o ruído aleatório, isso é se tivéssemos o valor exato
do polinômio desejado a solução seria simples. Bastaria pegar três pontos do plano e encontrar o (único)
polinômio do segundo grau que passa por eles. Veremos como fazer esse isso na próxima parte do curso.
Note que quaisquer que fosse os três pontos escolhidos iríamos obter o mesmo polinômio.
No caso de erros de medidas na função fica impossível usar a ideia acima já que o polinômio iria mudar
um pouco dependendo de quais três pontos fossem escolhidos. Nesse caso, uma ideia razoável nesse caso
é procurar qual o polinômio do segundo grau que melhor se "encaixa"a todos os dados ao mesmo tempo.
Como fazer isso matematicamente?
Dado um i = 1, . . . , m, o erro de uma função p com respeito ao valor medido em um único ponto xi é
|p(xi)− yi|.
1
Para levar em consideração todos os pontos, os erros em todos eles devem ser combinados de alguma
forma. Uma maneira natural é pensar em dois vetores do Rm. O primeiro representa os valores medidos
y = (y1, y2, . . . , ym) e o segundo a função p nos respectivos pontos p = (p(x1), p(x2), . . . , p(xm)). O que
gostaríamos é que y = p mas devido aos erros isso é impossível. Podemos então buscar os coeficientes
a2, a1 e a0 que minimizem a distância entre y e p, já que isso leva em conta todas as coordenadas. Ou seja
queremos encontrar os coeficientes para os quais
∥p − y∥2
é o menor possível. Uma coisa interessante é que minimizar a norma euclidiana ou o seu quadrado dá a
mesma resposta. Já a norma euclidiana ao quadrado tem a vantagem de ser diferenciável o que poderemos
explorar ao resolver o problema.
Chegamos finalmente ao seguinte problema:
Encontre coeficientes a2, a1, a0 de modo que ao definirmos a função p(x) = a2x2 + a1x + a0 o valor de
m
∑
i=1
(p(xi)− yi)
2 = ∥p − y∥2 (1)
seja o menor possível.
3 Ajuste linear de curvas
Uma observação importante é que o fato da função que desejamos ajustar aos dados medidos ser um po-
linômio de grau 2 desempenhou um papel secundário na dedução do problema final acima. A mesma
ideia poderia ser aplicada para ajustar um polinômio de grau 5, por exemplo. Mais ainda se função que
desejamos ajustar é uma soma de senos e cossemos como
c0 sin(x) + c1 cos(x) + c2 sin(2x) + c3 cos(2x)
o mesmo problema faz sentido e a modelagem é basicamente a mesma.
De uma maneira geral, considere que temos n funções ϕ1, ϕ2, . . . , ϕn que acreditamos que se se combi-
nam de forma linear para definir uma função ϕ. Isto é, acreditamos que ϕ tem a forma
ϕ(x) = c1ϕ1(x) + c2ϕ2(x) + . . . + cnϕn(x)
usando um vetor de coeficientes c = (c1, c2, . . . , cn) ∈ Rn. Considere também que temos um conjunto de
medidas (xi, yi), i = 1, . . . , m que representam, aproximadamente, avaliações da função ϕ com o formato
acima.
Podemos então desejar encontrar os coeficientes que definem a função ϕ que melhor se ajusta aos dados
no sentido de minimizar
∥ϕ − y∥2,
em que y é o vetor de medidas y = (y1, y2, . . . , ym) e ϕ o vetor da função ϕ(·) avaliado nos pontos xi, i =
1, 2, . . . , m:
ϕ =

ϕ(x1)
ϕ(x2)
...
ϕ(xm)
 =

c1ϕ1(x1) + c2ϕ2(x1) + . . . + cnϕn(x1)
c1ϕ1(x2) + c2ϕ2(x2) + . . . + cnϕn(x2)
...
c1ϕ1(xm) + c2ϕ2(xm) + . . . + cnϕn(xm)
 .
Escolhemos o critério de minimizar a norma euclidiana elevada ao quadrado por que assim a função que
desejamos minimizar é diferenciável.
Obtemos assim o problema de encontrar um vetor de coeficientes c = (c1, c2, . . . , cn) tal que ao definir-
mos a função ϕ(x) = c1ϕ1(x) + c2ϕ2(x) + . . . + cnϕn(x) o valor de
m
∑
i=1
(ϕ(xi)− yi)
2 = ∥ϕ − y∥2 (2)
2
seja o menor possível.
Para sermos capaz de realizar essa tarefa temos que entender como a expressão acima varia com os
coeficientes no vetor c. Para isso vamos usar um pouco de notação matricial. Se olharmos a segunda
fórmula acima e tentarmos entender o que ela diz colocando c em destaque temos:
ϕ =

ϕ(x1)
ϕ(x2)
...
ϕ(xm)

︸ ︷︷ ︸
ϕ
=

ϕ1(x1) ϕ2(x1) . . . ϕn(x1)
ϕ1(x2) ϕ2(x2) . . . ϕn(x2)
...
...
...
...
ϕ1(xm) ϕ2(xm) . . . ϕn(xm)

︸ ︷︷ ︸
Φ

c1
c2
...
cn

︸ ︷︷ ︸
c
. (3)
Como sugerido acima, podemos definir a matriz Φ como a matriz acima, ou seja a matriz que contém
na primeira coluna a função ϕ1 avaliada nos pontos x1, x2, . . . , xm, depois na segunda coluna a função ϕ2
avaliada nos mesmos pontos e assim sucessivamente. Com isso podemos escrever essa última expressão
na forma compacta
ϕ = Φc.
Assim, a função que queremos minimizar pode ser escrita tornando explícita sua dependência no vertor de
coeficientes como
f (c) = ∥Φc − y∥2 = (Φc − y)t(Φc − y) = ctΦtΦc − 2ctΦty − yty.
Como desejamos minimizar a função f : Rn → R podemos usar os nossos conhecimentos de cálculo.
Em primeiro lugar, observe que a função cresce de forma quadrática com c, logo para c grande ela passará
a crescer até infinito se a norma de c crescer indefinidamente. Logo ela atinge mínimo. Agora, sabemos de
cálculo que os candidatos a ponto de mínimo devem obedecer à esquação
∇ f (c) = 0.
Precismos descobrir como calcular a derivada de f . Para isso observe as seguintes regras de derivação que
são facilmente demonstradas a partir das respectivas fórmulas.
1. Se g(c) = ct Ac, para alguma matriz quadrada A, então ∇g(c) = 2Ac.
2. Se h(c) = ctb = btc, para algum vetor b, então ∇h(c) = b, um vetor constante.
3.1 Equações normais
Aplicando as regras de cálculo acima para obter o gradiente de f :
∇ f (c) = 2ΦtΦc − 2Φty.
Assim a equação que caracteriza o mínimo é
0 = ∇ f (c) = 2ΦtΦc − 2Φty.
Ou escrevendo de uma forma mais simples, eliminando a constante 2,
ΦtΦc = Φty.
Um sistema de equações lineares em c com n equações e n incógnitas que aprendemos a resolver ante-
riormente. Esse tipo de sistema é conhecido como equação normal associada ao problema de quadrados
mínimos. Em particular se a matriz do sistema ΦtΦ for inversível sabemos que haverá uma única solução
que será necessariamente então o mínimo da função f .
Nos problemas que vamos enfrentar é natural imaginar que a matriz ΦtΦ é inversível. Isso somente não
ocorrerá se uma das funções ϕ1, ϕ2, . . . , ϕm for redundante (poder ser vista como uma combinação linear das
outras, pelo menos nos pontos x1, x2, . . . , xn.). Para ver isso é necessário usar um pouco de álgebra linear,
fugindo um pouco do escopo de curso introdutório.
Estamos prontos para resolver o problema do início da discussão e usar o vetor medida para tentar
recuperar a funçãolei. Vamos começar criando a matriz Φ.
3
In [ ]: # Constroi uma matriz concatenando uma coluna de 1's de dimensão adequada, o vetor x e depois o vetor x^2.
= [ones(n_pontos) x x.^2];
In [13]: # Resolve a equação normal ' c = ' medida, usando o operador \ de Julia
c = ('*) \ ('*medida)
Out[13]: 3-element Array{Float64,1}:
1.71444
1.97455
3.13709
Veja como os valores se aproximam dos valores originais que eram 1.73 para o termos constante, 2 para
o coeficiente de x e π para o coeficiente associado a x2. Outra coisa interessante de fazer é mostrar o gráfico
da função final obtida e mostrar como ele se ajusta bem nos pontos medidos.
In [27]: # Define a função a partir dos coeficientes calculados e das funções originais 1, x, e x^2.
phi(x) = c[1].*1 + c[2].*x + c[3].*x.^2
# Plota os pontos medidos e o gráfico da função obtida.
scatter(x, medida)
plot(x, phi(x), linewidth=4, color="r", label="")
legend()
Out[27]: PyObject 
4
In [ ]:
In [ ]:
5
	Introdução
	Modelando o problema de ajuste polinômios de grau 2
	Ajuste linear de curvas
	Equações normais