(1)Metody Numeryczne Wykład 7 Aproksymacja Aproksymacja na dyskretnym zbiorze punktów Dany jest zbiór punktów {(xi, yi)}mi=0, należy znaleźć funkcj¸e v(x), która ”przybliża”te punkty

(1)

Metody Numeryczne Wykład 7 Aproksymacja Aproksymacja na dyskretnym zbiorze punktów

Dany jest zbiór punktów {(x_i, y_i)}^m_i=0, należy znaleźć funkcj¸e v(x), która ”przybliża”te punkty. Jak wprzypadku interpolacji zakładamy że funkcja ta ma postać

v(x) = c₀φ₀+ c₁φ₁+ . . . + c_nφ_n=

n

X

j=0

c_jφ_j(x)

Zbiór {φ_j(x)}ⁿ_j=0 jest zbiorem funkcji bazowych - jest wi¸ec liniowo niezależny. Poszu- kujemy nieznanych współczynników c_j, j = 0, . . . , n. Różnica mi¸edzy interpolacj¸a a aproksymacj¸a pollega na tym, że w przypadku interpolacji zakładaliśmy równość m = n.

Teraz możemy założyć, że n < m. Innymi słowy liczba punktów pochodz¸acych z analizy danych może przewyższać liczb¸e funkcji bazowych. Przypomnijmy z algebry liniowej, że norm¸e euklidesow¸a k k2 wektora x definiujemy jako

kxk₂ =

√ x^Tx =

v u u t

m

X

i=1

x²_i

Aproksymacja w sensie metody najmniejszych kwadratów

Jest to najprostszy, a zarazem najpi¸ekniejszy matematycznie rodzaj aproksymacji.

Definiujemy współczynniki ai,j = φ_j(x_i) dla j = 0, 1, . . . , n, i = 0, 1, . . . , m.

Wówczas

v(xi) =

n

X

j=0

cjφj(xi) =

n

X

j=0

cjai,j.

Zadanie najmiejszych kwadratów polega znalezieniu współczynników c_j, które minimali- zuj¸a funkcj¸e

v u u t

m

X

j=0

y_i−

n

X

j=0

a_i,jc_j

!2

co możemy zapisać wektorowo

min_cψ(c), gdzie ψ(c) = ky − Ack²₂

(2)

Warunkiem koniecznym istnienia minimum tej funkcji jest zerowanie si¸e pochodnych rz¸edu pierwszego:

ψ|k(c) = 0, k = 0, 1, . . . , n

Chc¸ac być ścisłym matematycznie należałoby wykazać, że macierz pochodnych cz¸askowych rz¸edu drugiego (Hesjan) φ|k|l(c), k, l = 0, 1, 2 . . . , n jest macierz¸a dodatnio określon¸a.

Wówczas na wektorze c istnieje minimum lokalne funkcji φ.

Prosz¸e sprawdzić, że macierz φ|k|l(c) jest macierz¸a dodatnio określon¸a.

Ponieważ

ψ(c) = krk² =

m

X

i=0

y_i −

n

X

j=0

a_i,jc_j

!2

, wi¸ec

ψ|k(c) = 2

m

X

i=0

y_i−

n

X

j=0

a_i,jc_j(−a_i,k)

!

= 0, dla k = 0, 1, . . . , n Układ ten to możemy zapisać w postaci

m

X

i=0

a_i,k

n

X

j=0

a_i,jc_j =

m

X

i=0

a_i,ky_i, dla k = 0, 1, . . . , n.

lub w postaci macierzowej

A^TAc = A^Ty.

Jest to tak zwany układ równań normalnych . Układ ten możemy zapisać w postaci:

Bc = b, gdzie B = A^TA, b = A^Ty

Zauważmy, że macierz B = A^TA jest macierz symetryczn¸a, dodatnio określon¸a tzn. dla każdego niezerowego wektora x

x^TBx = x^T(A^TA)x = (Ax)^T(Ax) = kAxk²₂ > 0 Przykład (Regresja liniowa )

Poszukujemy funkcji v(x) aproksymuj¸acej zbiór danych punktów w postaci linii prostej

(3)

v(x) = c₀φ₀(x) + c₁φ₁(x) = c₀+ c₁x,

W tym przypadku n = 1, funkcjami bazowymi s¸a jednomian stały φ0(x) = 1,i jednomian liniowy φ₁(x) = x. Obliczamy wartości współczynników c_i, i = 0, 1 za pomoc¸a układu równań normalnych.

Macierz układu B = A^TA ma postać

B_0,0 B_0,1 B_1,0 B_1,1

gdzie B_0,0 =Pm

i=0φ₀(x_i)φ₀(x_i) = m + 1;

B_0,1 =Pm

i=0φ₀(x_i)φ₁(x_i) =Pm i=0x_i; B_1,0 =Pm

i=0φ₁(x_i)φ₀(x_i) =Pm

i=0x_i = B_0,1; B_1,1 =Pm

i=0φ₁(x_i)φ₁(x_i) =Pm i=0x²_i;

Wektor kolumnowy niewiadomych współczynników c c = c₀

c₁

Wektor kolumnowy wyrazów wolnych b = A^Ty.

b = b₀ b₁

gdzie b0 =Pm

0 yiφ0(xi) = Pm i=0yi

b₁ =Pm

0 y_iφ₁(x_i) = Pm i=0x_iy_i

St¸ad otrzymujemy nast¸epuj¸acy układ równań dla zadania regresji linniowej:

m + 1 Pm i=0x_i Pm

i=0x_i Pm i=0x²_i

c₁ c₂

=

Pm i=0y_i Pm

i=0x_iy_i

Nie b¸edziemy wypisywali jawnych wzorów na nieznane współczynniki c_i, i = 0, 1. Po- każemy na przykład dla danych (xi, yi) : (0.0, 0.1), (1.0, 0.9), (2.0, 2.0), w jaki sposób

(4)

konstruuje si¸e ten układ.

m = 2, m + 1 = 3, P3

i=0x_i = 0.0 + 1.0 + 2.0 = 3.0, Pm

i=0y_i = 0.1 + 0.9 + 2.0 = 3.0, Pm

i=1x²_i = 0.0²+ 1.0²+ 2.0² = 5, Pm

i=0x_iy_i = 0.0 · 0.1 + 1.0 · 0.9 + 2.0 · 2.0 = 4.9.

St¸ad otrzymujemy układ

3.0 3.0 3.0 5.0

c₀ c₁

3.0 4.9

Rozwi¸azuj¸ac go w OCTAVE:

» A = [3 3; 3 5];

» b = [3 4.9]⁰;

» c = A\b Otrzymujemy c₀ = 0.05 c₁ = 0.95.

Prost¸a regresji dla tych danych jest prosta o równaniu y = 0.05 + 0.95x.

Zajmiemy si¸e teraz zadaniem najmniejszych kwadratów w OCTAVE.

Standardow¸a metod¸a aproksymacji w OCTAVE jest aproksymacja metod¸a najmniejszych kwadratów wielomianami wybranego stopnia n.

Polecenie

a = polyf it(x, y, n)

znajduje wektor a współczynników wielomianu stopnia n najlepiej dopasowanego w sen- sie aproksymacji średniokwadratowej do danych wektorów x i y.

Wartość wielomianu aproksymuj¸acego w dowolnym punkcie x0 można wyznaczyć, korzy- staj¸ac z polecenia OCTAVE polyval(a, x0).

Funkcje te wykorzystamy do rozwi¸azania zadań z listy 3 na laboratorium 4, teraz napi- szemy script f.m w OCTAVE do rozwi¸azywania układów równań normalnych o nazwie np. mnk.m od metody najmniejszych kwadratów.

f unction c = mnk(x, y, n) x = x⁰;

y = y⁰;

m = size(x, 1);

A = ones(m, n + 1);

f or j = 1 : n

A(1 : m, j + 1) = A(1 : m, j). ∗ x;

end

% Układ równań normalnych B = A ∗ A⁰;

b = A⁰∗ y;

c = B\b;

(5)

Interpretacja geometryczna metody najmniejszych kwadratów (patrz rysunek)

(6)

Rysunek 1: Metoda najmniejszych kwadratów