Aproksymacja i interpolacja

mgr inż.

Sebastian Dziarmaga-Działyński

Wprowadzenie

Aproksymacja oraz interpolacja, to dwa bardzo ważne zagadnienia przerabiane w ramach kursu badań operacyjnych. Nazwy wzbudzają u niektórych strach i niezrozumienie, a także są często mylone. Czy słusznie? W tym artykule postaram się przybliżyć te dwa zagadnienia.

Interpolacja

Interpolacja liniowa i przykład „życiowy”

Zaczniemy od interpolacji, jako od zagadnienia prostszego, a nawet bliższego codziennemu życiu. Czym zajmuje się interpolacja? Do czego służy?

Zacznijmy może od problemu życiowego. Załóżmy, że chcemy naszemu kotu, czy psu zaaplikować środek na pchły, czy kleszcze. Na ulotce mamy tabelkę z ilością preparatu, jaka powinna zostać zaaplikowana w zależności od wagi zwierzęcia. Załóżmy, że tabelka ta wygląda tak:

Oczywiście dane są całkowicie zmyślone. Ważymy naszego pupila – waga pokazuje 5,8 kg. Ile mililitrów preparatu zaaplikować? Czy potraktować to, jako 5 kg i zaaplikować 16 ml, a może już jako 7 kg i zaaplikować 21 ml? Co, jeśli 16 ml okaże się za mało i nie pokonamy robactwa a 21 ml to za dużo i zaszkodzimy pupilowi?

Co w takiej sytuacji zazwyczaj robimy? Waga 5,8 kg to waga pośrednia pomiędzy zawartymi w tabeli wartościami 5 kg i 7 kg, więc poszukujemy jakiejś pośredniej wartości mililitrów. Staramy się to robić zgodnie z jakąś logiką i kombinujemy tak:

Przy 5 kg dawka to 16 ml, przy wadze o 2 kg większej, dawka to 21 ml. A zatem tak, jakby w zakresie między 5 a 7 kg masy zwierzęcia, na 2 kg wagi przypadało 5 ml. Nasz pupil waży o 0,8 kg więcej, czyli jakby powinien dostać dawkę 16 ml powiększoną o część z tych 5 ml i to proporcjonalnie taką część, z tych 5ml, jaką częścią 2 kg jest nasze 0,8 kg.

Innymi słowy powinien dostać 16 ml + 0,8/2*5 ml, co daje nam 2 ml. Czyli optymalna dawka leku dla naszego zwierzaka, to 18 ml.

Brawo! Dokonaliśmy właśnie najprostszej interpolacji – interpolacji liniowej. Założyliśmy sobie, że pomiędzy sąsiednimi wartościami masy ciała w tabeli (5 kg oraz 7 kg), dawka leku od 16 ml do 21 ml zmienia się liniowo. To, co zrobiliśmy, na rysunku można przedstawić następująco.

Interpolacja wielomianowa

W przytoczonym przypadku, interpolację dokonaliśmy de facto na dwu punktach (5; 16) oraz (7; 21). Pozostałe punkty z tabeli w ogóle nas nie interesowały. Zrobiliśmy to po to, by poznać wartość zmiennej $y$ (dawki lekarstwa dla zwierzaka), dla wartości $x$ (masy jego ciała), która nie została wymieniona w tabeli.

Na zajęciach z metod numerycznych najczęściej interpolacji dokonuje się na większej ilości punktów. Mamy kilka punktów $(x_i; y_i)$ i chcemy poznać wartości $y$ dla $x$ niewymienionego w tabeli. W tym celu szukamy takiej funkcji która – to bardzo ważne, kluczowe dla interpolacji i odróżniające ją od aproksymacji – dla danych wartości $x_i$ da dokładnie takie same wartości $y_i$ i, niejako przy okazji, pozwoli na wspomniane wyliczenie $y$ dla wartości pośrednich.

Funkcją, która pod względem matematycznym dobrze nadaje się do spełnienia warunku idealnego dopasowania swych wartości do danych punktów, jest funkcja wielomianowa. Dla przypomnienia, wielomianem (stopnia $n$) nazywamy taką funkcję:

$$W_n(x) = a_0 + a_1 x + a_2 x^2 + ... + a_{n-1} x^{n-1} + a_n x^n \tag 1 \label {eq:1}$$

$a_0, a_1, ..., a_n$ są to współczynniki wielomianu. Interpolacja liniowa, o której mówiliśmy na początku, polegała na znalezieniu prostej przechodzącej przez dwa punkty. Matematycznie prostą, reprezentuje równanie prostej:

$$y = ax + b$$

poznane bodajże jeszcze w szkole podstawowej, które jest niczym innym, jak właśnie wielomianem. Wielomianem pierwszego stopnia, gdyż zmienna $x$ występuje tam w potędze pierwszej. W naszym przypadku nie wyznaczaliśmy tego równania prostej w sposób jawny, ale gdybyśmy chcieli, to jest to jak najbardziej możliwe. Póki co nie angażujmy do tego sposobów przerabianych na studiach, ale użyjmy wiedzy z poziomu szkoły średniej. Mamy dane dwa punkty (5; 16) oraz (7; 21). Korzystamy ze wzoru (zawartego w każdych tablicach „maturalnych”):

$$y - y_1 = (x - x_1) \cdot {{y_2 - y_1} \over {x_2 - x_1}}$$

u nas:

$$y - 16 = (x - 5) \cdot {{21 - 16} \over {7 - 5}}$$

po przekształceniach:

$$y = 2,5x + 3,5$$

Zatem nasz wielomian interpolacyjny, dla punktów $(5; 16)$ oraz $(7,21)$, wyznaczony metodą „szkolną” jest następujący:

$$W(x) = 2,5 x + 3,5$$

Sprawdźmy, czy się wszystko zgadza. Zobaczmy, jaki $y$ wyjdzie, dla $x = 5$:

$$W(5) = 2,5 \cdot 5 + 3,5 = 12,5 + 3,5 = 16$$

a dla $x = 7$:

$$W(7) = 2,5 \cdot 7 + 3,5 = 17,5 + 3,5 = 21$$

Wszystko się zgadza. Nasza „pośrednia dawka dla pośredniej wagi”:

$$W(5,8) = 2,5 \cdot 5,8 + 3,5 = 14,5 + 3,5 = 18$$

Tutaj też się zgadza.

Oczywiście nasz wielomian nie da wartości $y$ z tabelki dla innych $x$ z tej tabelki, bo oczywiście uwzględniliśmy tylko te dwa punkty.

Co by było, gdybyśmy chcieli, aby nasz wielomian „zgadzał się” dla wszystkich punktów z tabelki – a są one cztery? Zastanówmy się przede wszystkim, wielomianu jakiego stopnia poszukujemy. Zachodzi tutaj prosta zależność – jeśli mamy $n$ punktów, to potrzebujemy wielomianu stopnia $n-1$. No bo tak, jeśli są dwa punkty, to jak wiemy, można przez nie poprowadzić prostą (czyli wielomian stopnia 1). Jeśli mielibyśmy trzy punkty, to przez trzy punkty można jednoznacznie poprowadzić parabolę (wielomian stopnia drugiego). Z kolei wielomian stopnia zerowego, czyli funkcję stałą (jej wykres to prosta pozioma) można poprowadzić przez jeden punkt. I tak dalej.

Uwaga – w tym momencie trzeba zwrócić uwagę na bardzo ważną rzecz. Aby interpolacja była możliwa, dane punkty muszą odpowiadać funkcji, czyli jednej wartości zmiennej $x$ może odpowiadać tylko jedna wartość $y$. Gdybyśmy mieli np. punkty (3; 5) oraz (3; 6) to nie dałoby się dokonać interpolacji. Nadal dałoby się przez takie punkty poprowadzić prostą, ale byłaby to prosta pionowa, o równaniu $x = 3$, a takie równanie nie jest równaniem funkcji liniowej. Nic nie stoi natomiast na przeszkodzie, aby ten sam $y$ odpowiadał różnym wartościom $x$. Dla punktów (2; 4), (3; 5), (4; 4) bez problemu można znaleźć wielomian interpolacyjny.

U nas zatem, skoro mamy cztery punkty, potrzebujemy wielomianu stopnia trzeciego. Czyli takiej funkcji:

$$W_3(x) = a_0 + a_1 x + a_2 x^2 + a_3 x^3$$

Tak naprawdę, to właściwsze byłoby stwierdzenie, że potrzebujemy wielomianu stopnia co najwyżej trzeciego. Bo wyobraźmy sobie, że np. te cztery punkty akurat leżą na jednej prostej – wtedy wystarczyłoby równanie prostej. Matematyka sobie z taką sytuacją poradzi w taki sposób, że po prostu współczynniki $a_2$ oraz $a_3$ wyszłyby równe zero.

Spróbujmy wyznaczyć sobie ten wielomian interpolacyjny nie korzystając ze wzorów przerabianych na kursie metod numerycznych. Jak wiadomo, kluczowe dla interpolacji jest, by „dla podanych iksów wychodziły podane igreki”. Innymi słowy, uwzględniając punkty z naszej tabeli, muszą być spełnione warunki: $W_3(3) = 11$, $W_3(5) = 16$, $W_3(7) = 21$, $W_3(10) = 25$. Podstawiając te zależności do wzoru (1) dostaniemy układ równań:

$$\left\{ {\begin{aligned} a_0 + a_1 \cdot 3 + a_2 \cdot 3^2 + a_3 \cdot 3^3 = 11 \\ a_0 + a_1 \cdot 5 + a_2 \cdot 5^2 + a_3 \cdot 5^3 = 16 \\a_0 + a_1 \cdot 7 + a_2 \cdot 7^2 + a_3 \cdot 7^3 = 21 \\a_0 + a_1 \cdot 10 + a_2 \cdot 10^2 + a_3 \cdot 10^3 = 25 \end{aligned}}\right.$$

w którym niewiadomymi są współczynniki $a_0$, $a_1$, $a_2$, $a_3$. I o to chodzi! Będą współczynniki, będziemy mieć wielomian. Współczynniki są cztery, równania też cztery. Jest więc szansa na rozwiązanie.

Zapiszmy układ w „ładniejszy” sposób:

$$\left\{ {\begin{aligned} a_0 + 3 a_1 + 9 a_2 + 27 a_3 = 11 \\ a_0 + 5 a_1 + 25 a_2 + 125 a_3 = 16 \\a_0 + 7 a_1 + 49 a_2 + 343 a_3 = 21 \\a_0 + 10 a_1 + 100 a_2 + 1000 a_3 = 25 \end{aligned}}\right.$$

Do rozwiązania użyć można dowolnej metody. W tym wypadku użyjemy metody macierzowej.

Macierz współczynników ma postać:

$$\left[ \begin{matrix} 1 & 3 & 9 & 27 \\ 1 & 5 & 25 & 125 \\1 & 7 & 49 & 343 \\1 & 10 & 100 & 1000 \end{matrix} \right]$$

Jako, że tematem artykułu nie jest rozwiązywanie układów równań liniowych, przeto nie będziemy dokładnie opisywać metodyki rozwiązania układu – ograniczymy się tylko do wzoru na wektor niewiadomych:

$$\textbf A = \left[ \begin{matrix} a_0 \\ a_1 \\ a_2 \\ a_3 \end{matrix} \right] = \left[ \begin{matrix} 1 & 3 & 9 & 27 \\ 1 & 5 & 25 & 125 \\1 & 7 & 49 & 343 \\1 & 10 & 100 & 1000 \end{matrix} \right]^{-1} \cdot \left[ \begin{matrix} 11 \\ 16 \\ 21 \\25 \end{matrix} \right] $$

W rezultacie otrzymujemy:

$$\textbf A = \left[ \begin{matrix} 7 \\ 2 \over 15\\ 1 \over 2 \\ -{1 \over 30} \end{matrix} \right] $$

Poszukiwany wielomian zatem ma następująca postać:

$$W_3(x) = 7 + {2 \over 15} x + {1 \over 2} x^2 - {1 \over 30} x^3 \tag 2 \label {eq:2}$$

Bez problemu można sprawdzić, że:

$$W_3(3) = 7 + {2 \over 15} \cdot 3 + {1 \over 2} \cdot 9 - {1 \over 30} \cdot 27 = 11$$

$$W_3(5) = 7 + {2 \over 15} \cdot 5 + {1 \over 2} \cdot 25 - {1 \over 30} \cdot 125 = 16$$

$$W_3(7) = 7 + {2 \over 15} \cdot 7 + {1 \over 2} \cdot 49 - {1 \over 30} \cdot 343 = 21$$

$$W_3(10) = 7 + {2 \over 15} \cdot 10 + {1 \over 2} \cdot 100 - {1 \over 30} \cdot 1000 = 25$$

Na poniższym wykresie zaprezentowano przebieg wyliczonego wielomianu interpolacyjnego.

Jak widać, w zakresie wartości zmiennej $x$ pomiędzy najmniejszą a największą wartością z tabeli, nie można mieć do naszego wielomianu większych zastrzeżeń. Łatwo policzyć, że dla zakładanej wstępnie wartości 5,8 kg, nasz wielomian osiąga wartość (wynik zaokrąglono):

$$W_3(5,8) = 18,09$$

otrzymano więc wartość nieznacznie wyższą od uzyskanej poprzednio.

Jak widać na wykresie, poza zakresem zmienności zmiennej $x$ z danych, wykres zupełnie się „rozjeżdża” i bez jakiejkolwiek logiki, pomimo wcześniejszego wzrostu, natychmiast po przekroczeniu przez $x$ wartości 10, wartość funkcji gwałtownie spada… Nieco inaczej będzie to wyglądało w przypadku aproksymacji, ale o tym powiemy w dalszej części.

Interpolacja Lagrange’a

W ramach kursu badań operacyjnych zazwyczaj nie wyznacza się współczynników wielomianu interpolacyjnego poprzez rozwiązanie układu równań liniowych. Najczęściej korzysta się z wielomianów interpolacyjnych Lagrange’a oraz Newtona.

Wielomian interpolacyjny Lagrange’a wylicza się ze wzoru:

$$L_n(x) = \sum_{i=0}^n {{\prod_{j \neq i} \left(x - x_j \right)} \over {\prod_{j \neq i} \left(x_i - x_j \right)}}$$

We wzorze tym licznik sumy $i=0,1,2,...,n$ odpowiada kolejnym punktom interpolacyjnym – punktów tych jest $n+1$ (numerujemy od zera!) a wielomian jest stopnia $n$, co zgodne jest z poprzednią konkluzją, że stopień wielomianu interpolacyjnego jest o 1 niższy od ilości punktów.

Budowa wzoru jest prosta – ilość składników sumy równa jest ilości punktów interpolacyjnych. W liczniku od zmiennej $x$ odejmujemy wartości współrzędnej $x$ kolejnych punktów pomijając wartość współrzędnej $x$ punktu odpowiadającemu liczonemu składnikowi. Mianownik ma taką samą budowę, jak licznik, tyle, że w miejscu zmiennej $x$ znajduje się pominięta współrzędna.

W naszym przypadku ze wzoru skorzystamy w sposób następujący:

$$L_3(x) = {{(x - 5)(x-7)(x-10)} \over {(3 - 5)(3 - 7)(3 - 10)}} \cdot 11 + {{(x - 3)(x-7)(x-10)} \over {(5 - 3)(5 - 7)(5 - 10)}} \cdot 16 + $$
$$ + {{(x - 3)(x-5)(x-10)} \over {(7 - 3)(7 - 5)(7 - 10)}} \cdot 21 + {{(x - 3)(x-5)(x-7)} \over {(10 - 3)(10 - 5)(10 - 7)}} \cdot 25 \tag 3 \label {eq:3}$$

Wzór działa według prostej koncepcji. Dla danego punktu interpolacyjnego $x_i$ wyrażenia $x - x_i$ przy wszystkich pozostałych punktach (w liczniku) równe są zero a z kolei dla tego punktu licznik równy jest mianownikowi (ułamek staje się jedynką), dzięki czemu w prosty sposób osiąga się warunek $W_n(x_i) = y_i$.

Po wymnożeniu liczb oraz wyrażeń we wzorze (3) oraz zredukowaniu wyrażeń podobnych, otrzymujemy wielomian określony wzorem (2). Jest to normalne, gdyż może istnieć tylko jeden wielomian interpolacyjny stopnia co najwyżej $n$ dla danego zestawu punktów. 

Oczywiście można wykonać owe przekształcenia i doprowadzić wielomian do postaci ogólnej i często wymagają tego prowadzący zajęcia, ale generalnie idea wzoru Lagrange’a polega na tym, aby właśnie tego nie robić, tylko bez kłopotliwych przekształceń od razu skonstruować działający wielomian interpolacyjny.

Wielomian interpolacyjny Lagrange’a to nie jest nazwa jakiegoś specjalnego rodzaju wielomianu. To tylko nazwa sprytnego zapisu działającego wielomianu interpolacyjnego, bez wyznaczania jego współczynników. Gdybyśmy np. chcieli w jakimś języku programowania zapisać taki działający wielomian, to spokojnie wystarczyłoby zdefiniować funkcję wg wzoru (3). Oczywiście postać (2) jest ładniejsza i pewnie szybsza dla kompilatora/interpretera/parsera, ale obie postaci są matematycznie równoważne.

W zasadzie, moim zdaniem, nie powinno się mówić wielomian interpolacyjny Lagrange’a tylko wzór interpolacyjny Lagrange’a.

Interpolacja Newtona

Tutaj również mówi się o wielomianie interpolacyjnym Newtona, choć de facto jest to po prostu wzór na jeden i ten sam wielomian interpolacyjny. Skonstruowany według nieco innych zasad. Sam wzór posiada dość skomplikowany zapis, więc pokażemy od razu sposób konstrukcji wielomianu.

Skonstruowanie wielomianu według wzoru interpolacyjnego Newtona wymaga wyliczenia ilorazów różnicowych funkcji. Jeśli ktoś uważał na zajęciach z analizy matematycznej, to powinien pamiętać to pojęcie. Występuje ono przy definiowaniu pochodnej funkcji. Jest ona definiowana, jako granica ilorazu różnicowego, dla przyrostu funkcji dążącego do zera.

Jeśli mamy dwa punkty $x_0$ oraz $x_1$, gdzie $x_0 < x_1$ pomiędzy którymi funkcja nasza jest ciągła, to ilorazem różnicowym funkcji nazywamy stosunek zmiany wartości funkcji do przyrostu jej argumentu, tj.:

$$ {{\Delta f(x_0, x_1)} \over {\Delta x}} = {{f(x_1) - f(x_0)} \over {x_1 - x_0}} \tag 4 \label {eq:4}$$

Pojęcie ilorazu różnicowego można rozszerzyć na iloraz różnicowy wyższych rzędów. Iloraz różnicowy drugiego rzędu jest to po prostu iloraz różnicowy z ilorazu różnicowego pierwszego rzędu. Czyli mamy trzy punkty $x_0, x_1, x_2$, obliczamy ilorazy różnicowe między $x_0$ a $x_1$ oraz między $x_1$ a $x_2$ i dla tych ilorazów różnicowych obliczamy ilorazy różnicowe, czyli:

$${{\Delta^2 f(x_0, x_1,x_2)} \over {\Delta x^2}} = {{{\Delta f(x_1,x_2)} \over {\Delta x}} - {{\Delta f(x_0,x_1)} \over {\Delta x}} \over {x_2 - x_0}}$$

$$ {{\Delta^2 f(x_0, x_1,x_2)} \over {\Delta x^2}} = {{{f(x_2) - f(x_1)} \over {x_2 - x_1}} - {{f(x_1) - f(x_0)} \over {x_1 - x_0}} \over {x_2 - x_0}} \tag 5 \label {eq:5} $$

Analogicznie wyprowadza się ilorazy różnicowe kolejnych rzędów.

Kolejne ilorazy obliczamy w tabeli:

Czy ta $-{1 \over 30}$ coś Wam przypomina?

Teraz, to co najfajniejsze, to współczynniki we wzorze interpolacyjnym Newtona czytamy tak jakby „po trójkącie” i mamy dwie możliwości: po górnym, lub dolnym trójkącie. Najpierw wypiszemy wzór a potem będzie komentarz:

Pierwsza możliwość – „po górnym trójkącie”:

$$N_3(x) = 11 + {5 \over 2} \cdot (x - 3) + 0 \cdot (x - 3) (x-5) - {1 \over 30} \cdot (x - 3) (x-5) (x - 7)$$

Wartość 2,5 przedstawiono, jako $5 \over 2$, aby nie mieszać różnych postaci ułamka w jednym wzorze, a 0 zapisano, by było wiadomo, jak wzór powstał. Po jego wyeliminowaniu wzór przybiera postać:

$$N_3(x) = 11 + {5 \over 2} \cdot (x - 3) - {1 \over 30} \cdot (x - 3) (x-5) (x - 7) \tag 6 \label {eq:6}$$

Ale możemy też zapisać wzór „po dolnym trójkącie”:

$$N_3(x) = 25 + {4 \over 3} \cdot (x - 10) - {7 \over 30} \cdot (x - 7) (x - 10) - {1 \over 30} \cdot (x - 5) (x - 7) (x - 10) \tag 7 \label {eq:7} $$

Zauważamy ten schemat? Wzór jakby narasta wraz z kolejnym składnikiem o wyrażenie $x-x_i$, gdzie $x_i$ jest współrzędną $x$ kolejnego punktu: licząc od początku, jeśli wybieramy górny trójkąt, albo od końca, jeśli wybieramy dolny.

Czytelnicy ze smykałką do programowania od razu zauważą, że wartości takiego wielomianu świetnie oblicza się w pętli, wykorzystując namnażanie przyrostowe – w Pascalu byłoby to:

iloczyn := iloczyn*(x - x[i])

czy też w C prościej:

iloczyn *= x - x[i];

Oczywiście wzory (6) oraz (7) po przekształceniu doprowadzą do wielomianu (2), co Czytelnik sprawdzić może samodzielnie. Z tego też powodu zarówno w postaci (6), jak i (7), występuje ten sam współczynnik $-{1 \over 30}$. Jak już bowiem powiedziano, wielomian interpolacyjny Newtona nie jest żadnym innym, specjalnego rodzaju wielomianem, ale jest to po prostu wzór interpolacyjny na wyznaczenie wielomianu interpolacyjnego bez de facto wyznaczania jego współczynników w postaci ogólnej (1).

Wzór interpolacyjny Newtona posiada także swoją specjalną wersję, przeznaczoną dla argumentów równoodległych. U nas nie miałby on zastosowania, gdyż pierwsze trzy punkty są wprawdzie równoodległe (3; 5; 7) – różnica między nimi wynosi 2 (tworzą więc ciąg arytmetyczny), ale kolejny punkt ma współrzędną $x$ równą 10 i z tego ciągu się wyłamuje. W tym artykule nie będziemy jednak omawiać tego wzoru – być może poświęcimy mu odrębny artykuł.

Istnieje wiele metod interpolacji i związanych z interpolacją, jak interpolacja funkcjami giętkimi (splajnami), interpolacja Hermite’a czy wybór węzłów interpolacyjnych w oparciu o wzór Czebyszewa. Również zasługują na inny, obszerny artykuł.

Aproksymacja

O co chodzi w aproksymacji?

Głównym tematem niniejszego artykułu nie jest szczegółowe omawianie zagadnień interpolacji czy aproksymacji, ale przede wszystkim uzmysłowienie podstawowej różnicy między tymi metodami.

W odróżnieniu bowiem od interpolacji, dla której najważniejszą kwestią było, aby wyznaczyć taką funkcję, która pozwoli dla zadanych współrzędnych $x_i$ otrzymać dokładne wartości $y_i$ – czyli (dla „wzrokowców”) aby otrzymać taką krzywą, która przejdzie przez wszystkie punkty, aproksymacja szuka takiej krzywej, która jak najlepiej dopasuje się do ułożenia tychże punktów niekoniecznie przez nie przechodząc.

Jak zauważyliśmy na wykresie sporządzonym dla danych z naszej tabeli, punkty układały się „mniej więcej” wzdłuż linii prostej – co jest logiczne, bo im więcej waży zwierzę, tym więcej preparatu należy mu zaaplikować. Wielomian interpolacyjny świetnie się do tych punktów dopasował, ale poza obszarem objętym interpolacją, zupełnie z owym spodziewanym liniowym wzrostem się rozjechał.

Aproksymacja tutaj, jak za chwilę zobaczymy, okaże się znacznie lepszym narzędziem. Polega ona bowiem na wyznaczeniu takiej funkcji aproksymującej $A(x)$, że łączna, szeroko pojęta, różnica (odległość) pomiędzy wartościami funkcji aproksymującej a aproksymowanej (czyli naszymi punktami $y_i$) będzie jak najmniejsza:

$$\sum_{i=0}^n d \left( f(x_i), A(x_i) \right) \rightarrow \min$$

Owa różnica, czy też odległość może być bardzo różnie zdefiniowana. Nie może być to po prostu różnica: $d \left( f(x_i), A(x_i) \right) = f(x_i) - A(x_i)$ gdyż wówczas nawet bardzo duże odległości ujemne znosiłyby się (kompensowały) z odległościami dodatnimi i dopasowanie wyszłoby fatalne. Już całkiem niegłupia byłaby wartość bezwzględna: $d \left( f(x_i), A(x_i) \right) = |f(x_i) - A(x_i)$|, ale ona z kolei jest funkcją „ciężką” obliczeniowo, choćby dlatego, że jej wykresem jest litera v (co oznacza, że nie jest różniczkowalna w zerze)¹.

Aproksymacja średniokwadratowa

Znacznie lepszą odległością jest kwadrat różnicy: $d \left( f(x_i), A(x_i) \right) = \left( f(x_i) - A(x_i) \right)^2$. Przyjęcie takiej funkcji prowadzi do najpopularniejszego wariantu aproksymacji, czyli aproksymacji średniokwadratowej. W aproksymacji tej funkcję $A(x)$ dopasowujemy w taki sposób, aby:

$$\sum_{i=0}^n \left( f(x_i) - A(x_i) \right)^2 \rightarrow \min \tag 8 \label {eq:8}$$

Aby troszkę wszystko „udziwnić” i utrudnić, często wprowadza się dodatek w postaci tzw. funkcji kary $w(x)$, wówczas wzór (8) przyjmuje postać:

$$\sum_{i=0}^n w(x_i) \left( f(x_i) - A(x_i) \right)^2 \rightarrow \min \tag 9 \label {eq:9}$$

O co chodzi z funkcją kary? Otóż możemy sobie zdecydować, że nie wszystkie punkty są jednakowo ważne. Dokładniejsze dopasowanie do pewnych punktów uznać możemy za znacznie ważniejsze niż dopasowanie do innych. Aby to osiągnąć, po prostu tym ważniejszym punktom przypisujemy wyższa wartość funkcji kary niż tym mniej ważnym i aparat matematyczny metody, o którym wkrótce napiszemy, minimalizując wyrażenie (9) zadba o to, by uczynić zadość naszym preferencjom.

Póki co jednak, dla uproszczenia, załóżmy, że nie uwzględniamy funkcji kary, czyli przyjmujemy $w(x_i) = 1$, a zatem wszystkie punkty są dla nas jednakowo ważne. Jak wyznaczyć funkcję $A(x)$?

Aby wyznaczyć funkcję aproksymującą $A(x)$, należy w pierwszej kolejności ustalić postać tej funkcji. W przypadku aproksymacji zakładaną postać funkcji aproksymującej definiuje się w oparciu o tzw. funkcje bazowe. Definicja taka wygląda następująco:

$$A(x) = a_0 \cdot \varphi_0(x) + a_1 \cdot \varphi_1(x) + ... + a_m \cdot \varphi_m(x) \tag {10} \label {eq:{10}}$$

gdzie $\varphi_j(x)$ są przyjętymi funkcjami bazowymi, natomiast $a_j$ to poszukiwane wartości współczynników.

Wygląda to groźnie, ale w praktyce jest znacznie łatwiejsze. Przypuśćmy, że chcemy przyjąć $A(x)$ jako funkcję liniową (czyli wielomian stopnia pierwszego), tj:

$$A(x) = a_0 + a_1 x \tag {11} \label {eq:{11}}$$

wówczas, porównując wyrażenia (10) oraz (11) zauważamy, że funkcjami bazowymi będą: $\varphi_0 (x) = 1$ oraz $\varphi_1 (x) = x$. Możemy przyjmować najróżniejsze postaci funkcji bazowych. Wybór uzasadniony powinien być albo układem punktów empirycznych, a najlepiej powinien być uzasadniony merytorycznie. Przykładowo jeśli punkty empiryczne są np. wynikami pomiaru wydłużenia próbki materiału pod wpływem działającej siły, spodziewamy się zależności liniowej i  taką też funkcję dobieramy. Jeśli jest to np. pomiar temperatury stygnącego obiektu, dokonywany w określonych odstępach czasu, wówczas należałoby dobrać funkcję wykładniczą, itd.

W przypadku danych przedstawionych w naszej wyjściowej tabeli, wydaje się, że najwłaściwszą funkcją będzie właśnie funkcja liniowa.

Współczynniki $a_0, a_1, ..., a_m$ funkcji aproksymującej oblicza się najczęściej wzorami macierzowymi.

Zdefiniujmy następujące macierze:

$\mathbf F = \left[ \begin{matrix} \varphi_0(x_0) & \varphi_1(x_0) & ... & \varphi_m(x_0) \\ \varphi_0(x_1) & \varphi_1(x_1) & ... & \varphi_m(x_1) \\ ... & ... & ... & ... \\ \varphi_0(x_n) & \varphi_1(x_n) & ... & \varphi_m(x_n) \end {matrix} \right]$    $\boldsymbol y = \left[ \begin {aligned} y_0 \\y_1 \\... \\y_n \end {aligned} \right]$   

$\mathbf F$ jest macierzą wartości funkcji bazowych dla naszych punktów empirycznych, natomiast wektor $\boldsymbol y$ jest wektorem współrzędnych $y$ naszych punktów. 

Wektor $\boldsymbol a$ współczynników funkcji aproksymującej $A(x)$ obliczamy ze wzoru:

$$\boldsymbol a = \left[ \begin {aligned} a_0 \\ a_1 \\ ... \\ a_m \end {aligned} \right] = \left( {\mathbf F}^{\top} {\mathbf F} \right)^{-1} \cdot {\mathbf F}^{\top} {\boldsymbol y} \tag {12} \label {eq:{12}}$$

Jeśli ktoś zetknął się kiedyś z ekonometrią, wzór (12) może wydawać mu się dziwnie znajomy. Skojarzenie jest jak najbardziej trafne! Aproksymacja średniokwadratowa i estymacja parametrów równania regresji wielorakiej metodą najmniejszych kwadratów, to – z matematycznego punktu widzenia – praktycznie jedno i to samo. Różnice są przede wszystkim natury merytorycznej. Przy aproksymacji nie analizuje się właściwości reszt (błędów aproksymacji) pod kątem stałości wariancji czy autokorelacji. Wyznaczonych wartości współczynników funkcji aproksymującej $A(x)$ nie traktuje się jak estymatory, gdyż nie traktujemy punktów empirycznych jako próby losowej.

Regresja wieloraka kładzie większy nacisk na zależność przyczynowo-skutkową pomiędzy zmienną zależną a zmiennymi objaśniającymi a w aproksymacji liczy się przede wszystkim dopasowanie matematyczne. W regresji wielorakiej kolumny macierzy $\mathbf F$, zwanej tam zazwyczaj macierzą $\mathbf X$ są wartościami kolejnych zmiennych objaśniających, podczas gdy w aproksymacji są to wartości funkcji bazowych wyliczone dla tego samego zestawu współrzędnych $x$. Można powiedzieć, że aproksymacja to taka regresja wieloraka, w której zmiennymi objaśniającymi są funkcje bazowe. Jak powiedziano, w zastosowaniach praktycznych czasem jednak się zdarza, że dobór poszczególnych funkcji bazowej wynika z przesłanek merytorycznych.

Wyliczmy zatem postać liniowej funkcji aproksymującej średniokwadratowo dane z naszej tabeli. Funkcje bazowe, jak wspomniano, mieć będą postać $\varphi_0 (x) = 1$ oraz $\varphi_1 (x) = x$. Nasze macierze równe będą:

$\mathbf F = \left[ \begin{matrix} 1 & 3 \\ 1 & 5 \\ 1 & 7 \\ 1 & 10 \end {matrix} \right]$    $\boldsymbol y = \left[ \begin {aligned} 11 \\ 16 \\ 21 \\ 25 \end {aligned} \right]$  

Obliczamy wartości poszczególnych iloczynów macierzy:

$${\mathbf F}^{\top}{\mathbf F} = \left[ \begin{matrix} 1 & 1 & 1 &1 \\ 3 & 5 & 7 & 10 \end {matrix} \right] \cdot \left[ \begin{matrix} 1 & 3 \\ 1 & 5 \\ 1 & 7 \\ 1 & 10 \end {matrix} \right] = \left[ \begin{matrix} 4 & 25 \\ 25 & 183 \end{matrix} \right]$$

$${\mathbf F}^{\top}{\mathbf y} = \left[ \begin{matrix} 1 & 1 & 1 &1 \\ 3 & 5 & 7 & 10 \end {matrix} \right] \cdot \left[ \begin{matrix} 11 \\ 16 \\ 21 \\ 25 \end {matrix} \right] = \left[ \begin{matrix} 73 \\ 510 \end{matrix} \right]$$

Stosujemy wzór (12) :

$$\boldsymbol a = \left[ \begin{matrix} 4 & 25 \\ 25 & 183 \end{matrix} \right]^{-1} \cdot \left[ \begin{matrix} 73 \\ 510 \end{matrix} \right]$$

Obliczenie macierzy odwrotnej do macierzy 2×2 jest bardzo proste, więc możemy je zaprezentować:

$$\left[ \begin{matrix} 4 & 25 \\ 25 & 183 \end{matrix} \right]^{-1} = {1 \over {4 \cdot 183 - 25 \cdot 25}} \cdot \left[ \begin{matrix} 183 & -25 \\ -25 & 4 \end{matrix} \right] $$

$$ = {1 \over 107} \cdot \left[ \begin{matrix} 183 & -25 \\ -25 & 4 \end{matrix} \right] $$

$${\boldsymbol a} = {1 \over 107} \cdot \left[ \begin{matrix} 183 & -25 \\ -25 & 4 \end{matrix} \right] \cdot \left[ \begin{matrix} 73 \\ 510 \end{matrix} \right] $$

$${\boldsymbol a} = \left[ \begin{matrix} 609 \over 107 \\ 215 \over 107 \end{matrix} \right] = \left[ \begin{matrix} 5,692 \\ 2,009 \end{matrix} \right] $$

Funkcja aproksymująca ma postać:

$$A(x) = 5,692 + 2,009x$$

Dopasowanie wyznaczonej linii pokazuje poniższy wykres:

Tym razem, poza zakresem zmienności zmiennej $x$, funkcja nie zmienia tendencji, co jest oczywiste w przypadku funkcji liniowej.

Obliczmy wartość funkcji aproksymującej dla $x = 5,8$:

$$A(5,8) = 5,692 + 2,009 \cdot 5,8 = 17,34$$

Wartość ta nieco niższa niż otrzymana w wyniku interpolacji liniowej przeprowadzonej między punktami $(5; 16)$ oraz $(7; 21)$. Na wyraźnie widać przyczynę takiego stanu rzeczy. Odcinek łączący wskazane punkty znajduje się powyżej wykresu wyznaczonej funkcji aproksymującej.

Dla porównania wyliczmy wartość funkcji aproksymującej w postaci wielomianu drugiego stopnia. Funkcje bazowe mają w tym wypadku postać: $\varphi_0(x) = 1$, $\varphi_1(x) = x$, $\varphi_2(x) = x^2$.

Macierz $\mathbf F$ ma w tym wypadku postać:

$$\mathbf F = \left[ \begin{matrix} 1 & 3 & 9 \\ 1 & 5 & 25 \\ 1 & 7 & 49 \\ 1 & 10 & 100\end {matrix} \right]$$

Wyliczona postać wielomianu aproksymującego:

$$A(x) = 0,495 + 3,888x - 0,143x^2$$

Wykres dopasowania obu funkcji pokazano na poniższym rysunku. 

Widać tu podobną „słabość”, jak w przypadku interpolacji. Wielomian drugiego stopnia świetnie dopasowuje się do punktów empirycznych, niemalże przez nie przechodząc (wyliczone wartości funkcji aproksymującej wynoszą $A(3) = 10,87$, $A(5) = 16,36$, $A(7) = 20,70$, $A(10) = 25,06$ czyli są bardzo zbliżone do wartości $y_0 = 11$, $y_1 = 16$, $y_2 = 21$, $y_3 = 25$.

Aproksymacja – interpolacją

Ale do najciekawszych konkluzji dochodzimy, gdy próbujemy znaleźć funkcję aproksymującą w postaci wielomianu trzeciego stopnia. Funkcjami bazowymi będą $\varphi_0(x) = 1$, $\varphi_1(x) = x$, $\varphi_2(x) = x^2$, $\varphi_3(x) = x^3$.

Macierz $\mathbf F$ ma wówczas postać:

$$\mathbf F = \left[ \begin{matrix} 1 & 3 & 9 & 81 \\ 1 & 5 & 25 & 125 \\ 1 & 7 & 49 & 343 \\ 1 & 10 & 100 & 1000\end {matrix} \right]$$

Gdzieś taką macierz już chyba widzieliśmy… Tak! Jest to przecież macierz współczynników z zadania interpolacji wielomianowej! Dochodzimy tu bowiem do najważniejszego odkrycia. Otóż w pewnej szczególnej sytuacji, a mianowicie: gdy ilość wyznaczanych parametrów funkcji aproksymującej równa jest ilości punktów empirycznych, wtedy aproksymacja staje się interpolacją.

To niesamowite. Przekonywaliśmy, że to dwa odrębne zagadnienia, że się różnią, a jednak. Jednak posiadają one pewien „wspólny mianownik” w postaci tego, że czasem po prostu aproksymacja staje się interpolacją. Dlaczego tak? Zarówno funkcja interpolująca, jak i aproksymująca posiadają parametry – tutaj, były nimi współczynnik wielomianu. Działają one, jak pokrętła sterujące. Zmieniając wartości poszczególnych współczynników można, coraz bardziej dopasować funkcję do przebiegu punktów empirycznych, tak samo, jakbyśmy kręcili pokrętłami skomplikowanej maszyny.

Gdy pokręteł jest mniej niż punktów, zazwyczaj nie damy rady idealnie dopasować do wszystkich punktów. Kręcąc dwoma pokrętłami, gdy punkty są cztery, zawsze coś poprawiając, gdzie indziej coś zepsujemy. Mamy wówczas aproksymację – dopasowanie jak najlepsze, ale nieidealne. 

Gdy jednak ilość pokręteł jest równa ilości punktów można sobie wyobrazić, że jednym pokrętłem dopasowujemy do jednego punktu, innym do innego – dając w ten sposób radę dopasować do wszystkich. Mamy interpolację.To oczywiście spore uproszczenie, ale dobrze oddaje matematyczny sens opisywanej sytuacji.

Z drugiej strony idealność dopasowania nie jest celem samym w sobie, bo najczęściej płacimy za nią całkowicie bezsensownym przebiegiem funkcji poza oknem interpolacji. Można powiedzieć, że aproksymacja, starając się kilkoma tylko parametrami dopasować do wielu punktów, robi coś uniwersalnego – dopasowanie na tyle dobre, że jeśli punktów dołożymy, to nadal efekt jej działania będzie miał sens. Interpolacja natomiast kreuje „fachowca od jednej śrubki”. Funkcję, która tak koncentruje się na dokładnym dopasowaniu do wszystkiego, co jest, że zupełnie się kompromituje, gdy coś nowego przybędzie.

Omawiając interpolację zaznaczyliśmy, że interpolować możemy tylko wówczas, gdy wśród punktów empirycznych nie ma dwu (lub więcej) takich, które maja jednakową wartość współrzędnej $x$ a różne wartości zmiennej $y$. Aproksymacji to w żaden sposób nie przeszkadza, o ile nie jest to ten przypadek, gdy aproksymacja przeszła w interpolację.

Inne aproksymacje

Warto na zakończenie zasygnalizować, że korzystając z odpowiednich przekształceń, można aproksymować najróżniejszymi funkcjami. Przykładowo można zastosować potęgową funkcje aproksymującą:

$$A(x) = a_0 \cdot x^{a_1} \tag {13} \label {eq:{13}}$$

jakie będą wówczas funkcje bazowe? Aby je wyznaczyć, logarytmujemy obie strony równania (13):

$$\ln A(x) = \ln \left( a_0 \cdot x^{a_1} \right)$$

$$\ln A(x) = \ln a_0 + a_1 \cdot \ln x \tag {14} \label {eq:{14}}$$

Jeśli teraz podstawimy $B(x) = \ln A(x)$, $z = \ln x$, $b_0 = \ln a_0$, $b_1 = a_1$, otrzymamy:

$$B(z) = b_0 + b_z \cdot z$$

czyli najzwyklejszą aproksymację funkcją liniową. Po wyznaczeniu wartości parametrów $b_0$ oraz $b_1$ dokonujemy zwrotnego podstawienia $a_0 = e^{b_0}$, $a_1 = b_1$ i w ten sposób mamy wyznaczoną postać (13).

Innym ciekawym zagadnieniem związanym z aproksymacją, jest aproksymacja ciągła. Dotychczas omawiana aproksymacja byłą aproksymacją dyskretną, polegającą na jak najlepszym dopasowaniu funkcji aproksymującej do pewnej skończonej ilości punktów empirycznych. W aproksymacji ciągłej dopasowujemy funkcję nie do punktów a do innej ciągłej funkcji. Czyli łączna odległość liczona jest nie od skończonej ilości punktów, ale od całego kontinuum punktów tamtej funkcji. Ciągłym odpowiednikiem sumy jest całka, zatem odpowiednikami wzorów (8) i (9) dla interpolacji ciągłej na przedziale $\left \langle a, b \right \rangle$. są:

$$\int_a^b \left( f(x) - A(x) \right)^2 \,dx \rightarrow \min$$

oraz z funkcją kary:

$$\int_a^b w(x) \cdot \left( f(x) - A(x) \right)^2 \,dx \rightarrow \min $$

W ten sposób uogólnić można oczywiście tylko aproksymację. Ciągła wersja interpolacji nie miałaby sensu, gdyż po prostu funkcją interpolującą musiałaby być ta sama funkcja, co funkcja interpolowana – przynajmniej na wskazanym przedziale.

1. Niektórzy twierdzą, że wartość bezwzględna nie jest funkcją elementarną. Oczywiście powinni oni jak najszybciej wrócić do szkoły, albo zapisać się na nasze korepetycje, bo przecież $|x| = \sqrt {x^2}$, czyli jest jak najbardziej funkcją elementarną. ↩︎