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

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

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

2 ч теории. 120 минут, два перерыва, два live-coding блока.

Последняя лекция модуля: сегодня мы соберём из недель 1–3 первый ML-алгоритм — линейную регрессию.

Не открывая конспект: view или копия? Какой axis схлопывает строки?

Срез возвращает view или копию? Как получить настоящую копию?
m = np.arange(20).reshape(4, 5): какой shape у m.sum(axis=0)?
Что вернёт np.where(a > 3, a, 0) — массив той же формы или 1D-список значений?

Неделя 3 опирается на slicing (U[:, :k]) и broadcasting (U * S) — припомните приёмы недели 2.

СЛАУ, SVD-сжатие, Монте-Карло — и «регрессия руками» как первый шаг к ML

СЛАУ (solve) — основа Методик вычислений: любой численный метод сводится к решению линейной системы.
SVD — сжатие изображений (JPEG-подобная идея), PCA в Анализе данных и ML; «первое сингулярное число» — как «главный компонент».
Случайные числа + Монте-Карло — моделирование в статистике: оценка π, интегралы, имитация случайных процессов.
Проект модуля: «линейная регрессия руками» — нормальное уравнение — это первая ML-модель; «на ML-курсе это уже предполагается, а не проходится».

Какой из четырёх приёмов вы уже видите в своих будущих дисциплинах?

К концу лекции вы сможете…

различать * и @ и применять их по назначению;
решать СЛАУ через solve и проверять решение;
объяснять, почему не inv(A) @ b;
находить собственные значения/векторы и проверять A·v = λ·v;
сжимать изображение через усечённое SVD (k = 5, 20, 100);
генерировать воспроизводимые случайные данные и оценивать π методом Монте-Карло.

Самый эффектный результат недели — пункт 5: картинку увидите своими глазами.

`*` — поэлементно, `@` — матрично; путаница — самая частая ошибка недели

A = np.array([[3., 1.], [1., 2.]])
B = np.array([[1., 2.], [3., 1.]])
A * B        # поэлементно
array([[3., 2.],
       [3., 2.]])

A = np.array([[3., 1.], [1., 2.]])
B = np.array([[1., 2.], [3., 1.]])
A @ B        # матрично
array([[6., 7.],
       [7., 4.]])

Правило: * — каждый элемент с соответствующим (broadcasting); @ — строка × столбец, размерности согласуются. Путаница между ними — самая частая ошибка недели 3.

Откуда берутся 6 (3·1+1·3) и 7 (3·2+1·1) в A@B?

Вычислите вручную: `A@B` и `A*B` для малой пары матриц

A = np.array([[1., 2.], [3., 4.]])
B = np.array([[0., 1.], [1., 0.]])
A * B     # ?
A @ B     # ?

Не запуская: чему равны оба произведения? Посчитайте на бумаге — те же исходные матрицы, два разных результата.

`trace`/`diag`/`eye` — след, диагональ, единичная матрица

np.trace(A)      # след: сумма диагонали
np.diag(A)       # диагональ
np.eye(3)        # единичная матрица 3×3
5.0
array([3., 2.])
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

Правило: след понадобится для проверки собственных значений (Σλ = trace); eye — «нейтральный элемент» матричного умножения: A @ np.eye(n) = A.

Чему равен след у np.eye(3)?

`np.linalg.solve(A, b)` решает `Ax=b`; проверка — `A@x ≈ b`

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

Правило: система 3x + y = 9, x + 2y = 8 → x=2, y=3. Проверка подстановкой — всегда, одной строкой: A @ x (через @, а не *).

Почему проверка именно через @, а не через *?

Решение СЛАУ — через `solve`, а не `inv(A)@b`: устойчивость и скорость

# так НЕ надо
x = np.linalg.inv(A) @ b

# так надо
x = np.linalg.solve(A, b)

Объяснение: solve использует LU-разложение — численно устойчиво и быстрее; inv строит обратную матрицу целиком (лишние операции) и накапливает ошибки округления. Для решения СЛАУ — solve, inv оставьте для анализа.

Почему inv(A) @ b хуже, чем solve(A, b)? На собеседованиях этот вопрос — маркер.

`inv`/`det`/`cond` — когда нужны: анализ вырожденности и устойчивости

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

Правило: det ≈ 0 → матрица близка к вырожденной (решение неустойчиво); cond — число обусловленности: чем больше, тем чувствительнее решение к ошибкам входа (у хорошо обусловленной системы — единицы); inv нужен для анализа, не для решения.

Как оценить «устойчивость» системы? О чём говорит cond ≈ 2.6?

Стоп-вопрос блока 1

Как устойчиво решить СЛАУ и почему не через inv? Что такое число обусловленности?

Ответ-эталон: solve (LU), inv неустойчив и медленнее; cond — мера чувствительности решения к ошибкам входа.

`eig`: собственный вектор не меняет направление — `A·v = λ·v`

Собственный вектор v при умножении на A не меняет направление — только масштабируется на λ (собственное значение).

vals, vecs = np.linalg.eig(A)
vals
array([3.61803399, 1.38196601])
vecs
array([[ 0.85065081, -0.52573111],
       [ 0.52573111,  0.85065081]])

Проверка: A @ vecs[:, 0] ≈ vals[0] * vecs[:, 0]; контроль: vals.sum() == trace(A), vals.prod() == det(A). Векторы — столбцы vecs[:, i], не строки.

Что такое собственное значение геометрически? (коэффициент растяжения вдоль собственного вектора).

`eigh` для симметричных; на несимметричных `eig` даёт комплексные — это нормально

np.linalg.eigvals(np.array([[0., -1.], [1., 0.]]))
array([0.+1.j, 0.-1.j])   # комплексная пара — поворот на 90°

Правило: на несимметричной матрице eig может вернуть комплексные значения — это не ошибка, а свойство матрицы; для симметричных матриц используйте eigh — быстрее и гарантированно вещественные значения.

Честно о сложности: увидеть комплексные числа в ответе и испугаться — нормально; теперь вы знаете, что это значит.

Почему у поворота на 90° нет вещественной оси «растяжения»? Отсюда мнимость.

Сначала почувствуй: картинка при k=5, 20, 100

исходник 200×200

k=5 · mse=995.18

k=20 · mse=678.50

k=100 · mse=45.35

Как 200×200 чисел можно хранить компактнее, почти не потеряв качество? Что вы видите — чем k=100 лучше k=5?

SVD: `A = UΣVᵀ`; `U`, `S`, `Vt` — формы факторов

Любая матрица раскладывается A = U Σ Vᵀ, где Σ — диагональ из сингулярных чисел (в NumPy — вектор S), U и Vt — ортонормированные.

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

Правило: S — ВЕКТОР длины 200, а не матрица: сингулярные числа на «диагонали». Для восстановления по k компонентам берём U[:, :k], S[:k], Vt[:k, :].

Частая ошибка — np.diag(S[:k]). Почему это лишнее, если S уже вектор?

Сингулярные числа быстро убывают — первые компоненты несут главную информацию

S[:8]
array([42791.6,  5507.3,  1097.2,  1040.6,  1017.2,  1002.5,   998.5,    978.8, ...])

Сингулярные числа изображения: лог-шкала по оси y, ось x — номер компоненты. Резкий спад и «хвост».

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

Первое число — 42791, восьмое — уже 979. Ошибка от отбрасывания хвоста мала?

Усечённое SVD: восстановление по k компонентам, MSE падает с k

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

Правило: восстановление — (U[:, :k] * S[:k]) @ Vt[:k, :] (broadcasting недели 2!); ошибка падает с k; экономия хранения — k×(200+200) вместо 200×200 (при k=100 — ~2× сжатие почти без потерь).

Что делает U[:, :k] * S[:k] — как здесь работает broadcasting из недели 2?

Predict: почему MSE монотонно падает с ростом k?

В таблице MSE монотонно убывает с k: 995 → 678 → 298 → 45. Почему так? Может ли MSE вырасти при увеличении k?

С ростом k мы добавляем направления, которые были в исходнике — приближение только улучшается; MSE не может вырасти (это свойство SVD — наилучшее приближение ранга k).

Стоп-вопросы блока 2

Что такое сингулярные числа и как их используют для сжатия?
Как восстановить матрицу по k компонентам (формула)?
Почему MSE не может вырасти при росте k?

Если сомневаетесь в ответе — вернёмся к слайдам 16–17: «хвост» и формула (U[:, :k] * S[:k]) @ Vt[:k, :].

Live-coding: SVD-сжатие изображения

Задача: сгенерировать «изображение» 200×200 (градиент + шум), разложить через SVD и восстановить по k = 5, 20, 100 — сверьте MSE с таблицей.

# ожидаемый результат (сверьте с таблицей)
k=5:   mse=995.18
k=20:  mse=678.50
k=100: mse=45.35

Зачем в генерации картинки default_rng(0)? — скоро в блоке случайных чисел.

`default_rng(seed)` — современный генератор; seed фиксирует поток

rng = np.random.default_rng(42)
rng.random(3)
array([0.77395605, 0.43887844, 0.85859792])

Правило: default_rng(seed) — рекомендуемый генератор (вместо глобального np.random.seed); seed фиксирует начальное состояние → одинаковые последовательности при каждом запуске.

Зачем фиксировать seed? (воспроизводимость эксперимента — без неё тесты и ревью не повторят ваш результат).

Один seed → один поток: эксперимент воспроизводим; без seed — «плавает»

rng2 = np.random.default_rng(42)
rng2.random(3)
array([0.77395605, 0.43887844, 0.85859792])

Правило: два генератора с одним seed дают один поток; без seed (или default_rng() без аргумента) — новый поток при каждом запуске, результаты «плавают».

Совпали ли первые три числа с предыдущим слайдом? Воспроизводимость — требование проекта модуля.

Функции rng: random, normal, integers, choice, shuffle

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)             # перемешивание

np.random.default_rng(1).integers(1, 7, size=5)
[3 4 5 6 1]

Правило: random — равномерно [0,1); normal(μ,σ) / standard_normal — нормальное; integers — целые; choice(x, p=веса) — выбор с вероятностями; shuffle — перемешивание на месте.

Как сделать «нечестную монету»? — choice с весами.

Predict: два вызова `rng.random(3)` подряд — одинаковы?

rng = np.random.default_rng(42)
print(rng.random(3))
print(rng.random(3))

Одинаковы ли два вывода? Вспомните слайд 22: два генератора с одним seed — да; а два вызова одного генератора подряд?

Live-coding: нечестная монета и шум к данным

Задача: 1000 бросков нечестной монеты (choice с весами p=[0.3, 0.7]) и шум normal для синтетических данных регрессии (проект).

# ожидаемый результат
доля «р» ≈ 0.7    # в пределах случайной ошибки ~1/√N
noise.shape = (200,)   # шум к данным регрессии

Почему доля «р» не ровно 0.7? (случайность, сходимость ~1/√N — мостик к блоку Монте-Карло).

Монте-Карло: оценка π одной функцией без циклов

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

Точки в квадрате [0,1]², четверть круга. Доля точек с x²+y² ≤ 1 — это π/4.

Правило: доля точек с x²+y² ≤ 1 — площадь квадранта круга (π/4); × 4 → π. Векторизация недели 1: один вызов rng.random((n, 2)) и sum(axis=1) вместо циклов.

Почему оценка 3.1430 «почти π» — какой осью здесь снова пользуемся из недели 2?

Точность Монте-Карло растёт как ~1/√N

N	оценка π	ошибка
10³	3.0840	~0.06
10⁴	3.1600	~0.02
10⁵	3.1501	~0.008
10⁶	3.1430	~0.001

Правило: ошибка ~ 1/√N (ЦПТ): чтобы добавить один знак, нужно в 100 раз больше точек. Поэтому Монте-Карло — метод для «черновой» оценки, а не для высокой точности.

Почему точность растёт как 1/√N? В 100 раз больше точек — во сколько раз меньше ошибка?

Интеграл 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

Прямоугольник [0,1]², кривая y=x². Точки «под» кривой — доля площади = интеграл.

Правило: доля точек «под кривой» — площадь под графиком = интеграл; сравнение с точным 1/3 через допуск, не == (плавающая арифметика).

Как оценить ∫₀¹ x² dx? Одинакова ли идея с оценкой π?

Predict: `10³` или `10⁶` — какой результат ближе к 1/3?

mc_integral(10**3) или mc_integral(10**6) — какой результат ближе к 1/3? Во сколько раз точнее?

Ответ: 10⁶ — в √(10³) ≈ 32 раза точнее (ошибка ~1/√N). Больше точек — меньше ошибка, но цена растёт.

`X.T * X` вместо `X.T @ X` → `ValueError`; `*` поэлементно

X = np.random.default_rng(0).standard_normal((5, 3))
X.T * X        # хотели матричное произведение XᵀX
ValueError: operands could not be broadcast together
with shapes (3,5) (5,3)

Разбор: * требует поэлементной совместимости форм — (3,5) и (5,3) несовместимы (неделя 2!). Правильно: X.T @ X. Если бы X был квадратным — ошибки не было бы вовсе: * молча дал бы неверный результат. Поэтому @ для матриц — всегда осознанно.

Что говорит последняя строка traceback? Почему * не сработал?

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

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

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

Проект: нормальное уравнение `θ = (XᵀX)⁻¹Xᵀy` — всё сходится

rng = np.random.default_rng(0)
n = 200
x = rng.uniform(0, 10, n)                # размер жилья
y = 3.0 * x + 5.0 + rng.normal(0, 2, n)  # цена с шумом
X = np.column_stack([np.ones(n), x])     # + столбец единиц (bias)
theta = np.linalg.solve(X.T @ X, X.T @ y) # нормальное уравнение
theta
array([5.09, 2.95])   # ≈ [5.0, 3.0] — восстановленные параметры

Правило: нормальное уравнение θ = (XᵀX)⁻¹Xᵀy решается через solve, а НЕ inv (слайд 9); здесь работают @ (слайд 5), broadcasting и axis из недель 1–2; seed делает эксперимент воспроизводимым (слайд 21). Сравнение со sklearn — в проекте.

Шум rng.normal(0, 2, n) — тот же приём, что в слайде 25. Восстановились ли параметры ≈ 3.0 и 5.0?

Сегодня вы научились…

различать * и @;
решать СЛАУ через solve и проверять A@x ≈ b; оценивать cond;
находить собственные значения и понимать A·v = λ·v;
сжимать изображение через усечённое SVD (k=5/20/100);
делать эксперименты воспроизводимыми (default_rng(seed));
оценивать π и интегралы методом Монте-Карло (ошибка ~1/√N).

One-minute paper: главное, что поняли сегодня + один вопрос. Все шесть навыков нужны для проекта модуля — он у вас почти готов.

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

Семинар недели 3: «СЛАУ и проверка» (базовый), «Монте-Карло: π и интеграл» (стандартный), «SVD-сжатие + pinv через SVD» (продвинутый). Ноутбук w03/seminar-03.ipynb, ветка seminar/w03.
Проект модуля (сдаётся на неделе 3): линейная регрессия на NumPy без scikit-learn — linear_regression.py с fit_normal, fit_gd, predict, mse; сравнение со sklearn (MSE ~1e-2); тесты ≥5; Git + CI.
Этапы: реализация (слайд 32 — каркас), тесты, сравнение, оформление (README + CI). Рубрика 0–12 баллов — в 01-модуль-1.md.
Финал: «Модуль закрыт: от массивов до первой ML-модели. На следующей неделе — Pandas».

Готовы? Проект модуля — это не «ещё одно ДЗ», а ваш первый шаг в ML.