The Ground Truth

Первое направление - сравнить эффективность "классических" SVD алгоритмов против "продвинутых". Под классическими понимаются стандартные алгоритмы, используемые в LAPACK для полного SVD: dgesvd и dgesdd. dgesvd использует implicit QR-iteration with zero shift (LAWN #3, раздел 3). dgesdd использует рекурсивный "divide-and-conquer" алгоритм (статья 1 и статья 2). Это относительно старые алгоритмы, но все современные пакеты/среды типа numpy, R и всего такого для операций линейной алгебры используют вызовы к LAPACK и все сводится к этим двум процедурам.

Есть альтернативный алгоритм, который называется dqds (differential quotient-difference algorithm with shifts). Судя по литературе, этот алгоритм очень хорошо находит собственные значения, но только их. То есть не находит собственные векторы. Про основы dqds есть немного здесь (стр 26).

Задача: сравнить, насколько dqds лучше "классики". И классика, и dqds реализованы в LAPACK, причем если вызвать dgesvd с указанием не считать сингулярные векторы, то будет использоваться именно dqds (согласно документации).

Для оценки точности нахождения сингулярных значений, очевидно, нужно знать истинные сингулярные значения. Для этого нужно построить матрицу с заданным спектром (набором сингулярных значений). Для этого можно применить SVD "наоборот": имея заданный набор сингулярных значений $\sigma_1 \ldots \sigma_k$, делаем из них диагональную матрицу $\Sigma = \mathrm{diag}(\sigma_1 \ldots \sigma_k)$. Тогда матрицу A с нужным спектром можно построить как

\underset{m \times n}{A} = \underset{m \times m}{U}\  \underset{m \times n}\Sigma\ \underset{n \times n}V^{T},

где $\Sigma$ - диагональная, $U$ и $V$ - ортогональные. (если $k < m$ или $n$, то в $\Sigma$ просто будет строка/столбец из нулей). Получить случайные ортогональные матрицы $U$ и $V^T$ нужного размера можно таким методом:

1. Сгенерировать две случайных квадратных матрицы квадратную матрицу нужного размера $\underset{m \times m}{M_1}$ и $\underset{n \times n}{M_2}$ (стоит только следить за числом обусловленности, чтобы не было слишком большим - проблем в алгоритмах вылезает много, и мы хотим быть уверены хотя бы в наших исходных данных).
2. Выполнить их QR-разложение: $M_1 = Q_1R_1$, $M_2 = Q_2R_2$
3. Из QR разложение берем только $Q$: $\underset{m \times m}{U} = \underset{m \times m}{Q_1}, \quad \underset{n \times n}{V^T} = \underset{n \times n}{Q_2}$

Итого, план действий:

1. Создать матрицу с заданным спектром
2. Вычислить ее SVD стандартными алгоритмами
3. Вычислить сингулярные значения с помощью dqds
4. Сравнить, насколько лучше (или вдруг хуже) справляется dqds.

Нахождение сингулярных векторов по сингулярным значениям

Результаты исследования dqds дадут нам ответ на вопрос: можно ли достаточно точно найти сингулярные значения? Это важная и полезная информация, но ее результат не критичен: если dqds очень точно находит сингулярные значения, мы говорим, что у нас есть точные сингулярные значения, найденные с помощью dqds. Если dqds справляется не очень, мы говорим: предположим, у нас есть точные сингулярные значения, найденные как-нибудь (dqds не единственная наша надежда). Возникает задача: как можно найти сингулярные векторы, зная достаточно точные сингулярные значения? Здесь есть три идеи:

Идея 1: Уточнение поворотов Гивенса в implicit zero-shift QR

Алгоритм в dgesvd основан на итеративном сведении бидиагональной матрицы $B$ к диагональной $\Sigma$ с помощью поворотов Гивенса (процесс подробно описан в LAWN#3):

$$\underset{n \times n}{B_{i+1}} =J_{2n-2} \cdots J_6 J_4 J_2 B_i J_1 J_3 J_5 \cdots J_{2n-3}$$

(всего по $n-1$ поворотов с каждой стороны для $n \times n$ бидиагональной $B$) После некоторого количества итераций $B_i$ станет близка к диагональной с сингулярными значениями ($\hat{\Sigma}$ на главной диагонали, и получается следующее представление исходной бидиагональной $B$:

$$B = (\prod J_L) \hat{\Sigma} (\prod J_R),$$

где итеративно "накопленные" левые и правые произведения поворотов Гивенса и дают ортогональные $U$ и $V$ для сингулярного разложения. При плохих числах обусловленности матрицы, сингулярные векторы и значения найдены не очень точно. Теперь мы используем предположение, что нам достаточно точно, лучше чем из этого алгоритма, известны сингулярные значения. Тогда в выражение выше в центр вместо $\hat{\Sigma}$ мы ставим "истинную" $\Sigma$:

(\prod J_L) \Sigma (\prod J_R) = \tilde{B}

После такого действия $\tilde{B}$ не будет равно исходной $B$, но в этом и суть: мы пытаемся "подкрутить" каждый из поворотов Гивенса $J_k$ c целью минимизировать какую-либо норму (самое очевидное - Фробениусову, но есть другие опции) разницы матриц $\lVert \tilde{B} - B \rVert$. Уточнение одного $J_k$ вещь довольно простая: каждый поворот Гивенса параметризуется всего одним значением угла $\theta$, который мы и изменяем при "подкрутке". Итого у нас есть куча поворотов слева и справа, каждый из которых параметризуется одним углом:

\tilde{B} = J(\theta_1) J(\theta_2) \cdots J(\theta_N) \Sigma J(\theta_{N+1}) J(\theta_{N+2}) \cdots J(\theta_{2N})

и мы решаем оптимизационную задачу:

\min_{\theta_1 \ldots \theta_{2N}} \lVert \tilde{B} - B \rVert

Это можно попробовать сделать, например, с помощью градиентного спуска. Целевая функция очень громоздкая, но концептуально довольно простая (умножение кучи матриц друг на друга). Если хочется заморочиться - можно попробовать реализовать автоматическое дифференцирование этой конструкции в JAX. Итого, если способ рабочий, произведение всех уточненных поворотов Гивенса слева и справа даст нам достаточно приличные матрицы $U$ и $V$, предполагая, что мы знаем достаточно точные сингулярные значения. Однако, чтобы проверить эту идею, нужна реализация этого алгоритма, которую потом можно будет модифицировать нужным образом. Влезать в оптимизированный фортрановский код LAPACK'а представляется гиблым делом, поэтому базовую версию алгоритма нужно будет написать (или найти приличную и понятную реализацию). Здесь потенциальными проблемами мне видится громоздкость всей этой конструкции и большое количество матриц. Поэтому вначале лучше попробовать отработать эту идею на небольших матрицах, но с достаточно плохими числами обусловленности (чтобы был смысл вообще в подкрутке). Если результаты получатся перспективными - будем думать, как это обобщить на большие матрицы.

Идея 2: "обратный" поворот Якоби

Известные довольно давно, но менее распространенные (из-за, как правило, худшей производительности, хотя ситуация меняется) методы вычисления SVD -- односторонний и двусторонний методы Якоби (про них компактно рассказано здесь, раздел 12). Основная идея похожа на implicit zero-shift QR, но применяется к исходной матрице $A$ (не бидиагонализированной):

A_{i+1} = J^T_i A_i K_i, \quad A_k \to \Sigma\ \text{при}\ k \to \infty

Опять же, после достаточно количества таких итераций $A$ сходится к диагональной $\Sigma$, а произведение левых и правых поворотов Якоби (они аналогичны поворотам Гивенса, просто иногда используется другое название) дают ортогональные $U$ и $V$. Повороты $J^T_i$ и $K_i$ выбираются таким образом, чтобы на одной итерации сделать один $2 \times 2$ блок диагональным; как это работает, обсуждается в Golub, Van Loan Matrix Computations, раздел 8.5. Наша задача заключается в том, чтобы, зная сингулярные значения, получить сингулярные векторы. Серия ортогональных поворотов слева и справа сводит исходную $A$ к диагональной $\Sigma$; мы точно знаем исходную $A$; мы предполагаем, что знаем исходную $\Sigma$. Из метода Якоби появляется следующая идея: начать с $\Sigma$ и серией ортогональных преобразований свести ее к исходной $A$: тогда произведение этих ортогональных преобразований и даст нам искомые сингулярные векторы. И в этом основная теоретическая задача этой идеи: придумать способ строить ортогональные повороты, которые итеративно превращают диагональную $\Sigma$ в исходную $A$. Исходный метод Якоби итеративно уменьшает норму внедиагональных элементов матрицы; возможно, обратные повороты стоит выбирать так, чтобы каждый поворот приближал текущую матрицу к исходной $A$ по норме разницы некоторых элементов? Здесь есть некоторые идеи, как только они обретут более-менее приличную форму, я их добавлю.

Идея 3: "наивный" метод

"Классическая" связь между сингулярным и собственным разложением заключается в следующем: левые и правые сингулярные векторы $A$ -- это правые сингулярные векторы $AA^T$ и $A^TA$ соответственно, сингулярные значения $A$ -- квадратные корни собственных значений $A^TA$. Однако подобное сведение задачи сингулярного разложения к задаче нахождения обычных собственных векторов и значений плохо подходит из численных соображений: умножение матрицы саму на себя возводит в квадрат число обусловленности и приводит к значительной потери точности практически на ровном месте. Однако есть альтернативная формулировка задачи сведения проблемы сингулярного разложения к проблеме собственного разложения, избегающая умножения матрицы саму на себя. Она заключается в составлении матрицы вида

\begin{bmatrix}
0 & A^T\\
A & 0
\end{bmatrix}

для исходной $A$ или

\begin{bmatrix}
0 & B\\
B^T & 0
\end{bmatrix}

для бидиагональной $B$. Соотношение между сингулярными числами/векторами и собственными числами/векторами для первого случая описаны в Matrix Computations, раздел 8.6.1, для второго случая - здесь на странице 121. Я думаю, можно сфокусироваться на втором случае (с бидиагональной $B$), потому что в нем проблема формулируется в итоге достаточно красиво: используя бидиагональность $B$, мы можем матрицей перестановки $P$ (которая идейно делает то же самое, что и Faro shuffle с картами) поменять местами столбцы и строки $\begin{bmatrix} 0 & B\ B^T & 0 \end{bmatrix}$, чтобы свести ее к так называемой форме Голуба-Кахана:

T_{\mathrm{GK}} = P
\begin{bmatrix}
0 & B\\
B^T & 0
\end{bmatrix}
P^T = 
\begin{bmatrix}
0 & \alpha_1 & 0 & & & & \\
\alpha_1 & 0 & \beta_1 & & & & \\
 & \beta_1 & 0 & \alpha_2 & & & \\
  &  & \alpha_2 & 0 & & & \\
 & &  & & \ddots & & \\
  & &  & & & & \alpha_n\\
    & &  & & & \alpha_n & 0\\
\end{bmatrix}

($\alpha_i$ и $\beta_i$ - элементы исходной бидиагональной $B$ с главной и над-диагонали соответственно). Более подробно про эту конструкцию описано здесь (стр. 212-213). В итоге у нас получается симметричная тридиагональная матрица с нулями на главной диагонали (проще уже особо не бывает, наверное). Есть устойчивые методы решения тридиагональных систем уравнений, например, метод прогонки. Предполагая, что мы хорошо знаем сингулярные значения $B$ (а они связаны с собственными значениями $T_{\mathrm{GK}}$), мы можем попробовать восстановить собственные векторы $T_{\mathrm{GK}}$ (и, соответственно, сингулярные векторы $B$) просто решая систему уравнений $\lvert T_{\mathrm{GK}} - \sigma_i I \rvert = 0$. Может сломаться (и скорее всего сломается) ортогональность собственных векторов, но тогда мы уже подумаем, как это можно исправить.

"Продвинутые" методы SVD

Нам на данный момент известно о нескольких альтернативных подходах к вычислению SVD, сразу с ортогональными векторами:

1. Вариация на тему dqds, которая считает еще и сингулярные векторы;
2. MR3 алгоритм (начиная со стр. 124), который одновременно применяется к $B^TB$, $BB^T$ и форме Голуба-Кахана $T_\mathrm{GK}$ (см. выше).

Это достаточно сложные алгоритмы, я не уверен, реализованы ли они где-нибудь. Так как реализовывать их самостоятельно, скорее всего, достаточно сложно, то эти алгоритмы по приоритету меньше остальных идей. Тем не менее, неплохой идеей может быть посмотреть их идею, поискать их реализацию или возможно даже попробовать сделать самим, если что-то окажется не таким уж и страшным и затратным (в японском тезисе про dqds/oqds есть аккуратно написанные в псевдокоде алгоритмы).