Разложения матриц: LU, QR, SVD и собственные значения

За магией «решить систему» и «понизить размерность» стоят разложения матриц. Разберём, что считает каждое и зачем.

Разложение (факторизация) матрицы — представление матрицы A как произведения нескольких матриц с особыми свойствами, с которыми задачу решать проще.

Идея: разложить, чтобы упростить

Как число 12 раскладывают на множители 2·2·3, матрицу раскладывают на «удобные» сомножители. Удобные — значит треугольные, ортогональные или диагональные: с ними системы решаются мгновенно, а структура данных становится видимой. Каждое разложение — инструмент под свой класс задач.

Четыре главных разложения

LU: «решить один раз — переиспользовать»

LU разбивает A на нижнетреугольную L и верхнетреугольную U. Зачем? Если надо решить A·x = b для десяти разных b (например, та же конструкция под разными нагрузками), разложение делается один раз — O(n³), а каждое последующее решение стоит дёшево — O(n²). Именно LU использует scipy.linalg.solve из прошлого урока.

from scipy.linalg import lu, solve
P, L, U = lu(A)        # A = P @ L @ U (P — перестановка строк)
# теперь любую систему A x = b решаем быстро через L и U

SVD: королева разложений

Сингулярное разложение A = U·Σ·Vᵀ работает для любой матрицы (даже не квадратной) и раскрывает её «скелет»: сингулярные числа на диагонали Σ показывают, по каким направлениям данные растянуты сильнее всего. Отбросив маленькие сингулярные числа, получаем сжатие и понижение размерности — на этом стоит метод главных компонент (PCA) и сжатие изображений.

Собственные значения «руками»: метод степеней

Наибольшее по модулю собственное значение можно найти простой итерацией: умножаем вектор на матрицу много раз — он «поворачивается» к главной собственной оси. Вот метод степеней (power iteration) на чистом Python:

A = [[2.0, 1.0],
     [1.0, 2.0]]          # собств. значения: 3 и 1

v = [1.0, 0.0]            # стартовый вектор
lam = 0.0
for _ in range(50):
    # w = A v
    w = [A[0][0]*v[0] + A[0][1]*v[1],
         A[1][0]*v[0] + A[1][1]*v[1]]
    # длина и нормировка
    norm = (w[0]**2 + w[1]**2) ** 0.5
    v = [w[0]/norm, w[1]/norm]
    lam = norm

print("Наибольшее собств. значение ~", round(lam, 6))

Вывод:

Наибольшее собств. значение ~ 3.0

А в SciPy — одна строка

import numpy as np
from scipy.linalg import eig, svd

vals, vecs = eig(np.array([[2, 1], [1, 2]]))
print(vals)            # [3.+0.j  1.+0.j]

U, s, Vt = svd(np.array([[2, 1], [1, 2]]))
print(s)               # [3. 1.] — сингулярные числа

Как работает под капотом

Метод степеней сходится к наибольшему собственному значению, потому что при каждом умножении компонента вдоль главной собственной оси растёт быстрее остальных (в λ_max/λ_2 раз за шаг). Нормировка не даёт вектору «улететь» в бесконечность, а длина w и есть оценка λ. Настоящие алгоритмы (QR-алгоритм в LAPACK) находят сразу все собственные значения и устойчивее, но идея «итерируй и наблюдай, куда поворачивается вектор» — та же.

Частые ошибки

• Считать собственные значения определителем «вручную» для больших n. Это численно неустойчиво; используйте eig.
• Ждать вещественных собственных значений всегда. У несимметричных матриц они бывают комплексными (+0.j).
• Путать собственные и сингулярные числа. Они совпадают только для симметричных положительно определённых матриц.

Итог

• Разложения раскладывают матрицу на удобные сомножители под конкретную задачу.
• LU — для повторного решения СЛАУ; QR — для МНК; SVD — для понижения размерности.
• Собственные значения описывают «оси» и устойчивость системы; метод степеней даёт главное из них.
• В SciPy всё это — функции lu, qr, svd, eig из scipy.linalg.