6.13. Линейная алгебра

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

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

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() – векторное произведение трёхмерных векторов.

  • 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() – это векторное произведение двух трёхмерных векторов, 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, когда матрица вырождена (диагональный элемент становится нулём в ходе исключения). Затраты RAM примерно вдвое превышают размер входной матрицы.

Определитель использует то же самое исключение – время выполнения по сути такое же, как у обращения.

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

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. Функция возвращает кортеж из двух элементов (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 — Процедуры линейной алгебры.