$n$ tundmatuga lineaarne võrrandisüsteem on üldkujul järgmine: \begin{align*} a_{11}x_1+a_{12}x_2+\cdots+a_{1n}x_n &= b_1\\ a_{21}x_1+a_{22}x_2+\cdots+a_{2n}x_n &= b_2\\ &\vdots\\ a_{m1}x_1+a_{m2}x_2+\cdots+a_{mn}x_n &= b_m\\ \end{align*} Teada on kordajad $a_{ij}$ ja vabaliikmed $b_j$ ning tuleb avaldada tundmatud $x_i$. Ilmselt on kolm varianti. Kui võrrandeid on samapalju nagu tundmatuid ($m=n$) ning kõik võrrandid on sõltumatud (st ükski võrrand pole teiste võrrandite lineaarkombinatsioon), siis lahend eksisteerib ja on ühene. Kui võrrandeid on rohkem kui tundmatuid, siis on süsteem ülemääratud. Sellegipoolest võib mõttekas olla sellise parima lahendi leidmine, mis minimeerib võrrandite vasakute ja paremate poolte erinevuse (näiteks vähimruutude mõttes). Kui tundmatuid on rohkem kui võrrandeid, siis süsteemi lõpuni lahendada ei saa (parimal juhul saab lahendi avaldada sellisel kujul, kus mõned tundmatud jäävad parameetritena sisse).
Kordajatest $a_{ij}$ võime konstrueerida omaette tabeli täpselt samas järjestuses nagu need võrrandisüsteemis esinevad: $$\vec A=\begin{bmatrix}a_{11} & \cdots & a_{1n}\\ \vdots & \ddots & \vdots\\ a_{m1} & \cdots & a_{mn}\end{bmatrix}.$$ Sellise matemaatilise objekti nimi on maatriks. Analoogiliselt saab konstrueerida üheveerulised maatriksid $$\vec x=\begin{bmatrix}x_1\\ \vdots\\ x_n\end{bmatrix}, \quad \vec b=\begin{bmatrix}b_1\\ \vdots\\ b_m\end{bmatrix}.$$ Lineaarvõrrandisüsteemis esineva arvutusmustri realiseerib maatrikskorrutis, nii et võrrandisüsteemi saab maatrikskujul kirja panna hästi kompaktselt: $\vec A\vec x=\vec b$.
Oletagem nüüd, et $\vec A$ on $n\times n$ ruutmaatriks ja püüame ette kujutada, kuidas süstemaatiliselt lahendada sellist võrrandisüsteemi. Selleks tuleb asuda tundmatuid järjekindlalt elimineerima. Ilmselt eksisteerib kolm elementaaroperatsiooni, mis ei mõjuta võrrandisüsteemi lahendit: (1) kaks võrrandit tohib ära vahetada, (2) iga võrrandi tohib läbi korrutada nullist erineva teguriga ning (3) igale võrrandile tohib otsa liita mõne teise võrrandi (mille võib enne veel läbi korrutada mingi teguriga). Pärast iga sellist tegevust kõik tundmatud $x_i$ jäävad ikka oma positsioonidele, seega neid elementaaroperatsioone on parem ette kujutada arvutabelis, kus algselt $\vec A$ ja $\vec b$ on asetatud üksteise kõrvale: $$\begin{bmatrix}a_{11} & \cdots & a_{1n} & b_1\\ \vdots & \ddots & \vdots & \vdots\\ a_{n1} & \cdots & a_{nn} & b_n\end{bmatrix}$$ Nimetatud elementaaroperatsioonid tähendavad nüüd lihtsalt vastavaid tegevusi selle maatriksi ridadega. Süstemaatiliselt muutujaid elimineerides (st maatriksisse nulle tekitades) teisendame selle maatriksi viimaks kujule: $$\begin{bmatrix}1 & 0 & \cdots & 0 & b_1'\\0 & 1 & \cdots & 0 & b_2'\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & b_n'\end{bmatrix}$$ Seega olemegi saanud kätte lahendi: $x_1=b_1'$, $x_2=b_2'$, …. Sellist algoritmi nimetatakse Gauss-Jordani elimineerimismeetodiks. Selle keerukus (aritmeetiliste tehete arv) on $\mathcal{O}(n^3)$.
On ilmne, et kõik saadud $b_i'$ väärtused kujutavad endast lineaarkombinatsioone erinevatest $b_i$ väärtustest. Maatrikskorrutis tekitabki sellised lineaarkombinatsioonid. Seega formaalselt võib võrrandi lahendi kirja panna kujul $\vec x=\vec{A}^{-1}\vec b$, kus $\vec{A}^{-1}$ on samuti teatav $n\times n$ maatriks. Maatriksite korrutamine (nagu ka arvude korrutamine) on assotsiatiivne, millest järeldub $\vec{A}\vec{A}^{-1}=\vec{A}^{-1}\vec{A}=\vec{I}$, kus ühikmaatriks $$\vec I=\begin{bmatrix}1 & 0 & \cdots & 0 \\0 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 \end{bmatrix}.$$ Nende omaduste tõttu maatriksit $\vec{A}^{-1}$ nimetatakse maatriksi $\vec A$ pöördmaatriksiks. Gauss-Jordani elimineerimismeetodit veidi laiendades suudaksime ilmselt ka pöördmaatriksi välja arvutada. Formaalselt on pöördmaatriksi kasutamine mugav, aga võrrandisüsteemi praktiline lahendamine pöördmaatriksi väljaarvutamist ei eelda (see oleks arvutuslikult ebaefektiivsem ja ebatäpsem).
Iga vabaliikmete vektori $\vec{b}$ puhul tuleb lahendi saamiseks kirjeldatud algoritm uuesti läbi teha (isegi kui maatriks $\vec{A}$ jääb samaks). Osutub, et maatriks $\vec{A}$ on võimalik lahutada kahe maatriksi korrutiseks: $\vec{A} = \vec{L}\vec{U}$, kus $\vec{L}$ (lower) on kolmnurkne maatriks nullidega pealpool peadiagonaali ning $\vec{U}$ (upper) on kolmnurkne maatriks nullidega allpool peadiagonaali. Seda nimetatakse LU-dekompositsiooniks. Kui kordajate maatriksi LU-dekompositsioon on korra ära tehtud, siis süsteemi lahendamine taandub kahe järjestikuse kolmnurkse süsteemi lahendamisele: esmalt tuleb lahendada süsteem $\vec{L}\vec{y}=\vec b$ ja seejärel süsteem $\vec{U}\vec{x}=\vec y$. Nende lahendamine on vaid keerukusega $\mathcal{O}(n^2)$.
On olemas kaks sarnast teeki lineaaralgebra vahenditega: numpy.linalg ja scipy.linalg. Viimane sisaldab veidi põhjalikumat komplekti vahendeid. Mõningaid piisavalt lihtsaid tegevusi (näiteks maatriksite liitmine või korrutamine) saab teostada ka ilma nende teekideta, kasutades vaid NumPy massiivi meetodeid.
Võrrandisüsteemi lahendamise LU-dekompositsiooniga realiseerib funktsioon numpy.linalg.solve. On olemas ka scipy.linalg.solve täiendavate optimeerimisvõimalustega. Näitena lahendame süsteemi \begin{align*}3x_1-x_2-2x_3 &= 1\\ x_1+4x_2+2x_3 &= 4\\ -5x_1+2x_2+6x_3 &= 6\end{align*}.
import numpy as np
A = np.array( ((3, -1, -2), (1, 4, 2), (-5, 2, 6)), dtype=np.float )
b = np.array( (1, 4, 6), dtype=np.float )
x = np.linalg.solve(A, b)
print('x =', x)
Nagu mainitud, saab lahendi avaldada põhimõtteliselt ka pöördmaatriksi kaudu: $\vec{x}=\vec{A}^{-1}\vec{b}$. Pöördmaatriksi annab funktsioon numpy.linalg.inv ning maatrikskorrutise saab leida numpy massiivi meetodiga dot. Seejuures on tõlgendamise küsimus, kas ühemõõtmelise massiivi b puhul on tegemist rea- või veeruvektoriga. A.dot(b) korral tõlgendatakse b veeruvektorina, b.dot(A) puhul reavektorina.
x = np.linalg.inv(A).dot(b)
print('x =', x)
NumPy massiivide otsene korrutamine (operaatoriga *) annab hoopis teistsuguse tulemuse (element-haaval korrutise, kuhu on vajadusel kaasatud broadcasting-mehhanism). Kui on siiski soov opereerida massiividega A ja b nii, nagu oleks tegemist maatriksi ja vektoriga lineaaralgebra mõistes, peaks mõlemad massiivid defineerima kahemõõtmelisena ja numpy.array asemel käsuga numpy.matrix. Nagu järgnevast näha, astendamine (täisarvulise astmenäitajaga) on maatriksite jaoks spetsiifilise tähendusega ja kooskõlas pöördmaatriksi kontseptsiooniga.
A = np.matrix(A, copy=False)
b = np.matrix(b, copy=False).T
x = A**(-1) * b
print('x =\n', x)
Eksisteerib lineaaralgebra probleeme, kus enamik elemente maatriksis on nullid. Arvutusressursi säästmise huvides pole mõttekas ei nende nullide säilitamine arvuti mälus ega ka nendega aritmeetiliste operatsioonide tegemine. Seega on mõistlik kasutada spetsiaalseid andmestruktuure ja algoritme, mis on olemas teegis scipy.sparse. Näiteks lil_matrix (linked list) võib olla sobilik hõreda maatriksi konstrueerimiseks (ükshaaval elementide lisamiseks) samas kui csc_matrix (compressed sparse column) või csr_matrix (compressed sparse row) on sobivam arvutamiseks. Mõningad funktsioonid (diags, identity, bmat, jt) abistavad sobiva struktuuriga hõredate maatriksite loomisel.
Konstrueerime võrdlemisi suure lineaarse võrrandisüsteemi kordajate maatriksi, kus vaid peadiagonaali lähedal on nullist erinevad elemendid:
from scipy.sparse import diags
n = 500
A = diags( (-1, 2, -1), (-1, 0, 1), shape=(n, n), format='csc')
b = np.arange(n)
A[:7,:7].toarray()
Nüüd võime veenduda, et scipy.sparse.linalg.spsolve lahendab selle süsteemi kiiremini kui numpy.linalg.solve:
from scipy.sparse.linalg import spsolve
%timeit x = spsolve(A, b)
A = A.toarray()
%timeit x = np.linalg.solve(A, b)