Линейная алгебра#

Подмодуль scipy.linalg в целом является надмножеством над numpy.linalg, т.е. кроме того, что он содержит все методы из numpy.linalg, в него включены еще ряд дополнительных методов.

Ряд примеров приведен в документации.

import numpy as np
from scipy import linalg

Вычисление норм матриц и векторов#

Существуют разные нормы определенные на пространстве векторов и матриц. Многие из них можно вычислить средствами функции scipy.linalg.norm.

Векторы#

Для одномерных массивов по умолчанию возвращается евклидова норма (она же \(l_2\)), которая для \(x\in\mathbb{R}^n\) (или \(x\in\mathbb{C}^n\)) определяется формулой

\[ \|x\|_2 = \sqrt{\sum_{i=1}^n |x_i|^2}. \]
x = np.array([3, 4])

print(f"Евклидова норма вектора x равняется {linalg.norm(x)}")

Евклидова норма вектора x равняется 5.0

Евклидова норма является частным случаем нормы \(l_p\), которая для любого положительно \(p\) определяется выражением

\[ \|x\|_p = \left(\sum_{i = 1}^n |x_i|^p\right)^{\tfrac{1}{p}}. \]

Параметр ord позволяет указать значение p необходимой нормы. Например, очень часто вместо \(l_2\) (евклидова, \(p=2\)) используется называемая иногда абсолютной норма \(l_1\) (\(p=1\)), которая определяется выражением

\[ \|x\|_1 = \sum_{i=1}^n|x|. \]
print(f"Абсолютная норма вектора x равняется {linalg.norm(x, ord=1)}")

Абсолютная норма вектора x равняется 7.0

На другом конце спектра нормированных пространств семейства \(l_p\) находится пространство с равномерной метрикой \(l_\infty\), определяемой выражением

\[ \|x\|_\infty = \max_{i=1}^n |x_i|. \]

Чтобы её вычислить, необходимо передать в качестве параметра ord значение numpy.inf.

print(f"Равномерная норма вектора x равняется {linalg.norm(x, ord=np.inf)}")

Равномерная норма вектора x равняется 4.0

Матрицы#

Если массив двухмерный, то по умолчанию вычисляется тоже евклидова норма (она же норма Фробениуса), которая для матрицы \(A\in\mathbb{R}^{n\times m}\) определяется выражением

\[ \|A\| = \sqrt{\sum_{i=1}^n\sum_{j=1}^m A_{ij}^2}. \]

Tip

Если необходимо вычислить векторную норму от множества векторов в многомерном массиве, то полезным окажется параметр axis этого метода.

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

print(f"Норма матрицы A равняется {linalg.norm(A)}")

Норма матрицы A равняется 5.477225575051661

В ряде задач, в которых матрицы используются для представления линейных операторов, важно следить за согласованностью и подчиненностью матричных норм с соответствующими нормами векторных пространств, в которых они оперируют.

Так, например, норма евклидово векторное пространство \(l_2\) порождает в матричном пространстве спектральную норму, которая равняется максимальному сингулярному числу этого оператора. Чтобы вычислить спектральную норму матрицы, необходимо передать 2 в качестве параметра ord.

singular_values = linalg.svdvals(A)

print(f"Максимальное спектральное число матрицы A равняется {singular_values.max()}")
print(f"Спектральная норма матрицы m равняется {linalg.norm(A, ord=2)}")

Максимальное спектральное число матрицы A равняется 5.464985704219043
Спектральная норма матрицы m равняется 5.464985704219043

За более полным списком поддерживаемых норм этим методом, автор направляет вас в документацию.

Обращение матриц#

Функция linalg.inv отвечает за обращение матрицы, т.е. за поиск такой матрицы \(A^{-1}\), что \(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]]

Чтобы убедиться, что она возвращает обратную матрицу, проверим выполнение равенства \(A A^{-1} = A^{-1} A = I\).

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

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

Сверять целый массив глазами не очень удобно. Упростим себе задачу и посчитаем норму невязки, т.к. нормы матриц

\[ R_1 = \| A A^{-1} - I \| \]

и

\[ R_2 = \| A^{-1} A - I \|. \]
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_2=}")

R_1=9.930136612989092e-16, R_2=4.965068306494546e-16

Решение систем линейных алгебраических уравнений#

Функция 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.

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"{x_approx=}. Погрешность: {error}")

x_approx=array([1., 1.]). Погрешность: 4.002966042486721e-16

Кроме метода solve, который пытается решить СЛАУ в общем случае, в NumPy предусмотрен набор методов для решения СЛАУ, если известны свойства матрицы в левой части. Например, для треугольной матрицы предпочтительно использовать метод solve_triangular, для циркулянтной матрицы — метод solve_circulant и т.п.

Определитель квадратной матрицы вычисляется методом det.

print(f"Определитель матрицы A равен {linalg.det(A)}")

Определитель матрицы A равен -2.0

Поиск собственных чисел и векторов#

В ряде задач требуется найти собственные числа и/или собственные значения.

Метод eigvals отвечает за поиск собственных чисел. Найдем с помощью него собственные значения матрицы

\[\begin{split} A = \begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}. \end{split}\]

Характеристический полином

\[\begin{split} |A-\lambda I| = \begin{pmatrix} 2 - \lambda & 1\\ 1 & 2 - \lambda \end{pmatrix} = \lambda^2 - 4 \lambda + 3 \end{split}\]

имеет два действительных решения \(\lambda_1=3\) и \(\lambda_2=1\).

A = np.array([
    [2, 1],
    [1., 2]
])

eigen_values = linalg.eigvals(A)
print(f"Собственные значения матрицы A: {eigen_values}")

Собственные значения матрицы A: [3.+0.j 1.+0.j]

Метод eig вычисляет сразу и собственные числа и собственные векторы. При этом собственные вектора возвращаются как столбцы матрицы второго возвращаемого значения.

eigen_values, eigen_vectors = linalg.eig(A)
print(f"Собственные значения матрицы A: {eigen_values}")

for vector in eigen_vectors.T:
    print(f"{vector} --- собственный вектор матрицы A")

Собственные значения матрицы A: [3.+0.j 1.+0.j]
[0.70710678 0.70710678] --- собственный вектор матрицы A
[-0.70710678  0.70710678] --- собственный вектор матрицы A

Проверим, что найденные собственные вектора и числа корректны. Для этого посчитаем невязки

\[ R_i = \| Av_i - \lambda_i v_i \|, \, i=1,2, \]

где \(v_1\) и \(v_2\) — собственные вектора соответствующие собственным числам \(\lambda_1\) и \(\lambda_2\).

for l, v in zip(eigen_values, eigen_vectors.T):
    R = linalg.norm(A@v - l*v)
    print(f"{R=}")

R=0.0
R=0.0

Tip

Если заранее известны определенные свойства матрицы \(A\), то лучше воспользоваться соответствующими методами.

Сингулярное разложение матрицы#

За поиск сингулярных чисел отвечает функция svdvalues. Функция svd (Singular Value Decomposition) вычисляет сингулярное разложение матрицы, которое определено и для прямоугольных матриц.

A = np.array([
    [1, 2, 3],
    [4, 5, 6],
])
left_vectors, singular_values, right_vectors = linalg.svd(A)

for v, s, u in zip(left_vectors.T, singular_values, right_vectors):
    R = linalg.norm(A.T @ v - s * u)
    print(f"{R=}")

R=1.9860273225978185e-15
R=1.1625080825424498e-15

Tip

Кроме сингулярного разложения в scipy есть и многие другие, такие как QR-разложение, LU-разложение и т.д.