Собственные значения, SVD и матричные разложения

В предыдущей статье трека — Матрицы, линейные отображения, определители и системы уравнений — матрица была таблицей чисел, которая умеет умножаться на вектор. Этого достаточно, чтобы считать, но недостаточно, чтобы понимать: матрица 1000×1000 — это миллион чисел, и, глядя на них, нельзя сказать ровным счётом ничего о том, что она делает.

Разложения — способ переписать матрицу произведением простых сомножителей, каждый с ясным геометрическим или вычислительным смыслом. A = U Σ Vᵀ — те же данные, но теперь видно: вот повороты, вот растяжения, вот сколько в матрице «настоящей» информации, а сколько шума. Как разложение числа на простые: 2⁴·3·17 говорит о числе 816 больше, чем его десятичная запись.

Это, пожалуй, самая окупаемая математика во всём треке: PCA, рекомендательные системы, PageRank, устойчивость численных схем, сжатие эмбеддингов, LoRA в дообучении LLM — буквально одна и та же теорема на разных данных.

Зачем это программисту: четыре сцены

Сцена 1. Регрессия сходит с ума. Вы решаете AᵀA x = Aᵀb, коэффициенты скачут на порядки от добавления одной строки данных. Диагноз: cond(A) огромно, а нормальные уравнения возводят его в квадрат. Лечение — SVD или QR, в одну строчку.

Сцена 2. Эмбеддинги весят 40 ГБ. Матрица 10⁷×768 не влезает в память ноды. Усечённое SVD до 128 компонент даёт 6× сжатие и теряет 3% качества поиска — потому что спектр реальных данных убывает быстро.

Сцена 3. Симуляция взрывается. Схема x_{k+1} = M x_k расходится. Спектральный радиус ρ(M) = 1.02 — всё, теорема сказала всё за вас: схема неустойчива при любом шаге и любом типе float.

Сцена 4. Надо разбить граф на кластеры. Считаете лапласиан L = D − W, берёте второй снизу собственный вектор, режете по знаку — получаете разрез, близкий к минимальному.

Во всех четырёх случаях вопрос один: какие направления в пространстве матрица трогает особым образом.

Часть 1. Собственные значения и векторы

Интуиция: оси, которые преобразование не сбивает

Возьмём растяжение плоскости в 3 раза по горизонтали. Почти любой вектор при этом и повернётся, и изменит длину: (1, 1) уйдёт в (3, 1) — это другое направление. Но есть исключительные направления: (1, 0) уходит в (3, 0) (то же направление, длина ×3), (0, 1) остаётся на месте (×1).

«Направления, которые преобразование не сбивает, а только растягивает» — собственные векторы; коэффициенты растяжения — собственные значения. Найти их значит найти естественную систему координат самой матрицы, в которой она выглядит максимально просто: не таблица из миллиона чисел, а «умножить каждую ось на свой коэффициент».

Строгое определение

Пусть A — квадратная матрица n×n над полем F (обычно R или C).

Число λ ∈ F называется собственным значением A,
если существует вектор v ≠ 0 такой, что

    A v = λ v.

Такой v — собственный вектор, отвечающий λ.
Множество всех собственных значений — спектр матрицы, spec(A).

Требование v ≠ 0 критично: иначе равенство выполнялось бы для любого λ и определение стало бы пустым.

Перепишем: (A − λI) v = 0. Ненулевое решение у однородной системы есть тогда и только тогда, когда A − λI вырождена, то есть det(A − λI) = 0. Слева — многочлен степени n от λ, характеристический многочлен χ_A(λ); его корни и есть собственные значения. Над C их ровно n с учётом кратности (основная теорема алгебры); над R их может не быть вовсе — поворот на 90° не оставляет ни одного направления на месте, его собственные значения равны ±i.

Собственное подпространство E_λ = ker(A − λI); его размерность — геометрическая кратность. Кратность λ как корня χ_A — алгебраическая кратность. Всегда 1 ≤ геом ≤ алг.

Ручной разбор 2×2

A = [[2, 1], [1, 2]]

det(A − λI) = (2−λ)² − 1 = λ² − 4λ + 3 = (λ − 1)(λ − 3)   ⇒  λ₁ = 3, λ₂ = 1

λ₁ = 3:  (A − 3I) = [[−1, 1], [1, −1]],  ядро: v₁ = v₂  ⇒  v = (1,  1)
λ₂ = 1:  (A − 1I) = [[ 1, 1], [1,  1]],  ядро: v₁ = −v₂ ⇒  v = (1, −1)

Геометрия: матрица растягивает диагональ (1,1) втрое и оставляет антидиагональ (1,−1) на месте. Собственные векторы ортогональны — не случайность, а следствие симметричности A (спектральная теорема ниже).

import numpy as np

A = np.array([[2.0, 1.0], [1.0, 2.0]])

# eigh — для симметричных/эрмитовых: вещественные λ и ортонормированный базис
w, V = np.linalg.eigh(A)
print(w)                      # [1. 3.] — по возрастанию
for i in range(2):            # проверка определения A v = λ v
    assert np.allclose(A @ V[:, i], w[i] * V[:, i])

Два бесплатных инварианта: след и определитель

По теореме Виета для χ_A:

tr(A) = λ₁ + … + λₙ      (сумма диагонали = сумма собственных значений)
det(A) = λ₁ · … · λₙ     (определитель = произведение)

Для примера выше: tr = 4 = 3 + 1, det = 3 = 3·1 — дешёвая проверка любого расчёта. Следствие: матрица вырождена ⟺ 0 является её собственным значением. И ещё: спектр не зависит от базиса — при подобии A → P⁻¹AP многочлен не меняется, так как det(P⁻¹(A − λI)P) = det(A − λI).

Спектральный радиус ρ(A) = max|λᵢ|. Теорема: Aᵏ → 0 ⟺ ρ(A) < 1. Это критерий устойчивости любой итерации x_{k+1} = A x_k + b — метода Якоби, Гаусса–Зейделя, разностной схемы, замкнутого контура управления (подробнее — Численные методы). Тонкость: ρ(A) < 1 гарантирует сходимость лишь асимптотически, а на первых шагах норма несимметричной матрицы может расти (transient growth). Поэтому смотрят и на ‖A‖₂ = σ₁.

Диагонализация: когда матрица «распадается» на числа

Если есть базис из n линейно независимых собственных векторов, соберём их в столбцы P и положим D = diag(λ₁, …, λₙ):

A P = P D    ⇒    A = P D P⁻¹        (спектральное разложение)
Aᵏ = P Dᵏ P⁻¹,   где Dᵏ = diag(λ₁ᵏ, …, λₙᵏ)

Это смена координат: P⁻¹ переводит вектор в собственный базис, D масштабирует каждую ось, P возвращает обратно. Возведение матрицы в степень превращается в возведение n чисел. Так же определяются функции от матрицы: exp(A) = P exp(D) P⁻¹. Матричная экспонента — решение линейного ОДУ x' = A x, то есть непрерывная динамика системы.

Критерий: A диагонализуема ⟺ для каждого λ геометрическая кратность равна алгебраической. Достаточное (не необходимое) условие: все n значений различны.

Дефектные матрицы: когда всё ломается

J = [[2, 1], [0, 2]]

χ_J(λ) = (2 − λ)²  ⇒  λ = 2 алгебраической кратности 2
(J − 2I) = [[0, 1], [0, 0]]  — ядро одномерно: v = (1, 0)

Геометрическая кратность 1, алгебраическая 2: базиса из собственных векторов нет, P вырождена. Такие матрицы называются дефектными; максимум, что достижимо, — жорданова форма с единицами над диагональю.

J = np.array([[2.0, 1.0], [0.0, 2.0]])
w, V = np.linalg.eig(J)
print(w)                            # [2. 2.]
print(np.linalg.matrix_rank(V))     # 1 — столбцы коллинеарны, базиса нет

Практический вывод: жорданова форма — прекрасный теоретический инструмент и негодный вычислительный. Она разрывна по элементам: сколь угодно малое возмущение делает матрицу диагонализуемой с чудовищно плохо обусловленной P. В численных библиотеках её не считают вообще — вместо неё разложение Шура (см. ниже).

Пример: Фибоначчи через спектр

Рекуррента F_{n+1} = F_n + F_{n−1} — линейное отображение с матрицей [[1,1],[1,0]]. Её собственные значения — корни λ² − λ − 1 = 0: золотое сечение φ ≈ 1.618 и ψ ≈ −0.618. Диагонализация мгновенно даёт формулу Бине F_n = (φⁿ − ψⁿ)/√5.

F = np.array([[1, 1], [1, 0]], dtype=float)
print(np.linalg.eigvals(F))         # [1.618, -0.618]

phi, psi = (1 + 5**0.5) / 2, (1 - 5**0.5) / 2
binet  = [round((phi**n - psi**n) / 5**0.5) for n in range(1, 11)]
matpow = [int(round(np.linalg.matrix_power(F, n)[0, 1])) for n in range(1, 11)]
print(binet == matpow, binet)       # True [1, 1, 2, 3, 5, 8, 13, 21, 34, 55]

|ψ| < 1, поэтому вклад второго слагаемого экспоненциально гаснет: F_n ≈ φⁿ/√5. Асимптотика любой линейной рекурренты определяется наибольшим по модулю собственным значением — тот же механизм работает при решении рекуррент в анализе алгоритмов (Дискретная математика и комбинаторика).

Пример: марковские цепи, PageRank и λ = 1

Пусть P — стохастическая матрица переходов (строки суммируются в 1), состояние — вектор-строка вероятностей: x_{k+1} = x_k P. Стационарное распределение π есть решение π P = π, то есть левый собственный вектор при λ = 1.

P = np.array([
    [0.90, 0.07, 0.03],   # Free
    [0.15, 0.80, 0.05],   # Trial
    [0.25, 0.25, 0.50],   # Churn
])
assert np.allclose(P.sum(axis=1), 1.0)

w, V = np.linalg.eig(P.T)              # левые собственные векторы = правые у Pᵀ
i = np.argmin(np.abs(w - 1.0))
pi = np.real(V[:, i]); pi /= pi.sum()
print(pi)                              # [0.6272 0.3047 0.0681]

x = np.array([1.0, 0.0, 0.0])          # тот же ответ степенным методом
for _ in range(200):
    x = x @ P
print(x)                               # [0.6272 0.3047 0.0681]
print(sorted(np.abs(w))[-2])           # 0.743 — второе по модулю λ

Теорема Перрона–Фробениуса гарантирует: для неразложимой стохастической матрицы с положительными элементами λ₁ = 1 — единственное максимальное по модулю значение, а его собственный вектор положителен. Скорость сходимости степенного метода задаёт спектральная щель |λ₂/λ₁|: здесь 0.743, ошибка падает вчетверо за 5 шагов.

Это буквально PageRank. Google-матрица G = d·P + (1−d)/n · J с d = 0.85 нужна затем, чтобы гарантировать неразложимость (нет висячих страниц-поглотителей) и одновременно зажать |λ₂| ≤ d — то есть обеспечить предсказуемые 50–100 итераций на графе из миллиардов вершин. Оригинал: Page, Brin et al., «The PageRank Citation Ranking» (http://ilpubs.stanford.edu:8090/422/).

Спектральная теорема: почему симметричные матрицы — рай

Если A = Aᵀ (вещественная симметричная), то:

1. Все собственные значения вещественны.
2. Собственные векторы для разных λ ортогональны.
3. Существует ОРТОНОРМИРОВАННЫЙ базис из собственных векторов:
       A = Q Λ Qᵀ,   QᵀQ = I,   Λ = diag(λᵢ)
4. Дефектных симметричных матриц не бывает.

Доказательство пункта 2 в одну строку: пусть Av = λv, Aw = μw, λ ≠ μ. Тогда λ(v·w) = (Av)·w = v·(Aᵀw) = v·(Aw) = μ(v·w), значит (λ−μ)(v·w) = 0, значит v·w = 0.

Ценность пункта 3 огромна: Q⁻¹ = Qᵀ, обращать нечего, а cond₂(Q) = 1 — идеально. Все ужасы предыдущего раздела относятся к несимметричному случаю. Поэтому в статистике и ML почти всё сводят к симметричным матрицам: ковариация XᵀX, гессиан, матрица Грама, лапласиан графа.

Положительная определённость и отношение Рэлея

A симметрична и положительно определена (SPD): xᵀAx > 0 для всех x ≠ 0  ⟺  все λᵢ > 0
Положительно полуопределена (PSD):             xᵀAx ≥ 0                 ⟺  все λᵢ ≥ 0

Отношение Рэлея R(x) = (xᵀAx)/(xᵀx) для симметричной A пробегает ровно отрезок [λ_min, λ_max], достигая концов на соответствующих собственных векторах. Это вариационная характеризация: старшее собственное значение — максимум квадратичной формы на единичной сфере. Отсюда растут и степенной метод, и PCA («направление максимальной дисперсии»), и оценки обусловленности.

Где SPD встречается в проде:

• Ковариационные матрицы — PSD по построению. Если численно вылезло отрицательное λ (регулярно бывает на почти вырожденных данных), это ошибка округления: лечится клиппингом спектра или регуляризацией C + εI.
• Гессиан в минимуме — PSD. Отрицательные λ означают седловую точку, а не минимум: главный сюжет оптимизации нейросетей.
• Ядровые матрицы в SVM и гауссовских процессах обязаны быть PSD, иначе задача теряет выпуклость.
• Матрицы жёсткости в МКЭ — SPD, что и позволяет применять Холецкого и сопряжённые градиенты.

Часть 2. Сингулярное разложение (SVD)

Проблема, которую решает SVD

Собственные значения есть только у квадратных матриц; у несимметричных они комплексные, базис может отсутствовать вовсе, а P — быть плохо обусловленной. Реальные же матрицы данных прямоугольные: 10⁶ пользователей × 10⁴ товаров, 50 000 документов × 200 000 слов. SVD — разложение, которое существует для любой матрицы вообще. Плата за универсальность: вместо одного базиса используются два.

Геометрия: круг всегда переходит в эллипс

Ключевая геометрическая теорема: образ единичной сферы при любом линейном отображении — эллипсоид (возможно, сплющенный). SVD просто выписывает его оси.

Читается справа налево по формуле A = U Σ Vᵀ: Vᵀ поворачивает пространство так, чтобы особые направления v₁, v₂ легли на координатные оси; Σ растягивает вдоль осей на σ₁, σ₂; U доворачивает результат. Любое линейное отображение — это «поворот → масштаб по осям → поворот», других возможностей нет.

Определение и существование

Для любой A размера m×n существуют
    U — ортогональная m×m,
    Σ — m×n «диагональная», σ₁ ≥ σ₂ ≥ … ≥ σ_min(m,n) ≥ 0,
    V — ортогональная n×n,
такие что  A = U Σ Vᵀ.

σᵢ — сингулярные числа; столбцы U — левые, столбцы V — правые сингулярные векторы.

Существование выводится из спектральной теоремы: AᵀA симметрична и PSD, значит AᵀA = V Λ Vᵀ с λᵢ ≥ 0. Положим σᵢ = √λᵢ и uᵢ = A vᵢ / σᵢ; ортонормированность uᵢ проверяется в две строки, а A vᵢ = σᵢ uᵢ верно по построению. Отсюда фундаментальные связи:

AᵀA = V ΣᵀΣ Vᵀ   ⇒  правые сингулярные векторы = собственные векторы AᵀA
AAᵀ = U ΣΣᵀ Uᵀ   ⇒  левые  сингулярные векторы = собственные векторы AAᵀ
σᵢ = √λᵢ(AᵀA)

Важная деталь: AᵀA считать не надо — возведение в квадрат удваивает потерю точности (см. ниже). Настоящие алгоритмы работают с A напрямую (бидиагонализация Голуба–Кахана).

Ручной разбор SVD 2×2

A = [[3, 0], [4, 5]],   Aᵀ = [[3, 4], [0, 5]]

AᵀA = [[3·3+4·4, 3·0+4·5], [0·3+5·4, 0·0+5·5]] = [[25, 20], [20, 25]]

χ(λ) = (25−λ)² − 400  ⇒  25−λ = ±20  ⇒  λ = 45, 5
σ₁ = √45 = 3√5 ≈ 6.708,   σ₂ = √5 ≈ 2.236

Собственные векторы AᵀA:  v₁ = (1, 1)/√2 (λ=45),  v₂ = (1, −1)/√2 (λ=5)

u₁ = A v₁ / σ₁ = (3, 9)/√2 / (3√5) = (1,  3)/√10 ≈ (0.316,  0.949)
u₂ = A v₂ / σ₂ = (3, −1)/√2 / √5   = (3, −1)/√10 ≈ (0.949, −0.316)

Проверка: σ₁σ₂ = 3√5 · √5 = 15 = |det A|.  ✓

A = np.array([[3.0, 0.0], [4.0, 5.0]])
U, s, Vt = np.linalg.svd(A)
print(s)                                    # [6.7082 2.2361] = [3√5, √5]
print(np.allclose(U @ np.diag(s) @ Vt, A))  # True
print(s.prod(), abs(np.linalg.det(A)))      # 15.0 15.0

Знаки столбцов U и строк Vt numpy может выдать противоположными — SVD определено с точностью до знака (и до поворота внутри собственного подпространства при совпадающих σ). Поэтому не сравнивайте сингулярные векторы поэлементно в тестах: сравнивайте A ≈ UΣVᵀ или подпространства.

Что SVD говорит о матрице бесплатно

rank(A)       = число ненулевых σᵢ
‖A‖₂          = σ₁                    (спектральная норма)
‖A‖_F         = √(σ₁² + … + σᵣ²)      (норма Фробениуса)
cond₂(A)      = σ₁ / σ_min            (число обусловленности)
|det A|       = σ₁ · … · σₙ           (для квадратной)
col-space(A)  = span(u₁, …, uᵣ)
null-space(A) = span(v_{r+1}, …, vₙ)

Последние две строки — «четыре фундаментальных подпространства» Стрэнга из одного разложения. Критический практический момент: численный ранг никогда не бывает точным. Матрица математического ранга 3 в float64 даст σ₄ ≈ 1e−16·σ₁ вместо нуля, поэтому np.linalg.matrix_rank считает не нули, а количество σᵢ выше порога σ₁ · max(m,n) · eps. Проверка вырожденности через det(A) == 0 — грубая ошибка; через сингулярные числа — правильный способ.

Малоранговое приближение: теорема Эккарта–Янга–Мирского

Вероятно, самая практически ценная теорема во всей прикладной линейной алгебре.

Пусть A = Σ σᵢ uᵢ vᵢᵀ (сумма r слагаемых ранга 1).
Обрежем на k членах:  Aₖ = Σ_{i≤k} σᵢ uᵢ vᵢᵀ.

Тогда для ЛЮБОЙ матрицы B ранга ≤ k:
    ‖A − B‖₂  ≥ ‖A − Aₖ‖₂  = σ_{k+1}
    ‖A − B‖_F ≥ ‖A − Aₖ‖_F = √(σ_{k+1}² + … + σᵣ²)

Усечённое SVD — доказуемо оптимальное приближение матрицей заданного ранга, и ошибка известна заранее по спектру. Ни один другой метод сжатия ранга не может быть лучше.

rng = np.random.default_rng(0)
M = rng.standard_normal((60, 40))
U, s, Vt = np.linalg.svd(M, full_matrices=False)

k = 10
Mk = (U[:, :k] * s[:k]) @ Vt[:k]              # эффективнее, чем U @ diag(s) @ Vt

print(np.linalg.norm(M - Mk, 2), s[k])                            # 9.4083  9.4083
print(np.linalg.norm(M - Mk, 'fro'), np.sqrt((s[k:]**2).sum()))   # 32.148  32.148

Совпадение до последнего знака — теорема выполняется точно. На случайной гауссовой матрице спектр убывает медленно и сжатие бессмысленно; на реальных данных (пользователь×товар, документ×слово, пиксели) спектр падает почти экспоненциально, и k = 50 из 10 000 объясняет 90% энергии. Это содержательное утверждение о мире: реальные матрицы почти малоранговые, потому что за ними стоит небольшое число скрытых факторов.

PCA — это SVD центрированной матрицы

1. Центрируем: Xc = X − mean(X, axis=0)
2. SVD: Xc = U Σ Vᵀ
3. Главные компоненты (направления) — строки Vᵀ
4. Дисперсия вдоль i-й компоненты = σᵢ²/(N−1)
5. Проекция на k компонент = Uₖ Σₖ

X = rng.standard_normal((200, 5)) @ rng.standard_normal((5, 5))   # коррелированные признаки
Xc = X - X.mean(axis=0)

U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
evr_svd = s**2 / (s**2).sum()                    # explained variance ratio

C = np.cov(Xc, rowvar=False)
w, _ = np.linalg.eigh(C)
evr_cov = w[::-1] / w.sum()

print(evr_svd)   # [0.6053 0.2224 0.114  0.0571 0.0011]
print(evr_cov)   # [0.6053 0.2224 0.114  0.0571 0.0011] — то же самое
print(np.allclose(s**2 / (len(X) - 1), w[::-1]))   # True

Математически пути эквивалентны, численно — нет: через ковариацию вы явно строите XᵀX и теряете вдвое больше значащих цифр. Правило: мало признаков и хорошая обусловленность — можно через eigh(cov); во всех остальных случаях SVD. Так и устроен sklearn.decomposition.PCA (https://scikit-learn.org/stable/modules/decomposition.html).

Три классические ошибки в PCA:

1. Забыли центрировать — первая компонента просто указывает на среднее, разложение бессмысленно.
2. Не отмасштабировали признаки — «зарплата в рублях» с дисперсией 10¹⁰ раздавит «возраст» с дисперсией 100. PCA не инвариантно к масштабу.
3. Думают, что PCA выбирает важные признаки. Нет: он выбирает направления максимальной дисперсии, а дисперсия не равна полезности. Направление, различающее классы, может иметь крошечную дисперсию и быть отброшенным. Для «важности относительно цели» нужны LDA, PLS или feature selection.

Псевдообратная и МНК: главная численная ловушка

Задача наименьших квадратов: минимизировать ‖Ax − b‖₂. Учебный ответ — нормальные уравнения AᵀA x = Aᵀb. Ловушка: cond(AᵀA) = cond(A)². При cond(A) = 10⁸ (рядовое дело для полиномиальных признаков) получаем 10¹⁶ — ровно машинная точность float64, от ответа не остаётся ни одной верной цифры.

n, t = 12, np.linspace(0, 1, 50)
A = np.vander(t, n, increasing=True)          # матрица Вандермонда — плохо обусловлена
x_true = np.ones(n)
b = A @ x_true

print(np.linalg.cond(A))        # 1.17e+08
print(np.linalg.cond(A.T @ A))  # 1.38e+16   ← квадрат

x_ne  = np.linalg.solve(A.T @ A, A.T @ b)     # нормальные уравнения
x_svd = np.linalg.lstsq(A, b, rcond=None)[0]  # SVD
Q, R  = np.linalg.qr(A)
x_qr  = np.linalg.solve(R, Q.T @ b)           # QR

print(np.abs(x_ne  - x_true).max())   # 0.229    ← катастрофа
print(np.abs(x_svd - x_true).max())   # 6.6e-09
print(np.abs(x_qr  - x_true).max())   # 1.6e-08

Разница в восемь порядков на одной задаче. Правило, которое стоит вытатуировать: никогда не собирайте AᵀA руками — есть lstsq (SVD), QR и cho_solve для SPD.

Псевдообратная Мура–Пенроуза A⁺ = V Σ⁺ Uᵀ, где Σ⁺ обращает ненулевые σᵢ и оставляет нули нулями. Она даёт решение МНК всегда, а при бесконечном числе решений выбирает минимальное по норме, то есть неявно регуляризует:

A = np.array([[1.0, 2.0], [2.0, 4.0], [3.0, 6.0]])   # ранг 1: столбцы пропорциональны
b = np.array([1.0, 0.0, 1.0])
x = np.linalg.pinv(A) @ b
print(x, np.linalg.norm(x))     # [0.0571 0.1143] 0.1278 — минимальное по норме решение
print(np.linalg.matrix_rank(A), np.linalg.svd(A, compute_uv=False))   # 1  [8.367 0.]

Параметр rcond в pinv/lstsq — отсечка «какие σ считать нулями». Слишком мал — усиление шума делением на почти-ноль; слишком велик — потеря сигнала. Такое отсечение родственно гребневой регрессии: Ridge с параметром α заменяет 1/σ на σ/(σ² + α), гася вклад малых σ плавно, а не ступенькой.

Где SVD работает в проде

Важная оговорка: «SVD» в контексте Netflix Prize (Funk SVD) — не настоящее SVD. Настоящее требует полной матрицы, а матрица рейтингов на 99% пуста. То, что зовут SVD в рекомендациях, — градиентная минимизация Σ_{(u,i) наблюдённые} (r_ui − pᵤᵀqᵢ)² + λ(…), то есть малоранговая факторизация по наблюдённым ячейкам. Идея та же, алгоритм другой. Каноническое изложение: Koren, Bell, Volinsky, IEEE Computer 2009 (https://ieeexplore.ieee.org/document/5197422).

Часть 3. Рабочая свита: LU, QR, Холецкий, Шур

Спектральное разложение и SVD отвечают на вопрос «что матрица делает». Остальные — на вопрос «как быстро и устойчиво посчитать».

LU: как на самом деле решают Ax = b

PA = LU, где P — перестановка строк, L нижнетреугольная с единицами на диагонали, U верхнетреугольная. Это метод Гаусса в матричной записи: 2n³/3 операций, после чего каждое новое b решается за O(n²). Перестановка P (частичный выбор ведущего элемента) — не украшение: без неё деление на крошечный ведущий элемент уничтожает точность.

import scipy.linalg as sla

A = rng.standard_normal((6, 6))
lu, piv = sla.lu_factor(A)                # факторизуем ОДИН раз
b = rng.standard_normal(6)
x = sla.lu_solve((lu, piv), b)            # решаем МНОГО раз за O(n²)
print(np.abs(A @ x - b).max())            # 1.1e-14

Практика: решаете систему для многих правых частей — факторизуйте один раз. И никогда не вычисляйте inv(A) ради решения системы: это и дороже (2n³ против 2n³/3), и менее точно. Пишите np.linalg.solve(A, b), а не np.linalg.inv(A) @ b.

Холецкий: разложение для SPD

A = L Lᵀ, L нижнетреугольная с положительной диагональю. Существует ⟺ A симметрична и положительно определена. Вдвое быстрее LU (n³/3), не требует выбора ведущего, численно образцово устойчив.

B = rng.standard_normal((5, 5))
S = B @ B.T + 5 * np.eye(5)               # гарантированно SPD
L = np.linalg.cholesky(S)

# 1. Устойчивый логарифм определителя (без переполнения)
print(2 * np.log(np.diag(L)).sum(), np.linalg.slogdet(S)[1])    # совпадают

# 2. Сэмплирование из многомерного нормального распределения
z = rng.standard_normal((200_000, 5))     # N(0, I)
samples = z @ L.T                         # ковариация теперь ровно S

# 3. Проверка положительной определённости: исключение — это ответ
try:
    np.linalg.cholesky(np.array([[1.0, 2.0], [2.0, 1.0]]))
except np.linalg.LinAlgError as e:
    print("не SPD:", e)

Третий приём стоит подчеркнуть: самый дешёвый способ проверить SPD — попробовать разложить по Холецкому, это в разы быстрее полного спектра. В гауссовских процессах и калмановской фильтрации именно падение Холецкого сигнализирует, что ковариация «поехала» от накопленных ошибок округления и пора добавить jitter + εI.

QR и Шур

A = QR: Q с ортонормированными столбцами, R верхнетреугольная. Строится отражениями Хаусхолдера (стандарт LAPACK) или вращениями Гивенса (хорошо для разреженных и для обновлений); наивный Грам–Шмидт численно неустойчив. QR — рабочая лошадка МНК (Ax ≈ b ⟹ Rx = Qᵀb) и основа QR-алгоритма поиска спектра.

A = Q T Qᵀ (Шур): Q ортогональна, T квазиверхнетреугольная — с блоками 2×2 для пар комплексно-сопряжённых значений; собственные значения стоят на диагонали.

A = rng.standard_normal((5, 5))
T, Z = sla.schur(A)
print(np.abs(Z @ T @ Z.T - A).max())      # 6.2e-15
print(np.abs(np.tril(T, -2)).max())       # 0.0 — квазиверхнетреугольная

Разложение Шура существует всегда и вычисляется устойчиво, потому что Q ортогональна. Это и есть ответ на вопрос «что делать с дефектными матрицами на практике»: не Жордан, а Шур. Матричные функции (scipy.linalg.expm, funm) реализованы через него.

Сводная таблица и выбор разложения

Отсюда универсальный совет по производительности: берите самое слабое разложение, решающее вашу задачу. SVD знает про матрицу всё, но платить за него ~20× от стоимости LU ради одной системы — расточительство.

Часть 4. Как это на самом деле считают

Почему НЕ через характеристический многочлен

Учебник говорит: «найдите корни det(A − λI) = 0». Численно так делать нельзя никогда: корни многочлена катастрофически чувствительны к коэффициентам (классический пример Уилкинсона), да и переход матрица → коэффициенты сам теряет точность.

d = np.arange(1.0, 21.0)
A = np.diag(d)                                   # спектр очевиден: 1, 2, …, 20
c = np.poly(A)                                   # коэффициенты χ_A
roots = np.sort(np.real(np.roots(c)))
print(np.abs(roots - d).max())                   # 0.0698  ← ошибка ~7% на λ = 20
print(np.abs(np.linalg.eigvalsh(A) - d).max())   # 0.0     ← точно

Ошибка почти в десятую долю на диагональной матрице, где ответ написан прямо в элементах. Ирония в том, что реальные алгоритмы идут в обратную сторону: чтобы найти корни многочлена, numpy строит сопровождающую матрицу и считает её спектр QR-алгоритмом — np.roots устроен именно так.

Степенной метод

Псевдокод:
    x ← случайный единичный вектор
    повторять:
        y ← A x
        x ← y / ‖y‖
        λ ← xᵀ A x            (отношение Рэлея)
    пока λ не стабилизируется

Сходимость: ошибка ~ |λ₂/λ₁|ᵏ — линейная, со скоростью спектральной щели.
Стоимость итерации: O(nnz(A)) — только умножение матрицы на вектор.

def power_iteration(A, iters=1000, tol=1e-12, seed=0):
    """Старшее по модулю собственное значение и вектор. O(iters * nnz(A))."""
    rng = np.random.default_rng(seed)
    x = rng.standard_normal(A.shape[0])
    x /= np.linalg.norm(x)
    lam_old = 0.0
    for k in range(iters):
        y = A @ x
        x = y / np.linalg.norm(y)
        lam = x @ (A @ x)                    # отношение Рэлея
        if abs(lam - lam_old) < tol:
            return lam, x, k + 1
        lam_old = lam
    return lam_old, x, iters

S = rng.standard_normal((200, 200)); S = S + S.T
lam, v, it = power_iteration(S)
print(lam, np.linalg.eigvalsh(S)[-1])        # 40.14314  40.14314

Огромное достоинство: матрица нужна лишь как «чёрный ящик, умеющий умножать на вектор». Для графа в миллиард рёбер вы её никогда не материализуете — вы делаете проход по рёбрам. Ровно так считают PageRank в Spark.

Недостатки и лечение: даёт только одно значение → дефляция (вычесть λ₁u₁v₁ᵀ) или блочные варианты; медленно при |λ₂| ≈ |λ₁| → обратная итерация со сдвигом (степенной метод для (A − μI)⁻¹ даёт значение, ближайшее к μ, тем быстрее, чем точнее угадан сдвиг); не сходится, если старшая пара комплексно-сопряжённая.

QR-алгоритм: то, что внутри LAPACK

Промышленный стандарт с 1961 года (Фрэнсис, Кублановская). Идея обманчиво проста:

A₀ ← A   (сначала приводим к хессенберговой форме за O(n³) — это критично для скорости)
повторять:
    Qₖ Rₖ ← Aₖ            (QR-разложение)
    A_{k+1} ← Rₖ Qₖ       (перемножили в обратном порядке!)

A_{k+1} = Qₖᵀ Aₖ Qₖ — подобие, спектр сохраняется, а последовательность сходится к форме Шура. С неявными сдвигами Уилкинсона сходимость кубическая: 2–3 итерации на собственное значение. Для симметричных есть более быстрые специализации — divide-and-conquer (dsyevd) и MRRR (dsyevr). Для SVD аналог — алгоритм Голуба–Кахана: бидиагонализация Хаусхолдера, затем QR-итерации по бидиагональной матрице, без вычисления AᵀA на любом этапе.

Из цепочки следуют две практичные вещи. Первая: разница между gesdd (divide-and-conquer, быстрее, больше памяти) и gesvd (QR-итерации, надёжнее на патологических матрицах) — реальный тумблер scipy.linalg.svd(A, lapack_driver='gesvd'). Вторая: если np.linalg.svd кидает LinAlgError: SVD did not converge, первым делом ищите NaN/inf в матрице (в 95% случаев дело в них), вторым — пробуйте gesvd.

Крылов, Ланцош, Арнольди и randomized SVD

Когда матрица разрежена и имеет размер 10⁷×10⁷, полное разложение невозможно в принципе: множители U, V плотные и не влезут в память никогда. Работают подпространства Крылова span(x, Ax, A²x, …):

• Ланцош — для симметричных: трёхдиагональная проекция, крайние собственные значения за десятки умножений на вектор. Под капотом scipy.sparse.linalg.eigsh.
• Арнольди — то же для несимметричных, хессенбергова проекция. Под капотом eigs (обёртка над ARPACK).
• LOBPCG — очень большие SPD-задачи с предобуславливателем.

Отдельно — randomized SVD (Halko, Martinsson, Tropp, 2011, https://arxiv.org/abs/0909.4061), за 15 лет ставший стандартом приближённого малорангового разложения:

def randomized_svd(A, k, p=10, q=2, seed=0):
    """Приближённое SVD ранга k. O(mn(k+p)) вместо O(mn·min(m,n))."""
    rng = np.random.default_rng(seed)
    Omega = rng.standard_normal((A.shape[1], k + p))   # оверсэмплинг p ≈ 5..10
    Q, _ = np.linalg.qr(A @ Omega)                     # эскиз пространства столбцов
    for _ in range(q):                                 # степенные итерации: прижимаем хвост спектра
        Q, _ = np.linalg.qr(A.T @ Q)
        Q, _ = np.linalg.qr(A @ Q)
    B = Q.T @ A                                        # маленькая (k+p)×n
    Ub, s, Vt = np.linalg.svd(B, full_matrices=False)
    return (Q @ Ub)[:, :k], s[:k], Vt[:k]

Abig = (rng.standard_normal((2000, 30))                # спектр убывает экспоненциально
        @ np.diag(np.exp(-np.arange(30) / 5))
        @ rng.standard_normal((30, 1000)))

_, s_rand, _ = randomized_svd(Abig, 10)
s_exact = np.linalg.svd(Abig, compute_uv=False)[:10]
print(np.abs(s_rand - s_exact).max() / s_exact[0])     # 4e-09 — практически точно
# замер: randomized ≈ 0.014 c, полное SVD ≈ 0.73 c → ~50× быстрее

Пятидесятикратное ускорение при относительной ошибке 10⁻⁹ — не магия, а следствие убывания спектра: случайная проекция с высокой вероятностью «ловит» доминирующее подпространство, и Halko–Martinsson–Tropp дают на это строгие вероятностные оценки. Готовая реализация — sklearn.utils.extmath.randomized_svd.

Часть 5. Практика

Шпаргалка по API

import numpy as np, scipy.linalg as sla, scipy.sparse.linalg as spla

# --- спектр ---
np.linalg.eigvalsh(S)          # только λ, СИММЕТРИЧНАЯ: быстрее и точнее
np.linalg.eigh(S)              # λ и ортонормированные векторы, симметричная
np.linalg.eig(A)               # общая матрица: ждите комплексные λ
sla.eigh(S, B)                 # обобщённая задача S v = λ B v
spla.eigsh(S_sparse, k=6)      # k крайних λ разреженной симметричной (ARPACK)

# --- SVD ---
np.linalg.svd(A, full_matrices=False)   # экономичное: U m×r, Vt r×n
np.linalg.svd(A, compute_uv=False)      # только σ — заметно дешевле
sla.svd(A, lapack_driver='gesvd')       # надёжнее на патологических случаях
spla.svds(A_sparse, k=50)               # k старших для разреженной

# --- решение систем ---
np.linalg.solve(A, b)               # LU. НЕ inv(A) @ b
sla.lu_factor / sla.lu_solve        # многократные правые части
sla.cho_factor / sla.cho_solve      # SPD, вдвое быстрее
np.linalg.lstsq(A, b, rcond=None)   # МНК через SVD, работает при неполном ранге

# --- диагностика ---
np.linalg.cond(A)              # σ₁/σ_min: >1e8 — тревога, >1e15 — приговор для float64
np.linalg.matrix_rank(A)       # численный ранг через σ и порог
np.linalg.slogdet(A)           # log|det| без переполнения (det переполняется уже при n~200)

Три правила, экономящие время и нервы:

1. eigh вместо eig для симметричных. Вдвое-вчетверо быстрее, гарантированно вещественный спектр, ортогональные векторы. eig на симметричной матрице вернёт λ с мнимой частью 1e−17 и заставит писать np.real.
2. compute_uv=False, если нужны только σ. Для оценки ранга или cond векторы не нужны, а стоят они большую часть времени.
3. Порядок не гарантирован у eig; eigh возвращает по возрастанию, svd — по убыванию. Сортируйте явно, если логика зависит от порядка.

Спектральная кластеризация: собственные векторы графа

Собственные векторы лапласиана L = D − W знают о графе поразительно много: кратность нуля равна числу компонент связности, а второй снизу вектор (вектор Фидлера) даёт разрез, близкий к минимальному нормализованному.

W = np.zeros((6, 6))
for i, j in [(0,1),(1,2),(0,2),(3,4),(4,5),(3,5),(2,3)]:   # два треугольника + мостик
    W[i, j] = W[j, i] = 1.0

L = np.diag(W.sum(1)) - W
w, V = np.linalg.eigh(L)
print(w)                  # [0.  0.438  3.  3.  3.  4.562]
print(np.sign(V[:, 1]))   # [-1 -1 -1  1  1  1]  ← ровно два кластера

λ₁ = 0 всегда (вектор из единиц лежит в ядре). λ₂ = 0.438 — алгебраическая связность: чем она меньше, тем уже узкое место графа. Знаки вектора Фидлера безошибочно разделили треугольники, разрезав мостик. На этом построены sklearn.cluster.SpectralClustering, сегментация изображений (normalized cuts) и разбиение расчётных сеток. Подробнее — Теория графов.

Функции от матрицы

B = rng.standard_normal((4, 4)); S = B @ B.T          # PSD
w, V = np.linalg.eigh(S)

expS  = (V * np.exp(w))  @ V.T                        # матричная экспонента
sqrtS = (V * np.sqrt(w)) @ V.T                        # квадратный корень
print(np.abs(expS - sla.expm(S)).max())               # 9e-11
print(np.abs(sqrtS @ sqrtS - S).max())                # 3.6e-15

Идиома (V * f(w)) @ V.T — умножение на диагональ через broadcasting, без построения diag. Работает только для симметричных; в общем случае берите scipy.linalg.expm (масштабирование и Паде). Писать матричную экспоненту рядом Тейлора самостоятельно не надо — это классический источник катастрофической потери точности («Nineteen Dubious Ways to Compute the Exponential of a Matrix», Moler & Van Loan, https://doi.org/10.1137/S00361445024180).

Типичные заблуждения

«Собственные векторы всегда ортогональны». Только для симметричных (шире — нормальных, AAᵀ = AᵀA) матриц. У общей матрицы они бывают почти коллинеарны, тогда P плохо обусловлена, а вычисленный спектр ненадёжен.

«Собственные значения и сингулярные числа — почти одно и то же». Совпадают только для симметричных PSD (там σᵢ = |λᵢ|). В общем случае связи нет: у [[0,1],[0,0]] оба собственных значения нулевые, а сингулярные числа — 1 и 0.

«Большое собственное значение = важное». Для PCA — да, по построению (дисперсия). В некорректно поставленных обратных задачах интересное живёт как раз в малых σ, и там же весь шум.

«det(A) = 0 — значит матрица вырождена». Численно определитель почти бесполезен: у 0.1·I размера 200×200 определитель равен 10⁻²⁰⁰ (машинный ноль), хотя матрица идеально обусловлена. Судите по σ_min и cond.

«Симметрична — значит положительно определена». Нет: [[1,2],[2,1]] симметрична, а спектр {3, −1}. Проверка — Холецкий или eigvalsh.

«SVD слишком дорогое для продакшена». Полное SVD плотной 10⁴×10⁴ — да, десятки секунд. Но есть усечённое, разреженное (svds), randomized и инкрементальные варианты: вопрос не «можно ли», а «какой вариант».

«PCA убирает корреляции, значит признаки стали независимыми». PCA даёт некоррелированные компоненты (диагональная ковариация). Независимость сильнее; за ней идут к ICA.

Мини-итог

• Собственный вектор — направление, которое матрица только масштабирует; собственное значение — коэффициент. Найти их = найти естественную систему координат матрицы.
• tr(A) = Σλ, det(A) = Πλ — бесплатные проверки. ρ(A) < 1 — критерий сходимости итераций.
• A = PDP⁻¹ превращает степени и функции матрицы в операции над числами, но существует не всегда: дефектные матрицы. Жордан — теория, Шур — практика.
• Спектральная теорема: симметричные матрицы всегда диагонализуемы ортогональным преобразованием A = QΛQᵀ. Это фундамент PCA, оптимизации и всей статистики.
• SVD A = UΣVᵀ существует для любой матрицы и выдаёт разом ранг, нормы, обусловленность, четыре подпространства и оптимальное малоранговое приближение (Эккарт–Янг).
• PCA = SVD центрированных данных. LSA, рекомендации, LoRA, сжатие эмбеддингов — тот же приём в разных декорациях.
• Никогда не собирайте AᵀA и не ищите корни характеристического многочлена: обусловленность возводится в квадрат. Есть lstsq, QR, eigh.
• Берите самое дешёвое разложение, решающее задачу: Холецкий < LU < QR < спектральное < SVD.
• Для огромных данных: Ланцош/Арнольди (eigsh/eigs) для разреженных, randomized SVD для плотных.

Источники

• Gilbert Strang, «Introduction to Linear Algebra» и курс MIT 18.06 — https://ocw.mit.edu/courses/18-06-linear-algebra-spring-2010/
• Trefethen & Bau, «Numerical Linear Algebra», SIAM 1997 — эталон по устойчивости и обусловленности; лекции 4–5 про SVD и 12 про cond обязательны.
• Golub & Van Loan, «Matrix Computations», 4-е изд., Johns Hopkins 2013 — справочник по каждому алгоритму из этой статьи.
• Higham, «Accuracy and Stability of Numerical Algorithms», SIAM 2002 — https://epubs.siam.org/doi/book/10.1137/1.9780898718027
• LAPACK Users’ Guide — https://netlib.org/lapack/lug/ (какой driver что делает и когда какой брать).
• numpy.linalg — https://numpy.org/doc/stable/reference/routines.linalg.html ; scipy.linalg — https://docs.scipy.org/doc/scipy/reference/linalg.html
• Halko, Martinsson, Tropp, «Finding Structure with Randomness», SIAM Review 2011 — https://arxiv.org/abs/0909.4061
• Hu et al., «LoRA: Low-Rank Adaptation of Large Language Models» — https://arxiv.org/abs/2106.09685
• von Luxburg, «A Tutorial on Spectral Clustering» — https://arxiv.org/abs/0711.0189

Что дальше

Мы всё время пользовались тем, что скаляры образуют поле, векторы — абелеву группу по сложению, а обратимые матрицы — группу по умножению, но ни разу не спросили, что эти слова значат. Следующая статья трека вводит алгебраические структуры явно и показывает, где они уже живут в вашем коде: моноиды в reduce и CRDT, группы в криптографии, конечные поля в кодах Рида–Соломона и erasure-кодировании.

Абстрактная алгебра: группы, кольца, поля и моноиды в коде