Рандомизированная линейная алгебра

Матрица данных 10 миллионов пользователей × 100 тысяч продуктов. Полное SVD — невозможно: O(mn·min(m,n)) ≈ 10¹⁷ операций. Рандомизированные алгоритмы снижают сложность до O(mnk) при теоретических гарантиях близости к точному решению. Это не «приближение из-за лени» — это математически строгая теория.

Рандомизированное SVD

Задача: вычислить лучшее ранг-k приближение матрицы A ∈ ℝ^{m×n}: min_{rank(B)≤k} ||A − B||_F.

Точный алгоритм: SVD A = UΣVᵀ, Aₖ = U_{1:k}Σ_{1:k}Vᵀ_{1:k}. Сложность O(mn²).

Randomized SVD (Halko, Martinsson, Tropp, 2011):

Шаг 1. Случайная матрица Ω ∈ ℝ^{n×k}: Ωᵢⱼ ~ N(0,1/k). Sketch Y = AΩ ∈ ℝ^{m×k}. Смысл: Y аппроксимирует «пространство столбцов» A — проецируем A на k случайных направлений.

Шаг 2. QR-разложение Y = QR, Q ∈ ℝ^{m×k} — ортонормированный базис пространства столбцов Y.

Шаг 3. B = QᵀA ∈ ℝ^{k×n}. Матрица B намного меньше: k × n.

Шаг 4. SVD(B) = ŪΣVᵀ, U = QŪ.

Итог: U ∈ ℝ^{m×k}, Σ ∈ ℝ^{k×k}, V ∈ ℝ^{n×k} — аппроксимация SVD. Сложность: O(mnk) vs O(mn²).

Теоретические гарантии: с вероятностью ≥ 1−δ:

||A − QQᵀA||F ≤ (1 + O(√(k/δ))) σ{k+1}(A)

Здесь σ_{k+1}(A) — (k+1)-е сингулярное значение, нижняя граница для любого ранг-k приближения. Практически: k_actual = k+p (p=10 «oversampling» итераций) даёт почти точное решение.

Degree-q power iteration: Вместо Y = AΩ используем Y = (AAᵀ)^q AΩ. При q=1–2: ошибка ~σ_{k+1}^{2q+1} — экспоненциально улучшается при q. Стоимость: q дополнительных умножений на A.

Sklearn random SVD по умолчанию: n_iter=4 — 4 power iterations.

Приближённый поиск ближайших соседей (ANNS)

Задача: для запроса q найти x* = argmin_{x∈X} d(q, x). Exact kNN: O(nd) за запрос. При n=10^8, d=768 (BERT embeddings) — невозможно.

LSH (Locality Sensitive Hashing): семейство хэш-функций h такое, что P(h(x)=h(y)) велика при d(x,y) мала и мала при d(x,y) велика.

Для L2 (Дог, 2004): h_{a,b}(x) = ⌊(aᵀx + b)/w⌋, где a ~ N(0,I), b ~ Uniform[0,w]. Близкие точки попадают в один «ведёрко» с высокой вероятностью.

Амплификация: L независимых таблиц × k хэш-функций в каждой. Запрос: смотрим в L таблицах, берём union of buckets. Время: O(d·n^ρ), ρ = 1/(2L²·w² + 1) < 1.

HNSW (Malkov & Yashunin, 2018): иерархический граф навигации. Каждый уровень = граф ближайших соседей с разной плотностью. Поиск: начинаем с верхнего (разреженного) уровня → быстро приближаемся к зоне запроса → спускаемся вниз. Время: O(log n). SOTA для большинства задач ANNS.

FAISS (Facebook AI Similarity Search): библиотека с GPU-ускорением. Product Quantization: сжимаем векторы в коды (4–8 байт вместо 768×4=3072 байт). IVF (Inverted File): кластеризуем пространство, ищем только в ближайших кластерах. Поиск 10^8 векторов за <1 мс на GPU.

Потоковые алгоритмы: большие данные малой памятью

Задача: входит поток из N элементов x₁, x₂,..., xₙ. Память O(N) недопустима. Нужно вычислить статистику потока за один проход.

Count-Min Sketch (Cormode & Muthukrishnan, 2004): оцениваем частоту f(x) при малой памяти.

Структура: матрица счётчиков C[d][w] (d рядов × w столбцов). d независимых хэш-функций h₁,...,h_d: U → {1,...,w}.

Обновление для элемента x: для i=1,...,d: C[i][hᵢ(x)] += 1.

Оценка f̂(x) = min_i C[i][hᵢ(x)] (минимум по рядам).

Гарантии: P(|f̂(x) − f(x)| ≤ εN) ≥ 1−δ при w = ⌈e/ε⌉, d = ⌈ln(1/δ)⌉. Память: O(1/(ε) · log(1/δ)) — не зависит от N!

HyperLogLog (Flajolet et al., 2007): оцениваем число уникальных элементов (cardinality). Идея: хэшируем x → равномерный [0,1]. Максимальный показатель последовательных нулей в начале двоичного представления → оценка log₂(n). Средний по m=2^b подрегистрам. Ошибка: ~1.04/√m. Память: O(log log N) бит. В Redis реализован, используется для Count Distinct в аналитике.

Morris Counter: считать n до ~2^{64} в O(log log n) бит. Храним k = ⌊log₂(n)⌋ + случайная ошибка. Инкремент: k ← k+1 с вероятностью 2^{-k}. Оценка: 2^k. Дисперсия = O(n²).

Численный пример

Рандомизированное SVD Netflix Prize (480K пользователей × 17K фильмов):

Истинный ранг матрицы ≈ 100
Точное SVD: 480K × 17K × 100 ≈ 816 × 10⁹ операций (~4 часа)
Random SVD (k=100, q=2): ~240 × 10⁹ → ~1 час, ошибка <1%

HNSW поиск по 10M BERT-эмбеддингов (768d): индексация 10 минут, поиск top-10 за 0.3 мс (vs 4.5 секунды для exact brute-force).

Задание: (1) Реализуйте рандомизированное SVD без sklearn для матрицы рейтингов MovieLens (100K). Сравните точность (Frobenius error) с truncated SVD для k=10, 50, 100 при q=0, 1, 2 power iterations. (2) Реализуйте Count-Min Sketch и HyperLogLog в Python. Протестируйте на потоке из 10M слов английского текста. Сравните оцениваемое число уникальных слов с истинным.