6.13. Лінійна алгебра

Лінійна алгебра на камері – це робота з малими матрицями: поворот 3x3 для злиття зразка IMU у світовий кадр, калібрувальна матриця для випрямлення лінзи, оновлення коваріації стану фільтра Калмана, поліноміальна апроксимація, нормальні рівняння якої виражаються як невелика лінійна система. Підмодулі numpy.linalg та scipy.linalg покривають саме такий масштаб.

Функції розміщені у двох модулях. Операції матричного добутку знаходяться на верхньому рівні numpy; розкладання та обернена матриця – у numpy.linalg; спеціалізовані розв’язувачі лінійних систем – у scipy.linalg. Наприклад, множення матриць 2x2:

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. Що доступно

  • dot() – матричний або векторний добуток.

  • cross() – векторний добуток двох 3-D векторів.

  • trace() – сума діагональних елементів.

  • inv() – обернена матриця.

  • det() – визначник.

  • cholesky() – розкладання Холецького (симетрична позитивно визначена матриця на вході).

  • eig() – власні значення та власні вектори дійсної симетричної матриці.

  • norm() – 2-норма вектора або матриці.

  • qr() – QR-розкладання з mode='reduced' (за замовчуванням) або mode='complete'.

  • solve_triangular() – розв’язання A @ x = b, коли A є трикутною.

  • cho_solve() – розв’язання A @ x = b за відомим множником Холецького A.

6.13.2. dot, cross, trace

dot() – це спосіб запису матричного множення на камері. Оператор @, який надає настільний numpy для тієї ж задачі, не реалізований для ndarray, тому кожен добуток матриці на вектор, матриці на матрицю та скалярний добуток векторів виконується через виклик dot()

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() завжди має dtype float. cross() – це векторний добуток двох 3-векторів, trace() – сума головної діагоналі квадратної матриці.

6.13.3. inv та 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))

Обернена матриця обчислюється методом Гаусса-Жордана, тому inv() підіймає ValueError, коли матриця вироджена (діагональний елемент стає нулем під час виключення). Вартість оперативної пам’яті приблизно вдвічі більша від розміру вхідних даних.

Визначник повторно використовує те саме виключення – час виконання по суті такий самий, як і для оберненої матриці.

Коли мета застосунку – розв’язати лінійну систему, не слід обертати та множити – надайте перевагу спеціалізованим розв’язувачам нижче. Обидва швидші та чисельно стабільніші.

6.13.4. cholesky

Для симетричної позитивно визначеної матриці A, cholesky() повертає нижньотрикутну L таку, що A = L @ L.T

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

Якщо вхідні дані не є позитивно визначеними або не є симетричними, підіймається ValueError.

Множник Холецького вимагає вдвічі менше роботи, ніж LU-факторизація, і є правильною відправною точкою для будь-якої задачі, матриця якої відома як симетрична позитивно визначена (оновлення коваріації, нормальні рівняння з апроксимації методом найменших квадратів).

6.13.5. eig

eig() працює лише на дійсних симетричних матрицях. Несиметричні матриці підіймають ValueError. Функція повертає 2-кортеж (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)

Примітки:

  • Власні значення повертаються без певного порядку. Застосуйте sort() (і ту саму перестановку до власних векторів через argsort()), коли потрібен відсортований порядок.

  • Власний вектор визначений лише з точністю до ненульового скаляра, тому знак окремих власних векторів не є однозначно визначеним. Два коректних запуски можуть давати вектори з протилежними знаками; це нешкідливо.

6.13.6. norm

Евклідова (Фробеніусова) норма вектора або матриці:

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...

Необов’язкове ключове слово axis= обчислює норму вздовж одного осі замість усього масиву.

6.13.7. qr

qr() розкладає прямокутну матрицю A (форма (M, N)) на ортонормальну Q та верхньотрикутну R таку, що 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)

Розкладання реалізовано через послідовні повороти Гівенса. Правильний вибір для задач найменших квадратів, де матриця не є симетричною.

6.13.8. Розв’язання систем

Два спеціалізованих розв’язувачі з ulab.scipy.linalg є одночасно швидшими та точнішими, ніж np.dot(np.linalg.inv(A), b):

  • solve_triangular(a, b, lower=False)() – розв’язати a @ x = b за умови, що a є трикутною:

    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)() – за відомим множником Холецького L, розв’язати A @ x = b, де A = L @ L.T

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

Використовуйте їх замість обертання, коли структура A це дозволяє – вони економлять обчислення виключення і явне обернення.

6.13.9. Розв’язання невеликої лінійної системи

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

Ту саму задачу можна виразити, факторизувавши A та викликавши відповідний розв’язувач – це швидше та точніше, коли A має правильну структуру.

Для повної довідки на рівні аргументів, дивіться numpy.linalg — Процедури лінійної алгебри та scipy.linalg — Процедури лінійної алгебри.