Cvičení č. 8: Lineární algebra

Obrátíme se k algoritmům lineární algebry (LA). Ve světě praktických počtů se řešiče algebraických rovnic s plnou maticí berou z knihovny LAPACK (LA Package) s podporou knihovny BLAS (Basic LA Subprograms). V Pythonu LAPACKovské řešiče poskytuje NumPy (SciPy také, a více), ve Fortranu připojíme příslušný knihovní podprogram ručně. Vyřešíme tedy soustavu lineárních algebraických rovnic a pro ověření pak vynásobíme matici soustavy s vektorem řešení – výsledkem by měla být původní pravá strana. Obrázky dnes nečekejme, ačkoliv… i matice mohou mít svou galerii.

Python: numpy.linalg

Jsme u kapitoly č. 1 učebnic numerické matematiky. Maticové násobení,

\[c_i=\sum_k A_{ik}b_k,\quad C_{ij}=\sum_k A_{ik}B_{kj}\,,\]

se v nich sice nutně neobjevuje, to je (samozřejmě jen na první pohled) příliš jednoduché, ale řešení soustav lineárních algebraických rovnic ano:

\[\sum_j A_{ij}x_j=b_i,\quad \sum_j A_{ij}X_{jm}=B_{im}\,.\]

My si pro obě operace najdeme pythonské funkce: pro maticové násobení numpy.matmul nebo ekvivalentní operátor @, pro řešení soustav numpy.linalg.solve nebo vybavenější scipy.linalg.solve. Řešiče se výslovně hlásí k LAPACKovským podprogramům, např. dgesv pro obecné plné matice (double-precision general-matrix solver). NumPy si svůj LAPACK a podpůrný BLAS nese s sebou, takže nic dalšího nepotřebujeme. Nachystáme čtvercovou matici A o nmax řádcích (1 na diagonále, 1/nmax mimo ni), do vektoru pravé strany b vložíme jedničky a vypočteme vektor řešení x. Jako zkoušku pak matici A vektorem x vynásobíme a měli bychom dostat původní pravou stranu, tedy samé jedničky. Zvolíme-li nmax=1000, výsledek obdržíme do sekundy, zvolíme-li 10_000, bude to za několik sekund. Řešení soustav s plnou maticí má časovou složitost kubickou, další zdvojnásobení nmax zosminásobí dobu běhu. Tedy: soudobá PC vracejí řešení těchto algebraických soustav v minutových časech pro matice o nízkých desetitisících řádků. Ve Windows v Task Manageru (Ctrl+Shift+Esc) si na kartě Performance můžeme povšimnout, že úloha zatíží 100 % CPU, tedy že běží paralelně na všech jádrech procesoru. I to umí BLAS a LAPACK.

# Python a NumPy: řešení soustavy lineárních algebraických rovnic a následná zkouška
import numpy
format='%s%20.15f%20.15f%20.15f'
nmax=10_000
A=numpy.full((nmax,nmax),1/nmax); numpy.fill_diagonal(A,1) # matice soustavy
b=numpy.ones(nmax)                                         # vektor pravé strany
print(format%('b =',b[0],b[1],b[2]))
x=numpy.linalg.solve(A,b)                                  # vektor řešení
print(format%('x =',x[0],x[1],x[2]))
b=A@x            # nebo numpy.matmul(A,x), lze i numpy.dot(A,x)
print(format%('b =',b[0],b[1],b[2]))
print(format%('d =',sum(abs(b-1)),numpy.sqrt(sum((b-1)**2)),max(abs(b-1))))  # normy rozdílu b

LAPACK obsahuje algoritmy pro řešení algebraických soustav, vlastních čísel, singulárních rozkladů ad., nízkoúrovňový BLAS provádí základní operace s vektory a maticemi. Je to právě BLAS, kdo dává LAPACKU a NumPy jejich rychlost. Ovšem BLAS není jen jeden: rozhraní jeho procedur (ve Fortranu a C/C++) je sice dáno pevně, ale úspěšnou implementaci provedlo hned několik týmů. Dobré jméno mají OpenBLAS (odvozený od GotoBLASu autora Kazushige Goto), ATLAS (Automatically Tuned …) a BLAS z Intel MKL (Math Kernel Library). Který BLAS je v akci v našem NumPy, to se dozvíme dotazem: import numpy; print(numpy.__config__.show()). V NumPy z pipu to patrně bude OpenBLAS a v Anacondě a Intel Pythonu BLAS z Intel MKL, ve Windows i v Linuxu. Rychlost a přesnost mají oba podobnou, v detailech se liší.

Fortran: BLAS a LAPACK

V klasických programovacích jazycích si pro lineární algebru programátor připojuje knihovny ručně. Je tedy třeba si BLAS a LAPACK opatřit. V Linuxu je to snadné, v repozitářích se najde OpenBLAS i ATLAS, obvyklou cestou se nainstalují tak, že je systém najde, a stačí pak přidat na překladový řádek volby -lblas nebo -llapack nebo obě. Jsou-li instalovány Intel oneAPI Base & HPC Toolkits, s Intel překladači se překládá s volbou -qmkl (Linux), resp. -Qmkl (Windows). Najít dobrý BLAS pro Windows nebývalo snadné, ale něco se změnilo a openblas.net nabízí Windows implementaci OpenBLAS-*-x64.zip ke stažení z github.com. Rozbalíme-li obsah archivu např. do adresáře c:\openblas, můžeme pomocí systempropertiesadvanced přidat do proměnné prostředí PATH adresář c:\openblas\bin a do proměnné LIBRARY_PATH uložit c:\openblas\lib a překládat příkazem gfortran -O3 file.f90 -lopenblas (nebo lze všechny soubory z adresářů bin a lib přikopírovat k našemu zdrojovému souboru a překládat s volbami -L. -lopenblas).

Pokud chceme dosáhnout stejné úrovně pohodlí jako v NumPy, musíme knihovní funkce dgesv (general-matrix solver) a dgemv (matrix-vector multiplication) s nechutně dlouhým rozhraním zabalit do funkcí (wrappers) w_dgesv a w_dgemv a ty raději ještě sdružit s operátory .gesv. a .gemv.. Pro maticové násobení (zahrnující i součin matice a vektoru) má fortranská array language funkci matmul, tu můžeme použít také.

! Fortran, BLAS a LAPACK: řešení soustavy lineárních algebraických rovnic a následná zkouška
! Překlad ve Windows: gfortran -O3 file.f90 -lopenblas; ifort -O3 file.f90 -Qmkl
module mLinalg
implicit none
interface operator (.gesv.); module procedure w_dgesv; end interface
interface operator (.gemv.); module procedure w_dgemv; end interface
contains
function w_dgesv(A,b) result (x)                     ! wrapper pro dgesv
   real(8),intent(in) :: A(:,:),b(:)
   real(8) x(size(b))
   real(8),allocatable :: AA(:,:)
   integer,allocatable :: ipiv(:)
   integer n,info
   n=size(b); AA=A; x=b; ipiv=b
   call dgesv(n,1,AA,size(AA,1),ipiv,x,size(x),info) ! explicitní LAPACK
end function
function w_dgemv(A,b) result (r)                     ! wrapper pro dgemv
   real(8),intent(in) :: A(:,:),b(:)
   real(8) r(size(b))
   integer n
   n=size(b)
   call dgemv('N',n,n,1._8,A,n,b,1,0._8,r,1)         ! explicitní BLAS
end function
end module

use mLinalg
real(8),allocatable :: A(:,:),b(:),x(:)
character(20) :: format='(a,3f20.15)'
nmax=10000
allocate (A(nmax,nmax),b(nmax))
A=1._8/nmax; do n=1,nmax; A(n,n)=1; enddo  ! matice soustavy
b=1                                        ! vektor pravé strany
print format,'b =',b(1:3)
x=A.gesv.b                                 ! vektor řešení
print format,'x =',x(1:3)
b=A.gemv.x      ! nebo standardní fortranské b=matmul(A,x)
print format,'b =',b(1:3)
print format,'d =',sum(abs(b-1)),sqrt(sum((b-1)**2)),maxval(abs(b-1))  ! normy rozdílu b
end program

Když se nám zdaří spustit Python i Fortran s OpenBLAS, uvidíme na řádku s výpisem prvních prvků vektoru x prakticky to samé, a stejné výsledky (ale jiné než OpenBLAS) by vrátil Python i Fortran s Intel MKL. Zkouška už to samé v Pythonu a Fortranu neukazuje, pythonský matmul operátor @ zřejmě nevolá BLASovské dgemv a výsledné b je proti dgemv trochu jiné. To by se ostatně stalo i ve Fortranu, kdybychom místo dgemv zavolali matmul, a podobný efekt jsme viděli i dříve při numerické integraci, kdy funkce sum a math.fsum také vracely různé výsledky. Na posledním výpisovém řádku obou ukázek jsou proto spočteny tři normy (L1, L2, L∞) rozdílu zkouškové a originální pravé strany a můžeme tak exaktněji porovnat, jak by pro kterou matici dopadla která kombinace řešičů a zkoušek.

Matlab/Octave: * a \

Vraťme se dnes i k Matlabu a Octavu, kde je lineární algebra nejvíc doma. Maticové násobení je zde reprezentováno operátorem * a řešení algebraických soustav nabízí operátor \. Skript ekvivalentní našim předchozím ukázkám je úžasně úsporný:

% Matlab/Octave: řešení soustavy lineárních algebraických rovnic a následná zkouška
format='%s%20.15f%20.15f%20.15f\n';
nmax=10000;
A=ones(nmax)/nmax; A=A-diag(diag(A))+eye(nmax);  % matice soustavy
b=ones(nmax,1);                                  % sloupcový vektor pravé strany
fprintf(format,'b =',b(1:3))
x=A\b;                                           % vektor řešení
fprintf(format,'x =',x(1:3))
b=A*x;                                           % zkouška
fprintf(format,'b =',b(1:3))
fprintf(format,'d =',sum(abs(b-1)),sqrt(sum((b-1).^2)),max(abs(b-1)))  % normy
fprintf(format,'d =',norm(b-1,1),norm(b-1,2),norm(b-1,Inf))            % totéž

DÚ: Inverzní Hilbertovy matice

Domácí úkol: úpravou některé z předchozích ukázek připravte program, který pro čtvercové Hilbertovy matice \(H_{ij}=1/(i+j-1),\) \(i,j=1,\dots,N,\) spočítá inverzní matice a vypíše součet všech prvků každé z nich, \(S(N)=\sum_{i,j=1}^N (H^{-1})_{ij}\). To vše podstupte pro počet řádků \(N\) od 1 do 12 anebo až do 20 a ze získaných výsledků uhodněte, jak lze \(S(N)\) vyjádřit pomocí \(N\) úsporněji. Inverzní matici \(A^{-1}\) lze získat řešením soustavy \(AX=I\), kde \(I\) je jednotková matice. Funkce numpy.linalg.solve() řeší i tuto úlohu, tedy s maticí na pravé straně. Domácí úkol zpracujte právě takto. Dobrovolně si můžete ověřit, co vracejí funkce numpy.linalg.inv(), scipy.linalg.invhilbert(n) a scipy.linalg.invhilbert(n,exact=True).

Z toho, co jsme výše viděli, lze soudit, že LAPACK svižně invertuje plné matice s desetitisíci řádků. Ne ovšem každou. Pro Hilbertovy matice numerická inverze v 8bytové reálné aritmetice selhává už u desítky řádků. Jak tušíme, je-li třeba Hilbertovy matice invertovat, chodí se jinými cestami než metodami LAPACKu.