Линейная алгебра
Contents
Линейная алгебра¶
Подмодуль scipy.linalg в целом является надмножеством над numpy.linalg, т.е. кроме того, что он содержит все методы из numpy.linalg, в него включены еще ряд дополнительных методов.
Ряд примеров приведен в документации.
Обращение матриц (функция linalg.inv) и решение СЛАУ (функция linalg.solve);
Разложение матрицы: собственные значения и вектора, сингулярное разложение, LU и QR разложения и другое;
Методы создания специальных матриц, блочно-диагональных и многих других;
Обращение матриц¶
Функция linalg.inv позволяет обратить матрицу.
Чтобы убедиться, что она возвращает обратную матрицу, проверим выполнение равенства \(A A^{-1} = A^{-1} A = I\), где \(I\) — единичная матрица.
import numpy as np
from scipy import linalg
A = np.array([[1., 2.], [3., 4.]])
A_inv = linalg.inv(A)
print(A, A_inv, sep="\n")
[[1. 2.]
[3. 4.]]
[[-2. 1. ]
[ 1.5 -0.5]]
print(A @ A_inv)
print(A_inv @ A)
[[1.0000000e+00 0.0000000e+00]
[8.8817842e-16 1.0000000e+00]]
[[1.00000000e+00 0.00000000e+00]
[2.22044605e-16 1.00000000e+00]]
Как видим, получилась не совсем единичная матрица, но т.к. все вычисления с числами с плавающей точкой производятся в компьютере не точно, большего ожидать нельзя.
Сверять целый массив глазами не очень удобно. Давайте упростим себе задачу и посчитаем норму невязки:
и
Для этого воспользуемся функцией scipy.linalg.norm.
R_1 = linalg.norm(A @ A_inv - np.eye(2))
R_2 = linalg.norm(A_inv @ A - np.eye(2))
print(f"R_1 = {R_1}, R_2 = {R_2}" )
R_1 = 9.930136612989092e-16, R_2 = 4.965068306494546e-16
Note
Функция linalg.norm позволяет считать разные нормы, в том числе и матричные.
Если эта функция встречает вектор, то по умолчанию вычисляется евклидова норма, если матрицу, то по умолчанию вычисляется норма Фробениуса.
Решение СЛАУ¶
Функция linalg.solve используется для решения СЛАУ вида \(Ax = b\), где \(A\) — квадратная матрица размера \(N \times N\). Правая часть уравнения \(b\) может быть формы (N, NRHS), т.е. если \(b\) — одномерный массив, то решается обычное СЛАУ, а если \(b\) — матрица из \(N\) строк и NRHS столбцов, то решается матричное уравнение (или сразу решается NRHS СЛАУ с общей матрицей \(A\), а в качестве векторов правых частей выступают столбцы матрицы \(b\)).
Для иллюстрации создадим матрицу A и вектор x. Сгенерируем вектор правой части b = A @ x и попробуем восстановить значение x, решив СЛАУ A x = b.
import numpy as np
from scipy import linalg
A = np.array([[1., 2.], [3., 4.]])
x = np.array([1., 1.])
b = A @ x
x_approx = linalg.solve(A, b)
error = linalg.norm(x_approx - x)
print(f"Погрешность: {error}")
Погрешность: 4.002966042486721e-16
Проиллюстрируем одновременное решение нескольких СЛАУ. Одновременно решим уравнение \(A x = b\) и \(A y = c\), где \(A\) — общая матрица, а правые части \(b\) и \(c\) — разные. Для этого составим матрицу \(B\), которая состоит из столбцов \(b\) и \(c\). Решение \(X\) матричного уравнения \(A X = B\) — матрица, столбцы которой — решения исходных СЛАУ \(x\) и \(y\).
y = np.array([42., 42.])
c = A @ y
print(f"Точные решения: {x} и {y}")
B = np.vstack([b, c]).T
X = linalg.solve(A, B)
x, y = X.T
print(f"Численные: {x} и {y}")
Точные решения: [1. 1.] и [42. 42.]
Численные: [1. 1.] и [42. 42.]