Семестр 3 · Модуль 1 · Неделя 3

Линейная алгебра и случайные числа

Матричные операции · solve · SVD и сжатие · воспроизводимые случайные числа · Монте-Карло

Мотивация недели: «линейная регрессия руками» (проект модуля) — первый шаг к ML. SVD сжимает изображения, Монте-Карло моделирует случайные процессы.

Матричные операции: @ vs *

A = np.array([[3., 1.], [1., 2.]])
B = np.array([[1., 2.], [3., 1.]])
A * B        # поэлементное произведение (Hadamard)
A @ B        # матричное умножение
A.T          # транспонирование
np.trace(A)  # след (сумма диагонали) → 5.0
np.diag(A)   # диагональ → [3. 2.]

Акцент: * — поэлементное, @ — матричное. Путаница между ними — самая частая ошибка недели 3.

Какой оператор даст «обычное» перемножение матриц?

Решение СЛАУ: np.linalg.solve

A = np.array([[3., 1.], [1., 2.]])
b = np.array([9., 8.])
x = np.linalg.solve(A, b)
x
array([2., 3.])
A @ x                    # проверка: должно быть ≈ b
array([9., 8.])

Предупреждение: НЕ решайте через np.linalg.inv(A) @ b — численно неустойчиво и медленнее. solve использует разложение LU — устойчиво.

Почему inv(A) @ b хуже, чем solve(A, b)?

inv, det — когда нужны, когда нет

np.linalg.inv(A)
array([[ 0.4, -0.2],
       [-0.2,  0.6]])
np.linalg.det(A)
5.0
np.linalg.cond(A)
2.618033988749896

det — оценка вырожденности: det ≈ 0 → матрица близка к вырожденной.
cond — число обусловленности: чем больше, тем неустойчивее решение.
Для решения СЛАУ — используйте solve, не inv.

Как оценить «устойчивость» системы линейных уравнений?

Собственные значения: np.linalg.eig

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

# A = [[3, 1], [1, 2]] (из слайда solve)
vals, vecs = np.linalg.eig(A)
vals
array([3.61803399, 1.38196601])
vecs
array([[ 0.85065081, -0.52573111],
       [ 0.52573111,  0.85065081]])

Для симметричных матриц — eigh (быстрее, гарантированно вещественные).

Что такое собственное значение геометрически?

SVD: A = U Σ Vᵀ

U, S, Vt = np.linalg.svd(img)
U.shape, S.shape, Vt.shape
((200, 200), (200,), (200, 200))

Сингулярные числа S — «важность» каждого направления матрицы; они быстро убывают.

# сингулярные числа реального «изображения»:
[42791.6, 5507.3, 1097.2, 1040.6, 1017.2, 1002.5, 998.5, 978.8, ...]

Первые числа огромны, дальше — «хвост»: значит, большую часть информации несут первые компоненты.

Что такое сингулярные числа и как они упорядочены?

Усечённое SVD — сжатие изображения

def compress(a, k):
    U, S, Vt = np.linalg.svd(a)
    return (U[:, :k] * S[:k]) @ Vt[:k, :]

for k in [5, 20, 50, 100]:
    rec = compress(img, k)
    mse = np.mean((img - rec)**2)
    print(f"k={k}: mse={mse:.2f}")

k=5:   mse=995.18
k=20:  mse=678.50
k=50:  mse=297.88
k=100: mse=45.35

Восстановление по k компонентам: качество растёт с k, но «вес» хранение — k × (200+200) вместо 200×200. При k=100 это уже 2× сжатие почти без потерь.

Почему ошибка восстановления падает с ростом k?

Воспроизводимые случайные числа

rng = np.random.default_rng(42)   # современный API
rng.random(3)
array([0.77395605, 0.43887844, 0.85859792])

Зачем seed: при фиксированном seed последовательность одинакова каждый раз — эксперимент воспроизводим.

rng2 = np.random.default_rng(42)
rng2.random(3)
array([0.77395605, 0.43887844, 0.85859792])  # тот же поток

Рекомендуется default_rng, а не глобальный np.random.seed. Без seed — новый поток при каждом запуске.

Как сделать случайный эксперимент воспроизводимым?

Функции rng

rng.random((2, 3))            # uniform [0, 1)
rng.standard_normal(5)      # N(0, 1)
rng.normal(5.0, 2.0, 3)     # N(μ=5, σ=2)
rng.integers(1, 7, size=5)  # целые 1..6
rng.choice(['о', 'р'], size=5, p=[0.3, 0.7])   # с весами
rng.shuffle(deck)           # перемешивание

# пример: integers(1, 7, 5) при seed=1
[3 4 5 6 1]

Чем rng.choice(x, p=...) отличается от rng.integers?

Монте-Карло: оценка π

Бросаем точки в квадрат [−1,1]²; доля попавших в круг × 4 → π.

def estimate_pi(n, seed=42):
    rng = np.random.default_rng(seed)
    pts = rng.random((n, 2))            # n точек в [0,1]²
    inside = np.sum(pts**2, axis=1) <= 1.0
    return 4 * inside.mean()

estimate_pi(10**6)
3.1430

Точность растёт как ~1/√N: чтобы получить ещё один знак, нужно в 100 раз больше точек.

Почему точность Монте-Карло растёт как 1/√N?

Монте-Карло: интеграл hit-and-miss

def mc_integral(n, seed=42):
    rng = np.random.default_rng(seed)
    pts = rng.random((n, 2))            # x, y в [0,1]
    return np.mean(pts[:, 1] <= pts[:, 0]**2)

mc_integral(10**6)
0.3327    # точное ∫₀¹ x² dx = 1/3 ≈ 0.3333

Доля точек «под кривой» — оценка площади = интеграла.

Как оценить ∫₀¹ x² dx методом Монте-Карло?

Типичные ошибки недели 3

A * B вместо A @ B — поэлементное вместо матричного.
Решение СЛАУ через inv — численно неустойчиво.
Нет seed — случайный эксперимент невоспроизводим.
eig на несимметричной матрице даёт комплексные значения — пугаются, хотя это нормально.
Путаница глобального np.random и default_rng.

Какую функцию использовать для устойчивого решения СЛАУ?

Вопросы для проверки понимания

Как устойчиво решить СЛАУ и почему не через inv?
Что такое сингулярные числа в SVD и как их используют для сжатия?
Как сделать случайный эксперимент воспроизводимым?
Как оценить π методом Монте-Карло?

Переход: на семинаре — решение СЛАУ, SVD-сжатие и Монте-Карло. А дальше — проект модуля: линейная регрессия на NumPy.

Что дальше: проект модуля

Семинар недели 3: «СЛАУ и проверка» (базовый), «Монте-Карло: π и интеграл» (стандартный), «SVD-сжатие + pinv через SVD» (продвинутый).
Проект модуля (нед. 3): линейная регрессия на NumPy без scikit-learn — нормальное уравнение или градиентный спуск, сравнение с sklearn, тесты, CI.
К концу недели вы соберёте из элементов модуля полноценную регрессию — ваш первый шаг в ML.

Готовы? Начните с задачи 3.1 — решение системы и проверка подстановкой.