- Переписывать код с циклами в векторизованные выражения NumPy и объяснять разницу в скорости
- Заранее определять форму результата операции над двумя массивами по правилам broadcasting
- Применять
@,np.linalg.solve,det,invиeigh, понимая их математический смысл
Мурад обрабатывает миллион продаж интернет-магазина: для каждой нужно посчитать сумму, применить скидку, найти изменение за день. Циклом for это занимает секунды, а в NumPy — миллисекунды. Секрет в векторизации: операцию записывают один раз для всего массива, а не для каждого элемента. В этом уроке мы также поработаем с линейной алгеброй — математическим языком науки о данных.
Векторизация: вычисления без циклов
Замена цикла Python по элементам операциями над целыми массивами. Цикл никуда не исчезает, но выполняется внутри NumPy в скомпилированном коде на C. Такие поэлементные функции называют ufunc (универсальные функции): +, *, np.sqrt, np.exp, np.where и т. д.
Почему цикл Python медленный? На каждой итерации интерпретатор читает инструкцию, проверяет тип элемента и создаёт новый объект float для результата. NumPy же один раз проверяет тип для чисел одного типа, лежащих в памяти подряд, и использует векторные (SIMD) инструкции процессора. Запусти код и убедись в разнице сам:
import time
import numpy as np
x = np.random.default_rng(0).random(1_000_000)
t0 = time.perf_counter()
total = 0.0
for v in x:
total += v * v
t_loop = time.perf_counter() - t0
t0 = time.perf_counter()
total_vec = np.dot(x, x)
t_vec = time.perf_counter() - t0
print(f'loop: {t_loop:.4f} s')
print(f'vectorized: {t_vec:.4f} s')
print('same result:', bool(np.isclose(total, total_vec)))
print(f'speed-up: about {t_loop / t_vec:.0f} times')| С циклом | Векторизованно |
|---|---|
| накопление суммы в переменной | a.sum(), a @ b |
if / else для каждого элемента | np.where(cond, x, y) |
| подсчёт элементов по условию | (a > 0).sum() |
| нарастающий итог | np.cumsum(a) |
| разности соседних элементов | np.diff(a) |
| максимум и его позиция | a.max(), a.argmax() |
import numpy as np
prices = np.array([12.0, 45.5, 8.9, 120.0, 64.0])
qty = np.array([3, 1, 10, 1, 2])
revenue = prices * qty
discount = np.where(revenue > 100, 0.10, 0.0)
to_pay = (revenue * (1 - discount)).round(2)
print(revenue)
print(to_pay, to_pay.sum())
daily = np.array([520.0, 610.0, 580.0, 700.0, 655.0])
print(np.diff(daily))
print((np.diff(daily) / daily[:-1] * 100).round(1))
print(np.cumsum(daily))▸ Ожидаемый результат
[ 36. 45.5 89. 120. 128. ] [ 36. 45.5 89. 108. 115.2] 393.7 [ 90. -30. 120. -45.] [17.3 -4.9 20.7 -6.4] [ 520. 1130. 1710. 2410. 3065.]
Broadcasting: массивы разной формы
Когда пишешь a * 10, число 10 как бы «растягивается» до массива той же формы, что и a. Этот механизм называется broadcasting (транслирование), и он работает не только с числами, но и с массивами разной формы. NumPy сравнивает формы справа налево:
- Если у одной формы меньше измерений, слева к ней дописываются единицы: (4,) → (1, 4).
- По каждой оси две длины должны быть равны, либо одна из них должна быть 1.
- Ось длины 1 «растягивается» до длины другого массива (без копирования памяти); форма результата — максимум по каждой оси.
- Если по какой-то оси длины разные и ни одна не равна 1 —
ValueError.
| A | B | Результат |
|---|---|---|
| (3, 4) | (4,) | (3, 4) |
| (3, 1) | (1, 4) | (3, 4) |
| (2, 1, 5) | (3, 1) | (2, 3, 5) |
| (3, 4) | (3,) | ошибка: 4 ≠ 3 |
Классическое применение broadcasting — стандартизация: из каждого столбца вычитаем его среднее и делим на стандартное отклонение, чтобы признаки разного масштаба (рост в см, вес в кг) стали сопоставимы. Из массива формы (4, 2) вычитается среднее формы (2,) — работают правила 1 и 3.
- xᵢⱼj-й признак i-го примера
- μⱼ, σⱼсреднее и стандартное отклонение j-го столбца
- zᵢⱼстандартизованное значение (z-оценка): среднее 0, стандартное отклонение 1
import numpy as np
X = np.array([[170.0, 65.0],
[160.0, 52.0],
[182.0, 80.0],
[175.0, 71.0]])
mean = X.mean(axis=0)
std = X.std(axis=0)
Z = (X - mean) / std
print(mean, std.round(3))
print(Z.round(2))
print(np.allclose(Z.mean(axis=0), 0), Z.std(axis=0).round(6))▸ Ожидаемый результат
[171.75 67. ] [ 8.012 10.173] [[-0.22 -0.2 ] [-1.47 -1.47] [ 1.28 1.28] [ 0.41 0.39]] True [1. 1.]
Если объединить вектор-столбец с вектором-строкой, получится таблица «каждый с каждым». a[:, np.newaxis] (или a[:, None]) превращает массив формы (3,) в столбец (3, 1):
import numpy as np
a = np.arange(1, 4)
b = np.arange(1, 6)
table = a[:, np.newaxis] * b
print(table.shape)
print(table)
M = np.ones((3, 4))
v = np.array([10.0, 20.0, 30.0])
try:
M + v
except ValueError as e:
print('Error:', str(e).strip())
print((M + v[:, None])[:, 0])▸ Ожидаемый результат
(3, 5) [[ 1 2 3 4 5] [ 2 4 6 8 10] [ 3 6 9 12 15]] Error: operands could not be broadcast together with shapes (3,4) (3,) [11. 21. 31.]
Найди форму результата каждой операции или покажи, что будет ошибка: а) (5, 1, 3) + (4, 3); б) (8, 1) · (1, 6); в) (2, 3) + (2,); г) (256, 256, 3) · (3,).
Показать решениеСкрыть решение
б) 8 и 1 → 8, 1 и 6 → 6. Результат (8, 6) — таблица «каждый с каждым».
в) Справа: 3 и 2 — не равны и ни одна не равна 1 → ValueError.
г) Последние оси 3 и 3 совпадают: результат (256, 256, 3). Так каналы R, G, B каждого пикселя изображения умножаются на три разных коэффициента.
Изменение формы
reshape показывает те же элементы в другой форме; число элементов не должно меняться, а одно из измерений можно задать как -1 — NumPy вычислит его сам. .T транспонирует, ravel() делает массив одномерным (по возможности возвращая представление), а np.vstack, np.hstack и np.concatenate склеивают массивы.
import numpy as np
a = np.arange(12)
print(a.reshape(3, 4))
print(a.reshape(2, -1).shape)
m = a.reshape(3, 4)
print(m.T.shape, m.ravel()[:5])
print(np.vstack([m[0], m[1]]))
print(np.concatenate([m, m], axis=1).shape)▸ Ожидаемый результат
[[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11]] (2, 6) (4, 3) [0 1 2 3 4] [[0 1 2 3] [4 5 6 7]] (3, 8)
Линейная алгебра: np.linalg
Произведение матриц — это не поэлементное умножение. A * B перемножает соответствующие элементы, а A @ B (или np.dot) вычисляет математическое произведение матриц: строка A скалярно умножается на столбец B. Для этого число столбцов A должно равняться числу строк B.
- aᵢₖэлемент A в строке i, столбце k
- bₖⱼэлемент B в строке k, столбце j
- nобщее число столбцов A и строк B
import numpy as np
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
print(A * B)
print(A @ B)
print(np.dot(A, B).shape, (np.ones((2, 3)) @ np.ones((3, 5))).shape)▸ Ожидаемый результат
[[ 5 12] [21 32]] [[19 22] [43 50]] (2, 2) (2, 5)
Система линейных уравнений в матричной форме записывается как Ax = b. Если det A ≠ 0, у системы единственное решение и у A есть обратная матрица: x = A⁻¹b. Для матрицы 2 × 2 всё можно посчитать вручную:
- det Aопределитель; если он равен 0, матрица вырождена и обратной нет
- A⁻¹обратная матрица: A · A⁻¹ = I (единичная матрица)
3 тетради и 2 ручки стоят 5,5 маната, 1 тетрадь и 4 ручки — 3,5 маната. Найди цены тетради и ручки матричным методом.
Показать решениеСкрыть решение
det A = 3 · 4 − 2 · 1 = 10 ≠ 0 — решение единственное.
A⁻¹ = (1/10) · [[4, −2], [−1, 3]].
x = A⁻¹b = (1/10) · (4 · 5,5 − 2 · 3,5; −5,5 + 3 · 3,5) = (1/10) · (15; 5) = (1,5; 0,5).
Тетрадь стоит 1,5 маната, ручка — 0,5 маната. Проверка: 3 · 1,5 + 2 · 0,5 = 5,5 ✓.
С тремя неизвестными (тетрадь, ручка, линейка) ручной расчёт становится длинным, а NumPy решает систему одной строкой:
import numpy as np
A = np.array([[2.0, 1.0, 3.0],
[1.0, 3.0, 1.0],
[4.0, 2.0, 0.0]])
b = np.array([9.5, 5.0, 7.0])
x = np.linalg.solve(A, b)
print(x.round(2))
print(np.allclose(A @ x, b))
print(round(np.linalg.det(A), 2))
x_inv = np.linalg.inv(A) @ b
print(np.allclose(x, x_inv))▸ Ожидаемый результат
[1.5 0.5 2. ] True -30.0 True
np.allclose проверяет, что решение удовлетворяет уравнениям; определитель выведен с округлением.Собственный вектор — это ненулевой вектор v, который матрица A только растягивает или сжимает, не меняя его направления; коэффициент растяжения λ — собственное значение. Значения λ находят из характеристического уравнения. Применения: метод главных компонент (PCA), алгоритм PageRank в Google, собственные частоты колебаний мостов и зданий, устойчивость динамических систем.
- vсобственный вектор (v ≠ 0)
- λсобственное значение
- Iединичная матрица
Второе равенство — характеристическое уравнение: система (A − λI)v = 0 имеет ненулевое решение, только если определитель равен 0.
Найди собственные значения и собственные векторы матрицы A = [[2, 1], [1, 2]].
Показать решениеСкрыть решение
λ = 3: (A − 3I)v = 0 → −v₁ + v₂ = 0 → v = (1; 1).
λ = 1: v₁ + v₂ = 0 → v = (1; −1).
Матрица растягивает направление (1; 1) в 3 раза, а направление (1; −1) оставляет без изменений.
import numpy as np
A = np.array([[2.0, 1.0], [1.0, 2.0]])
values, vectors = np.linalg.eigh(A)
print(values.round(4))
v = vectors[:, 1]
print(np.allclose(A @ v, values[1] * v))
print(round(np.linalg.det(A), 4), round(values.prod(), 4))
print(round(np.trace(A), 4), round(values.sum(), 4))▸ Ожидаемый результат
[1. 3.] True 3.0 3.0 4.0 4.0
eigh: он возвращает собственные значения по возрастанию. Векторы нормированы на длину 1, но их знаки могут отличаться в разных библиотеках.Линейная алгебра ведёт прямо к машинному обучению. Подбор наилучшей прямой через точки (метод наименьших квадратов) сводится к нормальному уравнению. Если добавить к матрице X столбец из единиц, свободный член тоже найдётся автоматически:
- Xматрица признаков n × d (первый столбец — единицы)
- yвектор целевых значений
- wкоэффициенты, минимизирующие сумму квадратов ошибок
import numpy as np
rng = np.random.default_rng(7)
x = rng.uniform(0, 10, 50)
y = 2.5 * x + 1.0 + rng.normal(0, 1, 50)
X = np.column_stack([np.ones_like(x), x])
w = np.linalg.solve(X.T @ X, X.T @ y)
print(w.round(3))
w2, *_ = np.linalg.lstsq(X, y, rcond=None)
print(w2.round(3))▸ Ожидаемый результат
[1.236 2.458] [1.236 2.458]
np.linalg.lstsq; LinearRegression в scikit-learn решает ту же задачу.Приведи каждый столбец матрицы X к диапазону [0; 1] без цикла, с помощью broadcasting (min–max нормализация): (X − минимум столбца) / (максимум столбца − минимум столбца). Выведи результат, округлённый до 2 знаков после запятой.
import numpy as np
X = np.array([[2.0, 50.0, 1.0],
[4.0, 80.0, 3.0],
[6.0, 65.0, 2.0],
[8.0, 95.0, 5.0]])
# min-max scale every column▸ Ожидаемый результат
[[0. 0. 0. ] [0.33 0.67 0.5 ] [0.67 0.33 0.25] [1. 1. 1. ]]
На рынке: 2 кг яблок + 1 кг груш + 1 кг винограда = 9,5 маната; 1 кг яблок + 2 кг груш + 1 кг винограда = 10 манатов; 1 кг яблок + 1 кг груш + 3 кг винограда = 13,5 маната. С помощью np.linalg.solve найди цену 1 кг яблок, груш и винограда и выведи её с округлением до 2 знаков, затем выведи определитель матрицы с округлением до 1 знака.
import numpy as np
# A: kilograms in each purchase, b: amounts paid (manat)
A = np.array([[2.0, 1.0, 1.0],
[1.0, 2.0, 1.0],
[1.0, 1.0, 3.0]])
b = np.array([9.5, 10.0, 13.5])
# solve A x = b and print the prices, then the determinant▸ Ожидаемый результат
[2. 2.5 3. ] 7.0
Главное
- Векторизация заменяет циклы Python на ufunc, работающие на C, и ускоряет код в десятки и сотни раз.
- Broadcasting сравнивает формы справа: длины должны совпадать, либо одна из них должна быть 1.
v[:, None]иkeepdims=Trueсогласуют формы, превращая вектор в столбец.*— поэлементное умножение,@— матричное произведение: (m × n) @ (n × p) → (m × p).- Решай Ax = b через
np.linalg.solve; при det A = 0 единственного решения нет. - Av = λv; сумма собственных значений равна следу, произведение — определителю.
Проверь себя
Вопросов: 10. Каждый правильный ответ приносит XP.