6.13. Lineaarialgebra

Lineaarialgebra kameralla on pienten matriisien työtä: 3x3-rotaatio, joka yhdistää IMU-näytteen maailmankehykseen, kalibrointimatriisi, joka oikaisee linssin, Kalman-suotimen tila-kovarianssipäivitys, polynomisovitus, jonka normaaliyhtälöt tulevat ulos pikkuruisena lineaarisena ratkaisuna. numpy.linalg- ja scipy.linalg-alimoduulit kattavat juuri tuon mittakaavan.

Funktiot sijaitsevat kahdessa moduulissa. Matriisitulo-operaatiot ovat numpy:n ylätasolla; hajotelmat ja matriisin käänteismatriisi ovat numpy.linalg:n alla; varsinaiset lineaaristen järjestelmien ratkaisijat ovat scipy.linalg:n alla. Esimerkiksi 2x2-matriisikertolasku:

from ulab import numpy as np

A = np.array([[1, 2], [3, 4]], dtype=np.float)
B = np.array([[5, 6], [7, 8]], dtype=np.float)
np.dot(A, B)
# array([[19.0, 22.0],
#        [43.0, 50.0]])

6.13.1. Mitä on saatavilla

  • dot() – matriisi- tai vektoritulo.

  • cross() – 3-ulotteisten vektorien ristitulo.

  • trace() – lävistäjän summa.

  • inv() – matriisin käänteismatriisi.

  • det() – determinantti.

  • cholesky() – Choleskyn hajotelma (symmetrinen positiivisesti definiitti syöte).

  • eig() – reaalisen symmetrisen matriisin ominaisarvot ja ominaisvektorit.

  • norm() – vektorin tai matriisin 2-normi.

  • qr() – QR-hajotelma vaihtoehdolla mode='reduced' (oletus) tai mode='complete'.

  • solve_triangular() – ratkaise A @ x = b, kun A on kolmiomatriisi.

  • cho_solve() – ratkaise A @ x = b, kun annetaan A:n Cholesky-tekijä.

6.13.2. dot, cross, trace

dot() on tapa, jolla matriisikertolasku kirjoitetaan kameralla. Operaattoria @, jonka työpöytä-numpy tarjoaa samaan tarkoitukseen, ei ole toteutettu ndarray-oliolle, joten jokainen matriisi-vektori-, matriisi-matriisi- ja vektoripistetulo kulkee dot()-kutsun kautta:

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

np.dot(a, b)             # 32.0 (scalar product)
np.cross(a, b)           # array([-3.0, 6.0, -3.0])

m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
np.dot(m, a)             # matrix-vector product
np.dot(m, m)             # matrix-matrix product
np.trace(m)              # 1 + 5 + 9 = 15.0

dot():n tulos on aina dtypeltään float. cross() on kahden 3-vektorin ristitulo, trace() neliömatriisin päälävistäjän summa.

6.13.3. inv ja det

m = np.array([[1, 2, 3, 4],
              [4, 5, 6, 4],
              [7, 9, 9, 4],
              [3, 4, 5, 6]])
print(np.linalg.inv(m))
print(np.linalg.det(m))

Käänteismatriisi lasketaan Gauss-Jordan-eliminoinnilla, joten inv() nostaa poikkeuksen ValueError, kun matriisi on singulaarinen (jokin lävistäjän alkio menee nollaksi eliminoinnin aikana). RAM-kustannus on suunnilleen kaksi kertaa syötteen koko.

Determinantti käyttää uudelleen samaa eliminointia – ajoaika on käytännössä sama kuin käänteismatriisilla.

Kun sovelluksen tavoitteena on ratkaista lineaarinen järjestelmä, älä käännä ja kerro – suosi alla olevia varsinaisia ratkaisijoita. Molemmat ovat nopeampia ja numeerisesti paremmin käyttäytyviä.

6.13.4. cholesky

Symmetriselle positiivisesti definiitille matriisille A cholesky() palauttaa alakolmiomatriisin L siten, että A = L @ L.T

a = np.array([[25, 15, -5],
              [15, 18,  0],
              [-5,  0, 11]])
L = np.linalg.cholesky(a)

Jos syöte ei ole positiivisesti definiitti tai ei ole symmetrinen, nostetaan poikkeus ValueError.

Cholesky-tekijä on puolet LU-tekijöinnin työstä ja on oikea lähtökohta mihin tahansa ongelmaan, jonka matriisi tiedetään symmetriseksi positiivisesti definiitiksi (kovarianssipäivitykset, pienimmän neliösumman sovituksen normaaliyhtälöt).

6.13.5. eig

eig() toimii vain reaalisilla symmetrisillä matriiseilla. Epäsymmetriset matriisit nostavat poikkeuksen ValueError. Se palauttaa 2-monikon (eigenvalues, eigenvectors)

a = np.array([[1, 2, 1, 4],
              [2, 5, 3, 5],
              [1, 3, 6, 1],
              [4, 5, 1, 7]], dtype=np.uint8)
x, y = np.linalg.eig(a)

Huomioita:

  • Ominaisarvot palautuvat ilman erityistä järjestystä. Käytä funktiota sort() (ja samaa permutaatiota ominaisvektoreihin funktion argsort() kautta), kun lajiteltu järjestys on tärkeä.

  • Ominaisvektori on yksikäsitteinen vain nollasta poikkeavaa skalaaria vaille, joten yksittäisten ominaisvektorien etumerkki ei ole yksikäsitteisesti määritelty. Kaksi oikeaa ajoa voivat tuottaa vastakkaismerkkisiä vektoreita; tämä on harmitonta.

6.13.6. norm

Vektorin tai matriisin euklidinen (Frobeniuksen) normi:

v = np.array([1, 2, 3, 4, 5])
m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

np.linalg.norm(v)        # 7.416...
np.linalg.norm(m)        # 16.881...

Valinnainen avainsana axis= ottaa normin yksittäisen akselin suunnassa koko taulukon sijaan.

6.13.7. qr

qr() jakaa suorakulmaisen matriisin A (muoto (M, N)) ortonormaaliksi matriisiksi Q ja yläkolmiomatriisiksi R siten, että A == Q @ R

A = np.arange(6).reshape((3, 2))

q, r = np.linalg.qr(A)
# mode='reduced' (default): q is (3, 2), r is (2, 2)

q, r = np.linalg.qr(A, mode='complete')
# q is (3, 3), r is (3, 2)

Hajotelma on toteutettu peräkkäisillä Givensin rotaatioilla. Oikea valinta pienimmän neliösumman ongelmiin, joissa matriisi ei ole symmetrinen.

6.13.8. Järjestelmien ratkaiseminen

Kaksi varsinaista ratkaisijaa ulab.scipy.linalg:n alla ovat molemmat nopeampia ja tarkempia kuin np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – ratkaise a @ x = b olettaen, että a on kolmiomatriisi:

    A = np.array([[3, 0, 0, 0],
              [2, 1, 0, 0],
              [1, 0, 1, 0],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
x = sp.linalg.solve_triangular(A, b, lower=True)

  • cho_solve(L, b)() – annettuna Cholesky-tekijä L, ratkaise A @ x = b, missä A = L @ L.T

    L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)

Käytä näitä kääntämisen sijaan aina, kun A:n rakenne sen sallii – ne säästävät eliminointityön ja eksplisiittisen käänteismatriisin.

6.13.9. Pienen lineaarisen järjestelmän ratkaiseminen

A = np.array([[3, 0, 1, 1],
              [0, 1, 0, 2],
              [1, 0, 1, 1],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])

x = np.dot(np.linalg.inv(A), b)
print(x)
print(np.dot(A, x))         # should equal b

Sama ongelma voidaan ilmaista jakamalla A tekijöihin ja kutsumalla sopivaa ratkaisijaa – nopeampaa ja tarkempaa, kun A:lla on oikea rakenne.

Täydellisen argumenttitason viitteen löydät kohdista numpy.linalg — Lineaarialgebran rutiinit ja scipy.linalg — Lineaarialgebran rutiinit.