Polynomická aproximace – cvičení

Dnešním tématem jsou numerické metody pro aproximaci funkcí – budeme tabelovanými body prokládat aproximační polynom, jednak přesně (polynomická interpolace), jednak přibližně (polynomická regrese).

Polynomická interpolace

Vstupními daty budou dvojice souřadnic bodů v rovině, \((x_i,y_i),\ i=0,\dots,n\). Cílem polynomické interpolace je sestrojit aproximační polynom proložený vstupními body přesně. Pro \(n+1\) bodů to bude polynom \(n\)-tého stupně,

\[P_n(x)=c_0+c_1x+\dots+c_nx^n\,.\]

Chceme tedy dosáhnout toho, aby \(P_n(x_i)=y_i\), což je vlastně \(n+1\) lineárních algebraických rovnic pro \(n+1\) neznámých koeficientů \(c_0,\dots,c_n\) aproximačního polynomu. Matici této soustavy se říká Vandermondova:

\[\begin{split}\left(\begin{array}{cccc} 1 & x_0 & \dots & x_0^n\\ 1 & x_1 & \dots & x_1^n\\ \vdots & & & \vdots\\ 1 & x_n & \dots & x_n^n \end{array}\right) \left(\begin{array}{c} c_0 \\ c_1 \\ \vdots \\ c_n \end{array}\right)= \left(\begin{array}{c} y_0 \\ y_1 \\ \vdots \\ y_n \end{array}\right).\end{split}\]

Vstupní data načteme ze souboru, přičemž jedním čtením zjistíme počet datových řádků n+1, poté alokujeme v potřebné velikosti pole x a y a druhým čtením do nich vložíme data. Koeficienty aproximačního polynomu získáme funkcí jménem polyfit. S názvem přišel Matlab, převzal ho pythonský modul NumPy a teď ho převezmeme i my. Nejprve mu dáme dvouargumentovou podobu polyfit2(x,y) pro výpočet koeficientů interpolačního polynomu stupně n, poté obvyklou podobu polyfit(x,y,m) pro výpočet koeficientů regresního polynomu stupně m. Aproximační polynom poté pomocí funkce polyval (jméno opět z Matlabu i NumPy) pro vyčíslení polynomu užitím Hornerova schématu (známe z cvičení o numerické integraci) tabelujeme na hustější síti xi překrývající interval se vstupními daty a spolu s nimi ho vykreslíme, a to jak v gnuplotu, tak pomocí pythonského modulu Matplotlib.

# Python a NumPy/mypoly: polynomická aproximace (interpolace a regrese)
import numpy               # standardní polyval a polyfit
import mypoly              # náš polyval, regresní polyfit a interpolační polyfit2
runMypoly=False            # False pro NumPy, True pro náš modul mypoly
m=1                        # stupeň aproximačního polynomu
f_in='data'                # vstupní soubor s (x,y)
f_out='fit'                # výstupní soubor s (xi,yi)

# Načtení dat z f_in
with open(f_in,'r') as f:  # zjištění počtu řádků ve vstupním souboru
  n=0
  for line in f: n+=1
n-=1                       # stupeň interpolačního polynomu
x=numpy.empty(n+1); y=numpy.empty(n+1)   # jednorázová alokace
with open(f_in,'r') as f:  # načtení vstupních dat
  nn=0
  for line in f: (x[nn],y[nn])=line.split(); nn+=1
print('Soubor',f_in,' - přečteno řádků:',nn)

if n==m:
# Polynomická interpolace
  if runMypoly: p=mypoly.polyfit2(x,y)   # náš interpolační polyfit2
  else:         p=numpy.polyfit(x,y,n)   # numpy.polyfit
else:
# Polynomická regrese
  if runMypoly: p=mypoly.polyfit(x,y,m)  # náš regresní polyfit
  else:         p=numpy.polyfit(x,y,m)   # numpy.polyfit

# Tabelování interpolačního polynomu
nmax=n*100; xmin=x[0]-1; xmax=x[n]+1
xi=numpy.linspace(xmin,xmax,nmax+1)
if runMypoly: yi=mypoly.polyval(p,xi)    # náš polyval
else:         yi=numpy.polyval(p,xi)     # numpy.polyval

# Zápis dat do f_out (pro gnuplot)
with open(f_out,'w') as f:
  for nn in range(nmax+1): print(f'{xi[nn]:10.6f} {yi[nn]:10.6f}',file=f)
print('Soubor',f_out,' - zapsáno řádků:',nmax+1)

# Kreslení dat pomocí Matplotlibu
import matplotlib.pyplot as plt
# fig=plt.figure(figsize=(16,8),layout='tight') # volitelná konfigurace obrázku
plt.plot(x,y,'o',xi,yi)           # symboly v (x,y), čára skrz (xi,yi)
plt.savefig('fit.png')            # uložení do PNG
# plt.show()                      # na obrazovku
# fig.savefig('fit.png')          # alternativní metoda pro uložení PNG

Zdrojový soubor mypoly.py s naší implementací funkce polyval a interpolační verzí polyfit2 řešící algebraickou soustavu s Vandermondovou maticí je zde:

# Modul (soubor mypoly.py) s naší implementací polyval a polyfit
import numpy
from numba import njit

# Vyčíslení polynomu p[0]*x**n + ... + p[n]
@njit
def polyval(p,x):
  result=numpy.empty(x.size)
  for nn in range(result.size):
    result[nn]=p[0]
    for n in range(1,p.size): result[nn]=result[nn]*x[nn]+p[n]
  return result

# Interpolační polynom p[0]*x**n + ... + p[n]
@njit
def polyfit2(x,y):
  n=x.size-1
  V=numpy.empty((n+1,n+1))         # Vandermondova matice
  # for i in range(n+1):           # pomocí cyklů
  #   V[i,0]=1
  #   for j in range(1,n+1): V[i,j]=V[i,j-1]*x[i]
  V[:,0]=1                         # totéž vektorově
  for j in range(1,n+1): V[:,j]=V[:,j-1]*x[:]
  c=numpy.linalg.solve(V,y)        # řešení algebraické soustavy
  return c[-1::-1]                 # zrcadlení vektoru

# ... druhá část souboru následuje v dalším odstavci

Funkce polyval i polyfit v Matlabu, Octavu i NumPy řadí koeficienty polynomů od členu s nejvyšší mocninou, což je opačně než řazení prvků vektoru řešení v naší soustavě. Výstupní vektor proto kvůli kompatibilitě zrcadlově převrátíme.

Polynomická regrese

Vstupní data mají často jen omezenou přesnost a není třeba nutit aproximační funkci jimi procházet přesně. Vyřešíme nyní frekventovaný požadavek proložit \(n+1\) body aproximační polynom stupně nižšího než \(n\), tedy vystihnout body jen přibližně. Cílem polynomické regrese je sestrojit v duchu metody nejmenších čtverců (LSQ) takový polynom \(P_m(x)\) stupně \(m\), aby součet druhých mocnin jeho odchylek od vstupních bodů byl minimální,

\[\sum_{k=0}^n\left[P_m(x_k)-y_k\right]^2 = \min\,.\]

V poznámkách k přednášce je tento vztah zderivován vzhledem k \(m+1\) koeficientům aproximačního polynomu a odvozena tak soustava \(m+1\) lineárních algebraických rovnic (tzv. normálních rovnic) pro hledané koeficienty. Matice soustavy se jmenuje Gramova:

\[\begin{split}\left(\begin{array}{llll} \sum_k 1 & \sum_k x_k & \dots & \sum_k x_k^m\\ \sum_k x_k & \sum_k x_k^2 & \dots & \sum_k x_k^{m+1}\\ \vdots & & & \vdots\\ \sum_k x_k^m & \sum_k x_k^{m+1} & \dots & \sum_k x_k^{2m} \end{array}\right) \left(\begin{array}{l} c_0 \\ c_1 \\ \vdots \\ c_m \end{array}\right)= \left(\begin{array}{l} \sum_k y_k \\ \sum_k y_kx_k \\ \vdots \\ \sum_k y_kx_k^m \end{array}\right) .\end{split}\]

Algoritmu dáme podobu tříargumentové funkce polyfit(x,y,deg), kde třetí argument stojí pro stupeň regresního polynomu m. K našemu modulu mypoly.py doplníme funkci:

# ... druhá část souboru mypoly.py s naší implementací funkce polyfit

# Regresní polynom p[0]*x**m + ... + p[m], m=deg
@njit
def polyfit(x,y,deg):
  m=deg; n=x.size-1
  # if m==n: return polyfit2(x,y)  # interpolační polynom
  xx=numpy.ones(n+1); sumx=numpy.empty(m*2+1); sumy=numpy.empty(m+1)
  for i in range(m*2+1):
    sumx[i]=sum(xx)
    if i<=m: sumy[i]=sum(y*xx)
    xx*=x
  G=numpy.empty((m+1,m+1))         # Gramova matice
  for i in range(m+1):
    for j in range(m+1): G[i,j]=sumx[i+j]
  c=numpy.linalg.solve(G,sumy)     # řešení algebraické soustavy
  return c[-1::-1]                 # zrcadlení vektoru

V Matlabu i Octavu, odkud funkce polyfit a polyval přišly, se to celé řekne docela stručně:

% Matlab/Octave: polynomická aproximace
m=1;                                                            % stupeň aproximačního polynomu
input='data'; output='fit';                                     % vstupní a výstupní soubor
[x,y]=textread(input,"%f %f");                                  % načtení dvousloupcových dat ze souboru
n=length(x);                                                    % počet dat
disp(['Soubor ',input,' - precteno radku: ',int2str(n)])        % výpis vektoru řetězců
p=polyfit(x,y,m);                                               % polynomická fitace dat
nmax=n*100; xi=linspace(x(1)-1,x(n)+1,nmax)'; yi=polyval(p,xi); % tabelování aproximačního polynomu do dvou sloupců
t=fopen(output,'w'); fprintf(t,"%f %f\n",[xi yi]'); fclose(t);  % zápis dvousloupcových dat do souboru
disp(['Soubor ',output,' - zapsano radku: ',int2str(nmax)])     % výpis vektoru řetězců
plot(x,y,'o',xi,yi)
saveas(gcf,'fit.png')

Přidáme i fortranskou ukázku, když už jsme se minule namáhali s přípravou modulu mLinalg s LAPACKovským řešičem algebraických soustav dgesv.

! Fortran a mypoly: polynomická aproximace (interpolace a regrese)
! f=polyval(p,x) pro vyčíslení polynomu p ve skaláru x nebo na poli x (Horner)
! p=polyfit(x,y,m) pro fitaci dat (x_i,y_i), i=0,..,n, polynomem stupně m (LSQ)
! polynom p: vektor počínající koeficientem členu s nejvyšší mocninou
! vstupní data: soubor 'data', výstupní data: soubor 'fit'
! překlad: gfortran -O3 -o testMypoly testMypoly.f90 -lopenblas

include 'mLinalg.f90'       ! modul s LAPACK solverem dgesv
include 'mypoly.f90'        ! modul s polyval a polyfit

program testMypoly
use mypoly
implicit none
real(8),allocatable :: x(:),y(:),p(:),xi(:),yi(:)
real(8) xmin,xmax
integer m,n,nn,nmax,f,ierr
character(:),allocatable :: f_in,f_out
m=1                         ! stupeň aproximačního polynomu
f_in='data'                 ! vstupní soubor s (x,y)
f_out='fit'                 ! výstupní soubor s (xi,yi)

! Načtení dat z f_in
open (newunit=f,file=f_in)  ! zjištění počtu řádků ve vstupním souboru
n=0; do; read (f,*,iostat=ierr); if (ierr/=0) exit; n=n+1; enddo
close (f)
n=n-1                       ! stupeň interpolačního polynomu
allocate (x(0:n),y(0:n))    ! jednorázová alokace
open (newunit=f,file=f_in)  ! načtení vstupních dat
do nn=0,n; read (f,*) x(nn),y(nn); enddo
close (f)
print '(3a,i0)','Soubor ',f_in,' - precteno radku: ',n+1

if (m==n) then
  p=polyfit(x,y)            ! polynomická interpolace
else
  p=polyfit(x,y,m)          ! polynomická regrese
endif

! Tabelování interpolačního polynomu
nmax=n*100; xmin=x(0)-1; xmax=x(n)+1
xi=[(xmin+(xmax-xmin)*nn/nmax,nn=0,nmax)]
yi=polyval(p,xi)

! Zápis dat do f_out (pro gnuplot)
open (newunit=f,file=f_out)
do nn=0,nmax; write (f,'(2(f10.6,x))') xi(nn+1),yi(nn+1); enddo
close (f)
print '(3a,i0)','Soubor ',f_out,' - zapsano radku: ',nmax+1
end program

Následuje fortranský modul mypoly, řádně využívající array language Fortranu. Najdeme v něm dvě funkce polyval(p,x), jednu pro skalární, druhou pro vektorový argument x. Polynomickou interpolaci realizuje funkce polyfit2(x,y) a regresi polynomem stupně deg funkce polyfit(x,y,deg).

! Modul s naší implementací polyval a polyfit
module mypoly
use mLinalg                ! modul s LAPACK solverem dgesv
implicit none
integer,parameter :: wp=8  ! working precision: 8 bytes
interface polyval          ! přetížení funkce polyval
  module procedure polyval,polyval_1d
end interface

contains

! Vyčíslení polynomu p[0]*x**n + ... + p[n], skalární x
function polyval(p,x) result (result)
real(wp),intent(in) :: p(0:),x
real(wp) result
integer n
result=p(0)
do n=1,ubound(p,1); result=result*x+p(n); enddo
end function

! Vyčíslení polynomu p[0]*x**n + ... + p[n], vektorové x
function polyval_1d(p,x) result (result)
real(wp),intent(in) :: p(0:),x(:)
real(wp) result(size(x))
integer nn
do nn=1,size(result); result(nn)=polyval(p,x(nn)); enddo
end function

! Interpolační polynom p[0]*x**n + ... + p[n]
function polyfit2(x,y) result (result)
real(wp),intent(in) :: x(0:),y(0:)
real(wp),allocatable :: result(:),V(:,:),c(:)
integer n,j
n=size(x)-1
allocate (V(0:n,0:n))      ! Vandermondova matice
V(:,0)=1
do j=1,n; V(:,j)=V(:,j-1)*x; enddo
c=V.gesv.y                 ! řešič algebraické soustavy z mLinalg
result=c(size(c):1:-1)     ! zrcadlení vektoru
end function

! Regresní polynom p[0]*x**m + ... + p[m], m=deg
function polyfit(x,y,deg) result (result)
real(wp),intent(in) :: x(0:),y(0:)
integer,intent(in),optional :: deg
real(wp),allocatable :: result(:),G(:,:),c(:),xx(:),sumx(:),sumy(:)
integer m,n,i,j
if (.not.present(deg)) then
  result=polyfit2(x,y)     ! interpolační polynom
  return
endif
m=deg; n=size(x)-1
allocate (xx(0:n),sumx(0:m*2),sumy(0:m),G(0:m,0:m))
xx=1
do i=0,m*2
  sumx(i)=sum(xx)
  if (i<=m) sumy(i)=sum(y*xx)
  xx=xx*x
enddo
do i=0,m; do j=0,m; G(i,j)=sumx(i+j); enddo; enddo  ! Gramova matice
c=G.gesv.sumy              ! řešič algebraické soustavy z mLinalg
result=c(size(c):1:-1)     ! zrcadlení vektoru
end function

end module

DÚ: Kreslení regresních polynomů

Domácím úkolem je spočítat pro \(n+1\) (třeba \(n=3\)) Vámi zvolených bodů koeficienty aproximačních polynomů stupňů \(0,1,\dots,n\), polynomy vyčíslit na husté síti a skreslit je všechny spolu se vstupními body do jednoho obrázku, podobného tomuto:

Tento obrázek těží ze schopnosti gnuplotu (příkaz fit) počítat metodou nejmenších čtverců parametry předepsané funkce aproximující data v souboru a klidně tak odstavit (nebo lépe: ověřit) naše snahy v Pythonu a Fortranu:

# gnuplot: polynomická aproximace (interpolace a regrese)
f0(x)=a0                          # konstantní funkce
f1(x)=b1*x+b0                     # polynom 1. stupně
f2(x)=(c2*x+c1)*x+c0              # polynom 2. stupně
f3(x)=((d3*x+d2)*x+d1)*x+d0       # polynom 3. stupně
fit f0(x) 'data' via a0           # metoda nejmenších čtverců
fit f1(x) 'data' via b1,b0
fit f2(x) 'data' via c2,c1,c0
fit f3(x) 'data' via d3,d2,d1,d0
set key top left
set samples 30
plot [0:5][-1:5] 'data' with points linestyle 7 pointsize 3,\
     'fit0' with lines linestyle 1,f0(x) with points pointtype 7 linecolor 1,\
     'fit1' w l ls 2,f1(x) w p pt 7 lc 2,\
     'fit2' w l ls 3,f2(x) w p pt 7 lc 3,\
     'fit3' w l ls 4,f3(x) w p pt 7 lc 4
set terminal pngc enhanced size 1024,768
set output 'fit.png'
replot
set output
unset terminal

Implementace (nelineární) metody nejmenších čtverců v gnuplotu je postavena na obecnějším algoritmu, než jsme předvedli zde, aproximační funkce nemusí být jen polynomem.

Zdrojové soubory .py, .m a .f90 s dnešními ukázkami: zip.