Линейные системы — практический урок AI Engineering

Решение Ax = b — старейшая математическая задача, которая до сих пор работает внутри вашей нейронной сети.

Тип: Сборка Язык: Python Предварительные требования: Фаза 1, уроки 01 (Интуиция линейной алгебры), 02 (Векторы и матрицы), 03 (Преобразования матриц) Время: ~120 минут

Цели обучения

Решать Ax = b методом исключения Гаусса с частичным выбором главного элемента и обратной подстановкой
Разлагать матрицы на LU-, QR- и разложение Холецкого и объяснять, когда каждое из них уместно
Вывести нормальные уравнения для метода наименьших квадратов и связать их с линейной и гребневой регрессией
Диагностировать плохо обусловленные системы по числу обусловленности и применять регуляризацию для их стабилизации

Задача

Каждый раз, когда вы обучаете линейную регрессию, вы решаете линейную систему. Каждый раз, когда вы вычисляете аппроксимацию методом наименьших квадратов, вы решаете линейную систему. Каждый раз, когда слой нейронной сети вычисляет y = Wx + b, он вычисляет одну сторону линейной системы. Когда вы добавляете регуляризацию, вы изменяете систему. Когда вы используете гауссовские процессы, вы факторизуете матрицу. Когда вы обращаете ковариационную матрицу для расстояния Махаланобиса, вы решаете линейную систему.

Уравнение Ax = b встречается повсюду. A — матрица известных коэффициентов. b — вектор известных выходов. x — вектор неизвестных, который вы хотите найти. В линейной регрессии A — матрица данных, b — целевой вектор, а x — вектор весов. Вся модель сводится к следующему: найти x так, чтобы Ax было максимально близко к b.

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

Концепция

Геометрический смысл Ax = b

Система линейных уравнений имеет геометрическую интерпретацию. Каждое уравнение задаёт гиперплоскость. Решение — точка (или множество точек), в которой пересекаются все гиперплоскости.

2x + y = 5          Две прямые в 2D.
x - y  = 1          Они пересекаются в x=2, y=1.

Диаграмма к уроку «Линейные системы»

Возможны три исхода:

Диаграмма к уроку «Линейные системы»

В матричной форме «одно решение» означает, что A обратима. «Нет решения» означает, что система несовместна. «Бесконечно много решений» означает, что у A есть нулевое пространство. Большинство задач машинного обучения относятся к категории «нет точного решения», потому что у вас больше уравнений (точек данных), чем неизвестных (параметров). Здесь вступает в игру метод наименьших квадратов.

Представление по столбцам и по строкам

У Ax = b есть два способа прочтения.

Представление по строкам. Каждая строка A задаёт одно уравнение. Каждое уравнение — гиперплоскость. Решение находится там, где они все пересекаются.

Представление по столбцам. Каждый столбец A — вектор. Вопрос становится таким: какая линейная комбинация столбцов A даёт b?

A = | 2  1 |    b = | 5 |
    | 1 -1 |        | 1 |

Представление по строкам: решить одновременно 2x + y = 5 и x - y = 1.

Представление по столбцам: найти x1, x2 такие, что:
  x1 * [2, 1] + x2 * [1, -1] = [5, 1]
  2 * [2, 1] + 1 * [1, -1] = [4+1, 2-1] = [5, 1]   проверка.

Представление по столбцам более фундаментально. Если b лежит в пространстве столбцов A, система имеет решение. Если b там не лежит, вы находите ближайшую точку в пространстве столбцов. Эта ближайшая точка и есть решение методом наименьших квадратов.

Исключение Гаусса

Исключение Гаусса преобразует Ax = b в верхнетреугольную систему Ux = c, которую вы решаете обратной подстановкой. Это наиболее прямой метод.

Алгоритм:

1. Для каждого столбца k (столбца главного элемента):
   a. Найти наибольший элемент в столбце k в строке k или ниже неё (частичный выбор главного элемента).
   b. Поменять эту строку местами со строкой k.
   c. Для каждой строки i ниже k:
      - Вычислить множитель m = A[i][k] / A[k][k]
      - Вычесть из строки i строку k, умноженную на m.
2. Выполнить обратную подстановку: решать от последнего уравнения к первому.

Пример:

Исходная система:
| 2  1  1 | 8 |       R2 = R2 - (2)R1     | 2  1   1 |  8 |
| 4  3  3 |20 |  -->  R3 = R3 - (1)R1 --> | 0  1   1 |  4 |
| 2  3  1 |12 |                            | 0  2   0 |  4 |

                       R3 = R3 - (2)R2     | 2  1   1 |  8 |
                                       --> | 0  1   1 |  4 |
                                           | 0  0  -2 | -4 |

Обратная подстановка:
  -2 * x3 = -4    -->  x3 = 2
  x2 + 2  = 4     -->  x2 = 2
  2*x1 + 2 + 2 = 8 --> x1 = 2

Исключение Гаусса требует O(n^3) операций. Для системы 1000x1000 это примерно миллиард операций с плавающей точкой. Быстро, но можно сделать лучше, если вам нужно решать несколько систем с одной и той же A.

Частичный выбор главного элемента: почему он важен

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

Плохой главный элемент:             С частичным выбором главного элемента:
| 0.001  1 | 1.001 |            Сначала поменять строки местами:
| 1      1 | 2     |            | 1      1 | 2     |
                                 | 0.001  1 | 1.001 |
m = 1/0.001 = 1000              m = 0.001/1 = 0.001
R2 = R2 - 1000*R1               R2 = R2 - 0.001*R1
| 0.001  1     | 1.001   |      | 1      1     | 2     |
| 0     -999   | -999.0  |      | 0      0.999 | 0.999 |

x2 = 1.000 (верно)              x2 = 1.000 (верно)
x1 = (1.001 - 1)/0.001          x1 = (2 - 1)/1 = 1.000 (верно)
   = 0.001/0.001 = 1.000        Устойчиво, потому что множитель мал.

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

LU-разложение

LU-разложение факторизует A в нижнетреугольную матрицу L и верхнетреугольную матрицу U: A = LU. Матрица L хранит множители из исключения Гаусса. Матрица U — результат исключения.

A = L @ U

| 2  1  1 |   | 1  0  0 |   | 2  1   1 |
| 4  3  3 | = | 2  1  0 | @ | 0  1   1 |
| 2  3  1 |   | 1  2  1 |   | 0  0  -2 |

Зачем факторизовать, а не просто выполнить исключение? Потому что, получив L и U, вы можете решить Ax = b для любого нового b всего за O(n^2):

Ax = b
LUx = b
Пусть y = Ux:
  Ly = b    (прямая подстановка, O(n^2))
  Ux = y    (обратная подстановка, O(n^2))

Стоимость O(n^3) оплачивается один раз при факторизации. Каждое последующее решение требует O(n^2). Если вам нужно решить 1000 систем с одной и той же A, но разными векторами b, LU экономит множитель 1000/3 в общем объёме работы.

При частичном выборе главного элемента получается PA = LU, где P — матрица перестановки, записывающая обмены строк.

QR-разложение

QR-разложение факторизует A в ортогональную матрицу Q и верхнетреугольную матрицу R: A = QR.

Ортогональная матрица обладает свойством Q^T Q = I. Её столбцы — ортонормированные векторы. Умножение на Q сохраняет длины и углы.

A = Q @ R

Q имеет ортонормированные столбцы: Q^T Q = I
R верхнетреугольна

Чтобы решить Ax = b:
  QRx = b
  Rx = Q^T b    (просто умножить на Q^T, обращение не требуется)
  Выполнить обратную подстановку, чтобы получить x.

QR численно устойчивее LU для решения задач наименьших квадратов. Процесс Грама — Шмидта строит Q по столбцам:

Даны столбцы a1, a2, ... матрицы A:

q1 = a1 / ||a1||

q2 = a2 - (a2 . q1) * q1        (вычесть проекцию на q1)
q2 = q2 / ||q2||                (нормировать)

q3 = a3 - (a3 . q1) * q1 - (a3 . q2) * q2
q3 = q3 / ||q3||

R[i][j] = qi . aj    для i <= j

На каждом шаге удаляется компонента вдоль всех предыдущих векторов q, оставляя только новое ортогональное направление.

Разложение Холецкого

Когда A симметрична (A = A^T) и положительно определена (все собственные значения положительны), её можно факторизовать как A = L L^T, где L — нижнетреугольная матрица. Это разложение Холецкого.

A = L @ L^T

| 4  2 |   | 2  0 |   | 2  1 |
| 2  5 | = | 1  2 | @ | 0  2 |

L[i][i] = sqrt(A[i][i] - sum(L[i][k]^2 for k < i))
L[i][j] = (A[i][j] - sum(L[i][k]*L[j][k] for k < j)) / L[j][j]    for i > j

Холецкий вдвое быстрее LU и требует вдвое меньше памяти. Он работает только для симметричных положительно определённых матриц, но они встречаются постоянно:

Ковариационные матрицы симметричны и положительно полуопределены (а с регуляризацией — положительно определены).
Ядерная матрица в гауссовских процессах симметрична и положительно определена.
Гессиан выпуклой функции в точке минимума симметричен и положительно определён.
A^T A всегда симметрична и положительно полуопределена.

В гауссовских процессах вы факторизуете ядерную матрицу K с помощью Холецкого, затем решаете K alpha = y, чтобы получить предиктивное среднее. Множитель Холецкого также даёт логарифм определителя для маргинального правдоподобия: log det(K) = 2 * sum(log(diag(L))).

Метод наименьших квадратов: когда у Ax = b нет точного решения

Если A имеет размер m x n при m > n (уравнений больше, чем неизвестных), система является переопределённой. Точного решения нет. Вместо этого вы минимизируете квадрат ошибки:

[!warning] Уточнение переводчика Само по себе условие m > n не исключает точного решения: оно возможно, если все уравнения совместны. В типичной задаче с шумными данными точного решения действительно нет, поэтому применяется метод наименьших квадратов.

minimize ||Ax - b||^2

Это сумма квадратов невязок:
  sum((A[i,:] @ x - b[i])^2 for i in range(m))

Минимизатор удовлетворяет нормальным уравнениям:

A^T A x = A^T b

Вывод: раскроем ||Ax - b||^2 = (Ax - b)^T (Ax - b) = x^T A^T A x - 2 x^T A^T b + b^T b. Возьмём градиент по x и приравняем его к нулю: 2 A^T A x - 2 A^T b = 0.

Исходная система (переопределённая, 4 уравнения, 2 неизвестных):
| 1  1 |         | 3 |
| 1  2 | x     = | 5 |       Ни один точный x не удовлетворяет всем 4 уравнениям.
| 1  3 |         | 6 |
| 1  4 |         | 8 |

Нормальные уравнения:
A^T A = | 4  10 |    A^T b = | 22 |
        | 10 30 |            | 63 |

Решение: x = [1.5, 1.7]

Это линейная регрессия. x[0] — свободный член, x[1] — наклон.

Нормальные уравнения = линейная регрессия

Связь точна. В линейной регрессии ваша матрица данных X имеет одну строку на образец и один столбец на признак. Целевой вектор y имеет один элемент на образец. Вектор весов w удовлетворяет:

X^T X w = X^T y
w = (X^T X)^(-1) X^T y

Это замкнутая форма решения линейной регрессии. Каждый вызов sklearn.linear_model.LinearRegression.fit() вычисляет это (или эквивалентный вариант через QR или SVD).

Добавьте к матрице член регуляризации lambda * I — и получите гребневую регрессию:

(X^T X + lambda * I) w = X^T y
w = (X^T X + lambda * I)^(-1) X^T y

Регуляризация делает матрицу лучше обусловленной (её проще точно обратить) и предотвращает переобучение, сжимая веса к нулю. Матрица X^T X + lambda * I всегда симметрична и положительно определена при lambda > 0, поэтому для решения можно использовать Холецкого.

Псевдообратная матрица (Мура — Пенроуза)

Псевдообратная A+ обобщает обращение матрицы на неквадратные и вырожденные матрицы. Для любой матрицы A:

x = A+ b

где A+ = V Sigma+ U^T    (вычисляется через SVD)

Sigma+ образуется взятием обратной величины каждого ненулевого сингулярного значения и транспонированием результата. Если A = U Sigma V^T, то A+ = V Sigma+ U^T.

A = U Sigma V^T        (SVD)

Sigma = | 5  0 |       Sigma+ = | 1/5  0  0 |
        | 0  2 |                | 0  1/2  0 |
        | 0  0 |

A+ = V Sigma+ U^T

Псевдообратная даёт решение наименьших квадратов с минимальной нормой. Если система имеет:

Одно решение: A+ b его даёт.
Нет решения: A+ b даёт решение методом наименьших квадратов.
Бесконечно много решений: A+ b даёт решение с наименьшей ||x||.

np.linalg.lstsq и np.linalg.pinv из NumPy используют SVD внутри.

Число обусловленности

Число обусловленности измеряет, насколько решение чувствительно к малым изменениям входных данных. Для матрицы A число обусловленности равно:

kappa(A) = ||A|| * ||A^(-1)|| = sigma_max / sigma_min

где sigma_max и sigma_min — наибольшее и наименьшее сингулярные значения.

Хорошо обусловлена (kappa ~ 1):       Плохо обусловлена (kappa ~ 10^15):
Малое изменение b -->                 Малое изменение b -->
малое изменение x                     огромное изменение x

| 2  0 |   kappa = 2/1 = 2          | 1   1          |   kappa ~ 10^15
| 0  1 |   можно безопасно решать   | 1   1+10^(-15) |   решение — бессмыслица

Практические правила:

kappa < 100: безопасно, решение точно.
kappa ~ 10^k: вы теряете примерно k цифр точности из-за арифметики с плавающей точкой.
kappa ~ 10^16 (для float64): решение бессмысленно. Матрица фактически вырождена.

В машинном обучении плохая обусловленность возникает, когда признаки почти коллинеарны. Регуляризация (добавление lambda * I) улучшает число обусловленности с sigma_max / sigma_min до (sigma_max + lambda) / (sigma_min + lambda).

Итеративные методы: сопряжённые градиенты

Для очень больших разреженных систем (миллионы неизвестных) прямые методы, такие как LU или Холецкий, слишком дороги. Итеративные методы приближают решение, улучшая предположение на множестве итераций.

Метод сопряжённых градиентов (CG) решает Ax = b, когда A симметрична и положительно определена. Он находит точное решение не более чем за n итераций (в точной арифметике), но обычно сходится намного быстрее, если собственные значения сгруппированы.

Схема алгоритма:
  x0 = начальное предположение (часто ноль)
  r0 = b - A x0           (невязка)
  p0 = r0                 (направление поиска)

  Для k = 0, 1, 2, ...:
    alpha = (rk . rk) / (pk . A pk)
    x_{k+1} = xk + alpha * pk
    r_{k+1} = rk - alpha * A pk
    beta = (r_{k+1} . r_{k+1}) / (rk . rk)
    p_{k+1} = r_{k+1} + beta * pk
    если ||r_{k+1}|| < tolerance: остановиться

CG применяется в следующих случаях:

Крупномасштабная оптимизация (метод Newton-CG)
Решение дискретизаций уравнений в частных производных
Ядерные методы, где ядерная матрица слишком велика для факторизации
Предобуславливание для других итеративных решателей

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

Полная картина: какой метод когда

Метод	Требования	Стоимость	Сценарий использования
Исключение Гаусса	Квадратная, невырожденная A	O(n^3)	Однократное решение квадратной системы
LU-разложение	Квадратная, невырожденная A	O(n^3) факторизация + O(n^2) решение	Несколько решений с одной и той же A
QR-разложение	Любая A (m >= n)	O(mn^2)	Метод наименьших квадратов, численная устойчивость
Холецкий	Симметричная положительно определённая A	O(n^3/3)	Ковариационные матрицы, гауссовские процессы, гребневая регрессия
Нормальные уравнения	Переопределённая (m > n)	O(mn^2 + n^3)	Линейная регрессия (малое n)
SVD / псевдообратная	Любая A	O(mn^2)	Системы с дефицитом ранга, решения с минимальной нормой
Сопряжённые градиенты	Симметричная положительно определённая, разреженная A	O(n * k * nnz)	Большие разреженные системы, k = число итераций

Связь с машинным обучением

Каждый метод этого урока встречается в промышленном машинном обучении:

Линейная регрессия. Решение в замкнутой форме решает нормальные уравнения X^T X w = X^T y. Это делается через Холецкого (если n мало), QR (если важна численная устойчивость) или SVD (если матрица может иметь дефицит ранга).

Гребневая регрессия. Добавляет lambda * I к X^T X. Регуляризованная система (X^T X + lambda * I) w = X^T y всегда разрешима через Холецкого, потому что X^T X + lambda * I симметрична и положительно определена при lambda > 0.

Гауссовские процессы. Предиктивное среднее требует решения K alpha = y, где K — ядерная матрица. Стандартный подход — разложение K по Холецкому. Логарифм маргинального правдоподобия использует log det(K) = 2 sum(log(diag(L))).

Инициализация нейронной сети. Ортогональная инициализация использует QR-разложение для создания матриц весов, чьи столбцы ортонормированы. Это предотвращает затухание сигнала в глубоких сетях.

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

Конструирование признаков. Число обусловленности X^T X показывает, коллинеарны ли ваши признаки. Если kappa велико, удалите признаки или добавьте регуляризацию.

linear-system-conditioning

Соберите это

Шаг 1: исключение Гаусса с частичным выбором главного элемента

import numpy as np

def gaussian_elimination(A, b):
    n = len(b)
    Ab = np.hstack([A.astype(float), b.reshape(-1, 1).astype(float)])

    for k in range(n):
        max_row = k + np.argmax(np.abs(Ab[k:, k]))
        Abk, max_row = Abmax_row, k

        if abs(Ab[k, k]) < 1e-12:
            raise ValueError(f"Matrix is singular or nearly singular at pivot {k}")

        for i in range(k + 1, n):
            m = Ab[i, k] / Ab[k, k]
            Ab[i, k:] -= m * Ab[k, k:]

    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (Ab[i, -1] - Ab[i, i+1:n] @ x[i+1:n]) / Ab[i, i]

    return x

Шаг 2: LU-разложение

def lu_decompose(A):
    n = A.shape[0]
    L = np.eye(n)
    U = A.astype(float).copy()
    P = np.eye(n)

    for k in range(n):
        max_row = k + np.argmax(np.abs(U[k:, k]))
        if max_row != k:
            Uk, max_row = Umax_row, k
            Pk, max_row = Pmax_row, k
            if k > 0:
                L[[k, max_row], :k] = L[[max_row, k], :k]

        for i in range(k + 1, n):
            L[i, k] = U[i, k] / U[k, k]
            U[i, k:] -= L[i, k] * U[k, k:]

    return P, L, U

def lu_solve(P, L, U, b):
    n = len(b)
    Pb = P @ b.astype(float)

    y = np.zeros(n)
    for i in range(n):
        y[i] = Pb[i] - L[i, :i] @ y[:i]

    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]

    return x

Шаг 3: разложение Холецкого

def cholesky(A):
    n = A.shape[0]
    L = np.zeros_like(A, dtype=float)

    for i in range(n):
        for j in range(i + 1):
            s = A[i, j] - L[i, :j] @ L[j, :j]
            if i == j:
                if s <= 0:
                    raise ValueError("Matrix is not positive definite")
                L[i, j] = np.sqrt(s)
            else:
                L[i, j] = s / L[j, j]

    return L

Шаг 4: метод наименьших квадратов через нормальные уравнения

def least_squares_normal(A, b):
    AtA = A.T @ A
    Atb = A.T @ b
    return gaussian_elimination(AtA, Atb)

def ridge_regression(A, b, lam):
    n = A.shape[1]
    AtA = A.T @ A + lam * np.eye(n)
    Atb = A.T @ b
    L = cholesky(AtA)
    y = np.zeros(n)
    for i in range(n):
        y[i] = (Atb[i] - L[i, :i] @ y[:i]) / L[i, i]
    x = np.zeros(n)
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - L.T[i, i+1:] @ x[i+1:]) / L.T[i, i]
    return x

Шаг 5: число обусловленности

def condition_number(A):
    U, S, Vt = np.linalg.svd(A)
    return S[0] / S[-1]

Используйте это

Объединим все части для линейной и гребневой регрессии на реальных данных:

np.random.seed(42)
X_raw = np.random.randn(100, 3)
w_true = np.array([2.0, -1.0, 0.5])
y = X_raw @ w_true + np.random.randn(100) * 0.1

X = np.column_stack([np.ones(100), X_raw])

w_ols = least_squares_normal(X, y)
print(f"OLS weights (ours):    {w_ols}")

w_np = np.linalg.lstsq(X, y, rcond=None)[0]
print(f"OLS weights (numpy):   {w_np}")
print(f"Max difference: {np.max(np.abs(w_ols - w_np)):.2e}")

w_ridge = ridge_regression(X, y, lam=1.0)
print(f"Ridge weights (ours):  {w_ridge}")

from sklearn.linear_model import Ridge
ridge_sk = Ridge(alpha=1.0, fit_intercept=False)
ridge_sk.fit(X, y)
print(f"Ridge weights (sklearn): {ridge_sk.coef_}")

Подготовьте к выпуску

Этот урок создаёт:

code/linear_systems.py, содержащий реализации с нуля исключения Гаусса, LU-разложения, разложения Холецкого, метода наименьших квадратов и гребневой регрессии
Рабочую демонстрацию того, что нормальные уравнения и LinearRegression из sklearn дают одинаковые веса

Упражнения

Решите систему [[1,2,3],[4,5,6],[7,8,10]] x = [6, 15, 27] с помощью своего исключения Гаусса, своего LU-решателя и np.linalg.solve. Убедитесь, что все три дают один и тот же ответ с точностью до допуска вычислений с плавающей точкой.
Сгенерируйте случайную матрицу X размера 50x5 и целевой вектор y = X @ w_true + noise. Найдите w, используя нормальные уравнения, QR (через np.linalg.qr), SVD (через np.linalg.svd) и np.linalg.lstsq. Сравните все четыре решения. Измерьте число обусловленности X^T X и объясните, как оно влияет на выбор метода, которому вы доверяете.
Создайте почти вырожденную матрицу, сделав два столбца почти одинаковыми (например, столбец 2 = столбец 1 + 1e-10 * noise). Вычислите её число обусловленности. Решите Ax = b с регуляризацией и без неё (добавьте 0.01 * I). Сравните решения и невязки. Объясните, почему регуляризация помогает.
Реализуйте алгоритм сопряжённых градиентов для случайной симметричной положительно определённой матрицы 100x100. Подсчитайте, сколько итераций нужно для сходимости с допуском 1e-8. Сравните с теоретическим максимумом в n итераций.
Измерьте время работы своего решателя Холецкого, своего LU-решателя и np.linalg.solve на симметричных положительно определённых матрицах размеров 10, 50, 200, 500. Постройте график результатов. Убедитесь, что Холецкий примерно в 2 раза быстрее LU.

Ключевые термины

Термин	Как обычно говорят	Что это на самом деле означает
Линейная система	«Решить относительно x»	Набор линейных уравнений Ax = b. Найти x означает найти вход, который даёт выход b при преобразовании A.
Исключение Гаусса	«Привести по строкам»	Систематически обнулять элементы ниже диагонали строковыми операциями, получая верхнетреугольную систему, решаемую обратной подстановкой. O(n^3).
Частичный выбор главного элемента	«Менять строки для устойчивости»	Перед исключением в столбце k поменять строку с наибольшим абсолютным значением в этом столбце на позицию главного элемента. Предотвращает деление на малые числа.
LU-разложение	«Разложить на треугольники»	Записать A = LU, где L нижнетреугольная (хранит множители), а U верхнетреугольная (матрица после исключения). Распределяет стоимость O(n^3) на множество решений.
QR-разложение	«Ортогональная факторизация»	Записать A = QR, где Q имеет ортонормированные столбцы, а R верхнетреугольная. Для наименьших квадратов устойчивее LU.
Разложение Холецкого	«Квадратный корень матрицы»	Для симметричной положительно определённой A записать A = LL^T. Вдвое дешевле LU. Используется для ковариационных и ядерных матриц, а также гребневой регрессии.
Метод наименьших квадратов	«Наилучшая аппроксимация, когда точное решение невозможно»	Минимизировать сумму квадратов невязок
Нормальные уравнения	«Сокращение через исчисление»	A^T A x = A^T b. Получаются при приравнивании нулю градиента
Псевдообратная	«Обращение для неквадратных матриц»	A+ = V Sigma+ U^T через SVD. Даёт решение наименьших квадратов с минимальной нормой для любой матрицы — квадратной или прямоугольной, вырожденной или нет.
Число обусловленности	«Насколько этому ответу можно доверять»	kappa = sigma_max / sigma_min. Измеряет чувствительность к возмущениям входа. Теряется примерно log10(kappa) цифр точности.
Гребневая регрессия	«Регуляризованный метод наименьших квадратов»	Решает (X^T X + lambda I) w = X^T y. Добавление lambda I улучшает обусловленность и сжимает веса к нулю. Предотвращает переобучение.
Сопряжённые градиенты	«Итеративное Ax=b для больших матриц»	Итеративный решатель для симметричных положительно определённых систем. Сходится не более чем за n шагов. Практичен для больших разреженных систем, где факторизация слишком дорога.
Переопределённая система	«Данных больше, чем параметров»	m > n в системе размера m на n. Точного решения может не существовать. Метод наименьших квадратов находит наилучшее приближение. Это характерно для большинства задач регрессии.
Обратная подстановка	«Решать снизу вверх»	Для верхнетреугольной системы сначала решить последнее уравнение, затем подставлять назад. O(n^2).
Прямая подстановка	«Решать сверху вниз»	Для нижнетреугольной системы сначала решить первое уравнение, затем подставлять вперёд. O(n^2). Используется на шаге L при LU-решении.

Дополнительное чтение

MIT 18.06: Линейная алгебра (Gilbert Strang) – определяющий курс по линейным системам и факторизациям матриц
Численная линейная алгебра (Trefethen & Bau) – стандартный источник для понимания численной устойчивости, обусловленности и причин отказа алгоритмов
Вычисления с матрицами (Golub & Van Loan) – энциклопедический источник по всем алгоритмам для матриц
3Blue1Brown: Обратные матрицы – визуальная интуиция геометрического смысла решения Ax = b

Источник: оригинальный урок на GitHub Навигация: 01.16 — Методы сэмплирования · Фаза 1 — Математические основы · Полный каталог · 01.18 — Выпуклая оптимизация