Linearni algebra

<p Python je efektivni, ucinny a relativne snadny nastroj na provadeni rozlicnych algebraickych operaci. Klicovym objektem v numericke knihovne numpy je numpy.array, coz je generalizovana matice. </p>

``It is not knowledge, but the act of learning, not possession but the act of getting there, which grants the the greatest enjoyment."

                                                                             Carl Friedrich Gauss

Definice matice

Definujme si nejprve nejakou obecnou realnou matici $\pmb{A}$ s prvky z $\mathbb{R}$


In [1]:
import numpy as np

A = np.array([[3,-5,8],
             [-1,2,3],
             [-5,-6,2]], dtype=np.float64)

A


Out[1]:
array([[ 3., -5.,  8.],
       [-1.,  2.,  3.],
       [-5., -6.,  2.]])

Vsimnete si, ze si definuji i datovy typ promennych v matici (np.float64) viz


In [2]:
A.dtype


Out[2]:
dtype('float64')

Duvodem je "programovaci hygiena" a snaha, aby si python nevylozil datovy typ promenne jinak, nez bychom zrovna chteli. Srovnejte napr. s timto


In [3]:
divne = np.array([[3,-5,8],
             [-1,2,'co?'],
             [-5,-6,2]])
divne


Out[3]:
array([['3', '-5', '8'],
       ['-1', '2', 'co?'],
       ['-5', '-6', '2']], 
      dtype='<U21')

Vratme se k matici $\pmb{A}$ podivejme se, co objekt nabizi. Nejprve jeho rozmery


In [4]:
np.shape(A)


Out[4]:
(3, 3)

a celkovy pocet jejich prvku


In [5]:
np.size(A)


Out[5]:
9

Pristup k jednotlivym prvkum matice je pak


In [6]:
print(A[0, 1])
print(A[1, 0])


-5.0
-1.0

Radky matice $\pmb{A}$ zase ziskame jako


In [7]:
print('A[0]', A[0], A[0, :])
print('A[1]', A[1])
print('A[2]', A[2])


A[0] [ 3. -5.  8.] [ 3. -5.  8.]
A[1] [-1.  2.  3.]
A[2] [-5. -6.  2.]

Zajimat by nas dale mohla treba stopa matice $\text{Tr}(\pmb{A})$


In [8]:
np.trace(A)


Out[8]:
7.0

In [9]:
np.trace?

ci stopa po libovolne jine vedlejsi diagonaly


In [10]:
np.trace(A, 0)


Out[10]:
7.0

Diagonaly se pocitaji od hlavni, ktera ma oznaceni 0


In [11]:
print(np.trace(A, 0))
print(np.trace(A, 0)==np.trace(A))


7.0
True

Stopu pod diagonalou ziskame jako


In [12]:
np.trace(A, -1)


Out[12]:
-7.0

hodnoty mimo matici se vrati jako 0.0


In [13]:
print(np.trace(A, 3))

print(np.trace(A, -4))


0.0
0.0

Podobne lze extrahovat primo jednotlive diagonaly


In [14]:
print(np.diag(A))
print(np.diag(A, -1))
print(np.diag(A, -2))
print(np.diag(A, -3))


[ 3.  2.  2.]
[-1. -6.]
[-5.]
[]

Dale lze snadno provest transpozici matice $\pmb{A}^T$ jako


In [15]:
np.transpose(A)


Out[15]:
array([[ 3., -1., -5.],
       [-5.,  2., -6.],
       [ 8.,  3.,  2.]])

nebo jeste uspornejsim zapisem


In [16]:
A.T


Out[16]:
array([[ 3., -1., -5.],
       [-5.,  2., -6.],
       [ 8.,  3.,  2.]])

Operace s maticemi

Prejdeme nyni k maticovym operacim, ktere jsou soucasti balicku numpy.linalg. Zacneme inverzi matice $\pmb{A}^{-1}$


In [17]:
np.linalg.inv(A)


Out[17]:
array([[ 0.08494208, -0.14671815, -0.11969112],
       [-0.05019305,  0.17760618, -0.06563707],
       [ 0.06177606,  0.16602317,  0.003861  ]])

nebo nalezenim jejich vlastnich hodnot a vlastnich vektoru $\pmb{A} \pmb{v} = \lambda \pmb{v}$


In [18]:
np.linalg.eig(A)


Out[18]:
(array([ 1.25656268+7.49299485j,  1.25656268-7.49299485j,  4.48687464+0.j        ]),
 array([[-0.68307576+0.j        , -0.68307576-0.j        ,  0.79215535+0.j        ],
        [-0.27522781-0.05458384j, -0.27522781+0.05458384j, -0.57277576+0.j        ],
        [-0.02315491-0.6739003j , -0.02315491+0.6739003j , -0.21075538+0.j        ]]))

Tento objekt ma dve casti


In [19]:
vl_hodnoty, vl_vektory = np.linalg.eig(A)

vlastni hodnoty $(\lambda_1, \lambda_2, \lambda_3)$


In [20]:
print(vl_hodnoty)


[ 1.25656268+7.49299485j  1.25656268-7.49299485j  4.48687464+0.j        ]

a vlastni vektory $\pmb{v}_1, \pmb{v}_2, \pmb{v}_3$


In [21]:
vl_vektory


Out[21]:
array([[-0.68307576+0.j        , -0.68307576-0.j        ,  0.79215535+0.j        ],
       [-0.27522781-0.05458384j, -0.27522781+0.05458384j, -0.57277576+0.j        ],
       [-0.02315491-0.6739003j , -0.02315491+0.6739003j , -0.21075538+0.j        ]])

Vlastni hodnoty a vektory teto matice jsou ale komplexni $\pmb{v} = \pmb{x} + i \pmb{y}$, komplexni jednotka je v pythonu oznacovana jako '1j'.


In [22]:
print(vl_hodnoty.dtype)
print(vl_vektory.dtype)


complex128
complex128

oba vysledne objekty jsou ale opet typu numpy.array a da se s nim tak i dale zachazet


In [23]:
print(type(vl_hodnoty))
print(type(vl_vektory))


<class 'numpy.ndarray'>
<class 'numpy.ndarray'>

numpy.array rovnez podporuje nasobeni, a to jak nasobeni matice skalarem $2 \pmb{A}$


In [24]:
A*2


Out[24]:
array([[  6., -10.,  16.],
       [ -2.,   4.,   6.],
       [-10., -12.,   4.]])

umocnovani $\pmb{A}^2$, logaritmovani $\ln(\pmb{A})$, odmocnovani $\pmb{A}^{1/2}$, atd. vsech prvku vektoru


In [25]:
print('umocneni na druhou:')
print(A**2)
print('prirozeny logaritmus:')
print(np.log(A))
print('odmocneni:')
print(np.sqrt(A))


umocneni na druhou:
[[  9.  25.  64.]
 [  1.   4.   9.]
 [ 25.  36.   4.]]
prirozeny logaritmus:
[[ 1.09861229         nan  2.07944154]
 [        nan  0.69314718  1.09861229]
 [        nan         nan  0.69314718]]
odmocneni:
[[ 1.73205081         nan  2.82842712]
 [        nan  1.41421356  1.73205081]
 [        nan         nan  1.41421356]]
/home/ziky/.virtualenvs/plasmasolve3/lib/python3.4/site-packages/ipykernel/__main__.py:4: RuntimeWarning: invalid value encountered in log
/home/ziky/.virtualenvs/plasmasolve3/lib/python3.4/site-packages/ipykernel/__main__.py:6: RuntimeWarning: invalid value encountered in sqrt

Pri logaritmovani ci odmocneni jsme ale pro nektere hodnoty dostali 'nan', coz je chybova hlaska (nan=Not A Number). Vskutku, odmocnina ci logaritmus ze zaporneho realneho cisla neni v oboru realnych cisel definovana (opet 1:0 pro python :-)). Pokud chceme ale ziskat reseni odmocniny v oboru komplexnich cisel, je treba pretypovat $\pmb{A}$ na komplexni cisla


In [26]:
print('odmocneni s resenim v oboru komplexnich cisel:')
print(np.sqrt(np.complex128(A)))


odmocneni s resenim v oboru komplexnich cisel:
[[ 1.73205081+0.j          0.00000000+2.23606798j  2.82842712+0.j        ]
 [ 0.00000000+1.j          1.41421356+0.j          1.73205081+0.j        ]
 [ 0.00000000+2.23606798j  0.00000000+2.44948974j  1.41421356+0.j        ]]

Skalarni a vektorovy soucin

Matice numpy.array lze rovnez mezi sebou nasobit. Definujme si proto nyni jednorozmerny vektor $\pmb{b}$ slozeny obecne z komplexnich cisel


In [27]:
b = np.array([4, 3, -1j])
b


Out[27]:
array([ 4.+0.j,  3.+0.j, -0.-1.j])

Tento vektor muzeme vynasobit s matici $\pmb{A}$ jednak prvek po prvku $A_{ij} b_j, j = 1,3$ viz


In [28]:
A


Out[28]:
array([[ 3., -5.,  8.],
       [-1.,  2.,  3.],
       [-5., -6.,  2.]])

In [29]:
A*b


Out[29]:
array([[ 12.+0.j, -15.+0.j,   0.-8.j],
       [ -4.+0.j,   6.+0.j,   0.-3.j],
       [-20.+0.j, -18.+0.j,   0.-2.j]])

tedy nasobeni $\pmb{A} \otimes \pmb{b}$ nyni probiha takto


In [30]:
print('A[0]*b:', A[0]*b)
print('A[1]*b:', A[1]*b)
print('A[2]*b:', A[2]*b)


A[0]*b: [ 12.+0.j -15.+0.j   0.-8.j]
A[1]*b: [-4.+0.j  6.+0.j  0.-3.j]
A[2]*b: [-20.+0.j -18.+0.j   0.-2.j]

Nyni ukrok stranou, definujme vektor $\pmb{c}$


In [31]:
c = np.array([1,2,3])

a provedme skalarni soucin s $\pmb{b} \cdot \pmb{c}$


In [32]:
b.dot(c)


Out[32]:
(10-3j)

Musi pochopitelne platit $\pmb{b} \cdot \pmb{c} = \pmb{c} \cdot \pmb{b}$


In [33]:
b.dot(c) == c.dot(b)


Out[33]:
True

Procedura .dot dale umi provadet nasobeni matice a vektoru. Vynasobme takto tedy $\pmb{A} \cdot \pmb{b}$


In [34]:
A.dot(b)


Out[34]:
array([ -3.-8.j,   2.-3.j, -38.-2.j])

nyni ovsem neplati $\pmb{A} \cdot \pmb{b} = \pmb{b} \cdot \pmb{A}$


In [35]:
(A.dot(b) == b.dot(A))**2


Out[35]:
array([False, False, False], dtype=bool)

In [36]:
(A.dot(b) == b.dot(A))==True


Out[36]:
array([False, False, False], dtype=bool)

Vskutku, vektor $\pmb{b} \cdot \pmb{A}$ je jiny, jak se probira jiz v linearni algebre v prvaku. Python tedy opet pocita dobre. :-)


In [37]:
b.dot(A)


Out[37]:
array([  9.+5.j, -14.+6.j,  41.-2.j])

Definujeme nyni matici $\pmb{B}$ a provedeme maticove nasobeni:


In [38]:
B = np.array([[0,0,1],[0,2,0],[1,0,0]])

$\pmb{A}\times\pmb{B}$ a naopak


In [39]:
np.matmul(A,B)


Out[39]:
array([[  8., -10.,   3.],
       [  3.,   4.,  -1.],
       [  2., -12.,  -5.]])

In [40]:
np.matmul(B,A)


Out[40]:
array([[-5., -6.,  2.],
       [-2.,  4.,  6.],
       [ 3., -5.,  8.]])

Vyzkousejme si komutacni relaci $[\pmb{A},\pmb{B}]=\pmb{A}\times\pmb{B}-\pmb{B}\times\pmb{A}$ (muze se Vam hodit v kvantovce)


In [41]:
np.matmul(A,B) - np.matmul(B,A)


Out[41]:
array([[ 13.,  -4.,   1.],
       [  5.,   0.,  -7.],
       [ -1.,  -7., -13.]])

Zkusme definovat vektory $\pmb{a}, \pmb{b}$ a provedme jejich vektorovy soucin


In [42]:
a = np.array([0,0,1])
b = np.array([1,0,0])

In [43]:
c = np.cross(a,b)
c


Out[43]:
array([0, 1, 0])

A presvedcme se, ze je kolmy na oba puvodni vektory...


In [44]:
np.dot(a,c)


Out[44]:
0

In [45]:
np.dot(b,c)


Out[45]:
0

Reseni soustav linearnich rovnic

Podtrida linalg v tride numpy dale (mimo jine) obsahuje radu algoritmu na reseni soustav rovnic. Jako priklad pouzijme nyni vektor $\pmb{b}$ jako pravou stranu soustavy rovnic $\pmb{A} \pmb{x} = \pmb{b}$ a vyresme ji. Dostaneme vektor


In [46]:
np.linalg.solve(A,b)


Out[46]:
array([ 0.08494208, -0.05019305,  0.06177606])

Operace s prvky matic

Nyni se zminme o nekolika dalsich uzitecnych procedurach v numpy. Zacneme s tzv. slicing. Nekolik prikladu nasleduje, uhadnete, co provadeji


In [47]:
A


Out[47]:
array([[ 3., -5.,  8.],
       [-1.,  2.,  3.],
       [-5., -6.,  2.]])

In [48]:
A[:, 0]


Out[48]:
array([ 3., -1., -5.])

In [49]:
A[0, :]


Out[49]:
array([ 3., -5.,  8.])

In [50]:
A[1:, :]


Out[50]:
array([[-1.,  2.,  3.],
       [-5., -6.,  2.]])

In [51]:
A[0, :-2]


Out[51]:
array([ 3.])

Dalsi uzitecna procedura je rollaxis


In [52]:
np.rollaxis(A, 1, 0)


Out[52]:
array([[ 3., -1., -5.],
       [-5.,  2., -6.],
       [ 8.,  3.,  2.]])

Procedura zpusobila prehozeni radku za sloupce a naopak, coz muzeme videt pri srovnani s puvodni matici $\pmb{A}$ (de facto tak provedla transpozici)


In [53]:
A


Out[53]:
array([[ 3., -5.,  8.],
       [-1.,  2.,  3.],
       [-5., -6.,  2.]])

Procedura roll zase umoznuje posunout vsechny prvky o danou konstantu


In [54]:
np.roll(A, 1, axis=0)


Out[54]:
array([[-5., -6.,  2.],
       [ 3., -5.,  8.],
       [-1.,  2.,  3.]])

In [55]:
np.roll(A, 1, axis=1)


Out[55]:
array([[ 8.,  3., -5.],
       [ 3., -1.,  2.],
       [ 2., -5., -6.]])

In [56]:
np.roll(A, 1)


Out[56]:
array([[ 2.,  3., -5.],
       [ 8., -1.,  2.],
       [ 3., -5., -6.]])

In [57]:
np.roll(A, 4)


Out[57]:
array([[ 3., -5., -6.],
       [ 2.,  3., -5.],
       [ 8., -1.,  2.]])

Vyznamne matice

Sikovne jsou metody na automaticke generovani matic, ktere se skladaji ze stejnych prvku


In [58]:
np.zeros(np.shape(A))


Out[58]:
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])

In [59]:
np.ones(np.shape(A))


Out[59]:
array([[ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.]])

In [60]:
np.eye(3)


Out[60]:
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])

In [61]:
np.eye(3, k=1)


Out[61]:
array([[ 0.,  1.,  0.],
       [ 0.,  0.,  1.],
       [ 0.,  0.,  0.]])

In [62]:
(1e4+1j)*np.eye(3, k=-1)


Out[62]:
array([[     0.+0.j,      0.+0.j,      0.+0.j],
       [ 10000.+1.j,      0.+0.j,      0.+0.j],
       [     0.+0.j,  10000.+1.j,      0.+0.j]])

Cviceni

Najdete vlastni hodnoty a vektory nasledujici matice

$\pmb{M}= \begin{pmatrix} 1& -5& 3& 0& 0& 0& 0& 0\\ 3& 2& -3& 5& 0& 0& 0& 0\\ 0& -2& 1& 5& 9& 0& 0& 0\\ 0& 0& 9& -1& 5& 4& 0& 0\\ 0& 0& 0& 0& 2& -3& -2+5\mathrm{i}& 0\\ 0& 0& 0& 0& 2& 0& 1& -6\\ 0& 0& 0& 0& 0& -3& 2& 7\\ 0& 0& 0& 0& 0& 0& 9& 1\\ \end{pmatrix}$

a dale vyreste soustavu rovnic $\pmb{M}\cdot \pmb{x} = \pmb{v}_3$, kde $\pmb{v}_3$ je treti vlastni vektor matice $\pmb{M}$. Dale

  • vypiste matici $\pmb{M}$ slozenou z cisel nalezejicich do $\mathbb{C}$ jako soucet dvou matic s koeficienty v $\mathbb{R}$,
  • zmente komplexni matici v realnou a naopak a vypiste jako nove matice tak, nasledne je sectete a vysledek ulozte do matice $\pmb{D}$,
  • provedte komplexni sdruzeni matice $\pmb{M}$,
  • matici $\pmb{M}$ Hermitovsky sdruzte a provedte $\pmb{M}\times\pmb{M}^H$ a porovnejte s vhodne vytvorenou jednotkovou matici $\pmb{I}$,
  • provedte $|\pmb{M}|^2$ a $|\pmb{D}|^2$ a otestujte $|\pmb{M}|^2=|\pmb{D}|^2$ a $\pmb{M}=\pmb{D}$,
  • provedte rovnez $|\pmb{M}^2|$ a otestujte $|\pmb{M}|^2=|\pmb{M}^2|$.

Vyzkousejte python z kvantove mechaniky:

Uvazte nasledujici matice:

$\sigma_1=\begin{pmatrix} 0& 1\\ 1& 0\\ \end{pmatrix}, \sigma_2=\begin{pmatrix} 0& -\mathrm{i}\\ \mathrm{i}& 0\\ \end{pmatrix}, \sigma_3=\begin{pmatrix} 1& 0\\ 0& -1\\ \end{pmatrix},$

a overte komutacni relace

$[\pmb{\sigma}_1,\pmb{\sigma}_2]=2\mathrm{i}\pmb{\sigma}_3\\ [\pmb{\sigma_2},\pmb{\sigma}_3]=2\mathrm{i}\pmb{\sigma}_1\\ [\pmb{\sigma_3},\pmb{\sigma}_1]=2\mathrm{i}\pmb{\sigma}_2$

a antikomutacni relace $\{\pmb{C},\pmb{D}\}=\pmb{C}\times\pmb{D}+\pmb{D}\times\pmb{C}$

$\{\pmb{\sigma}_1,\pmb{\sigma}_1\}=2\pmb{I}\\ \{\pmb{\sigma}_1,\pmb{\sigma}_2\}=0$

Konecne naleznete vlastni hodnoty a vlastni vektory matic $\pmb{\sigma}_1$, $\pmb{\sigma}_2$ a $\pmb{\sigma}_3.\\\\\\\\\\\\$

``As far as the laws of mathematics refer to reality, they are not certain; and as far as they are certain, they do not refer to reality.

                                                                                  Albert Einstein

In [ ]: