Научные вычисления
NumPy и линейная алгебра
Цели урока
- Понимать ndarray как однородный массив в непрерывной памяти и почему это дешевле списков
- Применять broadcasting правильно и предсказывать форму результата
- Использовать np.linalg для решения систем, SVD, собственных значений; знать когда НЕ брать inv()
- Заменять Python-циклы на vectorized операции и понимать ускорение через BLAS/LAPACK
Предварительные знания
Климатическая модель ECMWF обрабатывает 100 TB данных каждый день для прогноза погоды. NumPy - фундамент PyTorch, TensorFlow, pandas. scipy.linalg - Fortran LAPACK под капотом. CERN Large Hadron Collider: 15 петабайт данных в год через ROOT framework. Матрица 10 000 x 10 000? Одна строка NumPy - секунды. Три вложенных цикла Python - часы.
- **PyTorch и TensorFlow** - Tensor = ndarray + autograd; весь deep learning строится на numpy broadcasting и matmul semantics
- **ECMWF прогноз погоды** - 100 TB ежедневно, NumPy/SciPy для постобработки; точность прогноза 10 дней улучшена в 2x за 20 лет
- **CERN ROOT framework** - 15 ПБ/год; Python-биндинги через NumPy arrays для particle physics analysis
- **NumPy SVD в рекомендациях** - Netflix Prize 2009 выиграли через матричную факторизацию; np.linalg.svd за 3 строки
- **Scipy LAPACK** - тот же Fortran код, что в aerospace engineering и CFD симуляциях 1980-х, теперь в pip install scipy
Тревис Олифант и рождение NumPy
В 2005 году Тревис Олифант объединил два конкурирующих проекта - Numeric (1995, Джим Хугунин) и Numarray (2001, STScI) - в единый пакет NumPy. Под капотом NumPy вызывает BLAS (Basic Linear Algebra Subprograms, 1979) и LAPACK (1992, Андерсон и др.) - библиотеки на Fortran которым уже больше сорока лет. Именно поэтому np.dot работает в сотни раз быстрее Python-цикла: SIMD-инструкции и кэш-оптимизированный GEMM из 80-х. В 2020 публикация "Array programming with NumPy" в Nature официально признала NumPy фундаментом современной науки о данных.
ndarray - основа NumPy
Джон фон Нейман и Герман Голдстайн, 1947 - первые численные вычисления
В 1947 году Джон фон Нейман и Герман Голдстайн запустили обращение матрицы на ENIAC - первые серьёзные численные вычисления на современном компьютере. ENIAC не имел float64, не знал о broadcasting и векторизации. NumPy 2005 года несёт тот же принцип: однородный массив чисел в памяти. Но на порядки быстрее. PyTorch, TensorFlow, pandas - все они ndarray под капотом.
Климатическая модель ECMWF обрабатывает 100 TB данных каждый день для прогноза погоды. NumPy - фундамент PyTorch, TensorFlow, pandas. CERN Large Hadron Collider: 15 петабайт данных в год через ROOT framework, и большая часть анализа - в NumPy/SciPy. Python-списки для этого не годятся. **ndarray** - N-мерный массив фиксированного типа в непрерывном блоке памяти. Именно это обеспечивает скорость в сотни раз выше списков.
**shape** - кортеж размерностей (строки, столбцы, ...). **dtype** - тип элементов (float64, int32, complex128). Все элементы одного типа - это позволяет хранить данные компактно и обращаться без накладных расходов Python. PyTorch Tensor и NumPy ndarray разделяют одну память при `torch.from_numpy()` - нулевое копирование.
Slicing в NumPy возвращает **view** - не копию, а окно в тот же блок памяти. Изменение view меняет оригинал. Для независимой копии: `M[0, :].copy()`. Это не баг - это архитектурное решение для нулевого копирования в ML pipeline.
| Создание | Функция | Пример |
|---|---|---|
| Нули | np.zeros(shape) | np.zeros((3, 3)) |
| Единицы | np.ones(shape) | np.ones((2, 5)) |
| Единичная матрица | np.eye(n) | np.eye(4) |
| Случайные | np.random.randn(shape) | np.random.randn(3, 3) |
| Диапазон | np.linspace(a, b, n) | np.linspace(0, 1, 50) |
Что произойдёт при выполнении: row = M[0, :]; row[0] = 42?
Broadcasting
Нормализация 1000 сэмплов по 5 фичам. Вычесть среднее и поделить на стандартное отклонение. Наивный подход: цикл по 1000 строкам. NumPy подход: одна строка. Матрица (1000, 5) минус вектор (5,) - это **broadcasting**: NumPy автоматически растягивает вектор до матрицы без создания копии в памяти. Нулевой overhead.
**Правила broadcasting:** 1. Размерности сравниваются справа налево. 2. Размерности совместимы, если равны ИЛИ одна из них = 1. 3. Отсутствующие размерности считаются равными 1. Если правила не выполняются - ошибка.
Broadcasting - не просто удобство. NumPy не создаёт увеличенную копию в памяти - операция выполняется «на лету». PyTorch использует тот же механизм: градиенты, loss per sample, batch normalization statistics - всё через broadcasting. Понять broadcasting = понять половину PyTorch.
Какой результат даст операция np.ones((3, 4)) + np.ones((3, 1))?
numpy.linalg
scipy.linalg - Fortran LAPACK под капотом. Тот самый LAPACK, написанный в 1970-80-х, оптимизированный под каждую архитектуру процессора. NumPy вызывает ту же библиотеку. Решение систем уравнений, собственные значения, SVD - всё это не Python, это decades of numerical analysis в скомпилированном коде.
**np.linalg.solve(A, b)** - O(n^3), но намного быстрее и точнее np.linalg.inv(A) @ b. Никогда не вычисляйте обратную матрицу для решения системы - это численно нестабильно. solve() использует LU-разложение с partial pivoting: стандарт в научных вычислениях с 1960-х.
| Функция | Что делает | Сложность |
|---|---|---|
| np.linalg.solve(A, b) | Решает Ax=b | O(n^3) |
| np.linalg.eig(A) | Собственные значения/векторы | O(n^3) |
| np.linalg.svd(A) | Сингулярное разложение | O(min(m,n)*m*n) |
| np.linalg.inv(A) | Обратная матрица (избегайте!) | O(n^3) |
| np.linalg.det(A) | Определитель | O(n^3) |
| np.linalg.norm(x) | Норма вектора/матрицы | O(n) |
**SVD** (Singular Value Decomposition) - один из самых полезных инструментов в ML. PCA = SVD на центрированных данных. Рекомендательные системы Netflix (collaborative filtering) - матричная факторизация через SVD. Сжатие изображений (JPEG-подобное) - усечённый SVD. Псевдообратная матрица Moore-Penrose - через SVD. Одна функция, четыре killer applications.
Почему np.linalg.solve(A, b) предпочтительнее np.linalg.inv(A) @ b?
Vectorization
PyTorch обучает ResNet-50 на ImageNet за 76 часов на одном V100 GPU. То же обучение на чистом Python с циклами заняло бы годы. Разница - не GPU. Разница - **vectorization**: замена Python-циклов на операции над целыми массивами. NumPy выполняет их в скомпилированном C-коде с SIMD-инструкциями (SSE4, AVX-512). Ускорение: 50-300x типично, 500x для матричного умножения.
**Почему vectorization так быстра:** 1. Нет накладных расходов интерпретатора Python на каждую итерацию. 2. Данные в ndarray хранятся непрерывно - процессор загружает их в L1/L2 кэш за одну операцию. 3. SIMD-инструкции AVX-512 обрабатывают 8 float64 за такт процессора.
| Паттерн цикла | Vectorized эквивалент | Ускорение |
|---|---|---|
| sum(a[i]*b[i]) | np.dot(a, b) | ~300x |
| for + if > threshold | a[a > threshold] | ~100x |
| тройной цикл C=A*B | A @ B | ~500x |
| for + math.sin(x) | np.sin(arr) | ~100x |
Правило: если есть `for` в NumPy-коде - почти наверняка есть vectorized альтернатива. Исключение - рекуррентные вычисления, где шаг N зависит от шага N-1. Там нужен Numba (@jit) или Cython. Но таких случаев в ML - меньшинство.
Под капотом NumPy вызывает **BLAS** и **LAPACK** - библиотеки, оптимизированные десятилетиями под Intel MKL, OpenBLAS, ATLAS. Когда пишем `np.dot(a, b)`, Python лишь передаёт указатели в C-функцию. Вся работа - в Fortran/C assembly.
Python слишком медленный для научных расчётов
Python-цикл действительно медленный, но NumPy/SciPy вызывают оптимизированные C/Fortran библиотеки (BLAS, LAPACK). Vectorized код на Python работает со скоростью, сопоставимой с C.
NumPy - тонкая Python-обёртка над скомпилированным кодом. np.dot(a, b) передаёт указатели в C-функцию, которая выполняет всю работу с SIMD. Python занимает <1% времени выполнения.
Какая операция НЕ поддаётся прямой vectorization в NumPy?
Ключевые идеи
- **ndarray** - однородный массив фиксированного типа в непрерывной памяти; основа PyTorch, TensorFlow, pandas
- **Broadcasting**: сравнение справа налево, совместимость если равны или одна = 1; нулевое копирование в памяти
- **np.linalg.solve** через LU-разложение - численно стабильнее inv(); SVD = PCA, рекомендации, псевдообратная
- **Vectorization** даёт 100-500x vs Python-циклов: SIMD AVX-512, BLAS/LAPACK под капотом
- **View vs copy**: slice возвращает view - изменение меняет оригинал; .copy() для независимости
- **Исключение для vectorization**: рекуррентные зависимости x[i] = f(x[i-1]) требуют Numba или Cython
Связанные темы
NumPy - фундамент всего научного стека Python:
- Системы линейных уравнений — np.linalg.solve - начало; дальше LU, Cholesky, итеративные методы
- Введение в научные вычисления — NumPy - основной инструмент реализации моделей и симуляций
- LU-разложение — Математика под капотом np.linalg.solve
Вопросы для размышления
- Почему NumPy требует, чтобы все элементы массива были одного типа? Какие преимущества это даёт?
- В каких случаях broadcasting может привести к неожиданным результатам? Придумайте пример.
- Если vectorization так быстра, почему не все Python-программы используют NumPy?
Связанные уроки
- sci-01 — Введение в научные вычисления - контекст для NumPy
- sci-03 — np.linalg.solve - старт; дальше LU, Cholesky, итеративные методы
- la-21-lu — LU-разложение под капотом np.linalg.solve
- la-15-svd — SVD в NumPy - частный случай полного теоретического SVD
- ml-03-math-foundations — NumPy - фундамент матем. базиса всего ML стека
- calc-17-multivariable — Broadcasting обобщает операции над многомерными массивами как многомерный анализ
- la-04-matrix-ops