PCA на больших матрицах: практические решения, когда Sklearn падает на 40k измерений

Проблема: Sklearn PCA упирается в стену на больших масштабах

Матрица 40 000 × 40 000 значений float32 занимает примерно 6 ГБ оперативной памяти. Это звучит вполне подъёмно для машины со 128 ГБ. Тем не менее стандартный SVD-солвер sklearn падает — причём не с ошибкой нехватки памяти, а зачастую с загадочными сообщениями вроде free(): invalid size. Матрица помещается в память. Вычисление — нет.

Это распространённый камень преткновения в representation learning, где матрицы эмбеддингов регулярно достигают десятков тысяч измерений. Стандартное полное SVD-разложение, которое sklearn использует внутри, выделяет промежуточные матрицы, многократно увеличивающие потребление памяти. Для квадратной матрицы 40k фактическое пиковое потребление памяти может превышать 50–80 ГБ, когда подпрограммы LAPACK создают свои внутренние копии и рабочие массивы.

Проще говоря: сами данные помещаются в оперативную память, но алгоритму для работы нужно гораздо больше места, чем занимают сами данные.

Почему стандартный Sklearn PCA не справляется с большими квадратными матрицами

Sklearn PCA с параметром по умолчанию svd_solver='full' вычисляет полное сингулярное разложение через LAPACK. Вычислительная сложность составляет O(n² × m), где n и m — размерности матрицы, а потребление памяти масштабируется аналогично. Для матрицы 40k × 40k солверу необходимо выделить несколько плотных матриц тех же размерностей для промежуточных вычислений.

Есть и более тонкая проблема. Подпрограммы LAPACK в некоторых сборках имеют внутренние ограничения на размер или паттерны выделения памяти, которые приводят к сбоям на больших матрицах вне зависимости от доступного объёма RAM. Крах free(): invalid size, который возникает на матрицах около 30k × 28k — тогда как матрицы поменьше, вроде 20k × 28k, работают нормально — часто указывает на целочисленное переполнение или проблемы выравнивания памяти в базовых Fortran-подпрограммах, а не на реальную нехватку RAM.

Честно говоря: даже если бы полное SVD работало без ошибок, на матрице 40k × 40k оно было бы мучительно медленным. Даже на машине, где оно не падает, вычисление может занять часы. Для большинства задач representation learning вам не нужны все 40 000 компонент — а это открывает дорогу к гораздо более быстрым подходам.

Решение 1: Рандомизированное SVD — самый быстрый способ

Если вам нужны только верхние k компонент (а в representation learning это почти всегда так), рандомизированное SVD полностью пропускает вычисление полного разложения. Вместо этого оно аппроксимирует верхние сингулярные векторы с помощью случайных проекций.

Согласно бенчмаркам Амеде д'Абовиля, реализация рандомизированного SVD в sklearn значительно улучшилась после внедрения методов численной стабилизации из работы Halko и др. и исправления ключевых ошибок производительности. Sklearn теперь автоматически переключается на рандомизированный PCA для достаточно больших матриц при использовании солвера по умолчанию, но вы можете задать это явно:

from sklearn.decomposition import PCA

pca = PCA(n_components=256, svd_solver='randomized', iterated_power=2)
result = pca.fit_transform(matrix)

Реальные цифры: потребление памяти для рандомизированного SVD снижается примерно до 2 × n × n_components вместо n × m для полного метода, как отмечено в документации sklearn. Для матрицы 40k × 40k с 256 целевыми компонентами это означает примерно 80 МБ вместо 12+ ГБ для промежуточного хранения. Временная сложность снижается с O(n² × m) до O(n² × n_components).

Для матрицы 40k × 40k с 256 компонентами рандомизированное SVD обычно завершается за минуты, а не за часы — и использует лишь малую долю памяти.

Решение 2: IncrementalPCA — когда даже рандомизированное SVD слишком много

Если матрица слишком велика, чтобы поместиться в память целиком — или если вы работаете с потоковыми данными — IncrementalPCA из sklearn обрабатывает данные батчами. Его потребление памяти постоянно: O(batch_size × n_features). Он никогда не загружает весь датасет целиком.

from sklearn.decomposition import IncrementalPCA
import numpy as np

ipca = IncrementalPCA(n_components=256, batch_size=1000)

# Может работать с memory-mapped файлами — никогда не загружает полную матрицу
data = np.memmap('embeddings.dat', dtype='float32', mode='r', shape=(40000, 40000))

for i in range(0, 40000, 1000):
    ipca.partial_fit(data[i:i+1000])

result = ipca.transform(data)

Ключевое преимущество: IncrementalPCA работает с файлами np.memmap, то есть данные хранятся на диске, а в оперативную память попадают лишь небольшие фрагменты. Вычислительная нагрузка на каждый батч составляет O(batch_size × n_features²), но одновременно в памяти находятся только 2 × batch_size семплов.

Компромисс — скорость. IncrementalPCA выполняет n_samples / batch_size отдельных SVD-разложений вместо одного. Для матрицы 40k × 40k с batch_size=1000 это 40 последовательных SVD. С другой стороны, каждое из них маленькое и быстрое.

Главный вывод для бизнеса: IncrementalPCA обменивает время вычислений на гарантированные границы потребления памяти. Если продакшн-пайплайн должен работать на машине с фиксированным объёмом памяти без сюрпризов — это самый надёжный выбор.

Решение 3: Dask для распределённого и out-of-core PCA

Для команд, уже использующих Dask, dask-ml предоставляет собственный IncrementalPCA, который параллелизирует батчевую обработку по ядрам или машинам. API повторяет sklearn, поэтому переход требует минимальных изменений в коде.

Dask наиболее ценен, когда матрица целиком не помещается в память одной машины — сотни тысяч строк, или когда шаг PCA является частью более крупного распределённого пайплайна. Для матрицы 40k × 40k на одной машине со 128 ГБ он добавляет сложность без существенного выигрыша по сравнению с IncrementalPCA из sklearn. Вот наша рекомендация: используйте Dask только если вы уже работаете в экосистеме Dask или ожидаете значительного роста ваших матриц.

Решение 4: fbpca и cuML — специализированные библиотеки

fbpca (рандомизированный PCA от Facebook) исторически был быстрее и численно стабильнее, чем рандомизированное SVD в sklearn. Как отмечает д'Абовиль, sklearn с тех пор сократил этот разрыв, но fbpca остаётся надёжным вариантом — особенно если вы уже его используете:

import fbpca
U, s, Vt = fbpca.pca(matrix, k=256, n_iter=2)

cuML (RAPIDS) переносит всё вычисление на GPU. Матрица 40k × 40k float32 требует 6 ГБ VRAM, что помещается на большинстве современных GPU. Ускорение впечатляющее — то, что занимает минуты на CPU, может завершиться за секунды на GPU. Ограничение очевидно: вам нужен GPU с достаточным объёмом VRAM.

Выбор правильного подхода

По нашему опыту работы с продакшн ML-пайплайнами, дерево решений простое:

1. Сначала попробуйте рандомизированное SVD. Оно решает 90% проблем с большим PCA одним изменением параметра.
2. Если матрица всё ещё не помещается в память, используйте IncrementalPCA с memory-mapped файлами.
3. Если критична скорость, переходите на GPU с cuML.
4. Если матрица будет продолжать расти, инвестируйте в инфраструктуру Dask.

Практические советы для PCA на больших матрицах

Снизьте точность перед вычислением. Если ваши эмбеддинги в float64, приведение к float32 вдвое уменьшает потребление памяти. Большинство моделей representation learning и так выдают float32.

Сначала подвыборка, валидация потом. Запустите PCA на случайном подмножестве из 10k, чтобы определить, сколько компонент захватывают значимую дисперсию. Если 128 компонент объясняют 95% дисперсии на подвыборке, они, скорее всего, объяснят аналогичную долю на полной матрице. Это позволяет избежать траты вычислительных ресурсов на исследовательские запуски.

Проверьте вашу сборку LAPACK. Крах free(): invalid size часто специфичен для определённых сборок LAPACK/BLAS. Переключение с OpenBLAS на MKL (или наоборот) может полностью устранить падение, даже с полным солвером. Это не настоящее решение проблемы производительности, но оно устраняет загадочный крах.

Используйте copy=False для экономии памяти. По умолчанию sklearn копирует вашу матрицу перед обработкой. Установка copy=False в конструкторе PCA модифицирует входные данные на месте, экономя объём RAM, равный одной полной копии.

Часто задаваемые вопросы

Как выполнить PCA на очень больших матрицах, когда стандартный SVD sklearn падает, несмотря на достаточный объём RAM?

Переключитесь на svd_solver='randomized' для самого быстрого решения или используйте IncrementalPCA с батчевой обработкой. Оба подхода позволяют избежать полного SVD-вычисления, из-за которого LAPACK превышает лимиты памяти, даже когда сами данные помещаются в RAM.

Почему PCA в sklearn падает с ошибкой «free(): invalid size» на матрицах около 30k × 28k, хотя матрицы поменьше работают нормально?

Это обычно указывает на проблему в базовой библиотеке LAPACK/BLAS — либо целочисленное переполнение при выделении рабочего пространства, либо проблему выравнивания памяти при определённых размерах матриц. Смена реализации BLAS (например, с OpenBLAS на MKL) или использование svd_solver='randomized' полностью обходит проблемный путь выполнения кода.

Когда следует использовать IncrementalPCA вместо рандомизированного SVD?

Используйте рандомизированное SVD, когда матрица помещается в оперативную память — это быстрее и проще. Используйте IncrementalPCA, когда матрица не помещается в память, когда вам нужно гарантированное постоянное потребление памяти или когда вы обрабатываете потоковые данные. IncrementalPCA с np.memmap может обрабатывать матрицы любого размера на машинах с ограниченным объёмом RAM.

Могут ли рандомизированные методы разложения существенно сократить как время вычислений, так и потребление памяти для PCA на данных размерностью 40k+?

Да. Согласно деталям реализации sklearn, рандомизированное SVD снижает потребление памяти с O(n × m) до O(n × n_components) и временную сложность с O(n² × m) до O(n² × n_components). Для матрицы 40k × 40k, редуцированной до 256 компонент, это примерно 150-кратное сокращение промежуточной памяти и пропорциональное ускорение.