Посчитать правые сингулярные векторы у стоячей матрицы

[ilghiz]

Добрый день,

пусть у нас есть стоячая матрица $A \in {\mathbb R}^{N \times M}$, $N>>M$, заданная с машинной точностью $\epsilon$, и мне надо найти правые сингулярные векторы $V$ и сингулярные значения $\diag(\sigma_1,\dots,\sigma_M)=D$, сингулярного разложения $A = UDV^T$.

Неприятность в том, что памяти есть только на несколько копий $M^2$, а вот всю $A$ в память поместить не получается, да и $N$ на столько много больше $M$, что могут набегать ошибки округления. Матрица $A$ подается по строково, и вернуться и прочитать предыдущую строку нет. Можно хранить несколько строк, порядка $M$ в памяти.

Наивный способ решения строить $A^T A$ и считать $V$ как собственные вектора, а из собстенных значений извлекая из них корень получать сингулярные числа.

В реальности это приводит к тому, что:

$N = 2^{20}, M=10$ точность элементов $A$ - 24 бита, и для надежного вычисления $A^T A$ нам необходимо строить ее с точностью $24*2+20 = 68$ бит. Не мало, к сожалению, но ничего не поделать, в double-double влезет, но уже lapack не вызовешь.

Альтернативно была идея делать что-то типа ранг-ревеалинга у $A^T = QRP$, ($P$ - матрица перестановок от rank-revealing) и предполагать, что $QR R^T Q^T$, вернее $R R^T$ будет строиться как-то лучше, но что-то подсказывает, что это подход в никуда.

Скажите, пожалуйста, есть ли какой-то красивый способ найти правые сингулярные вектора не возведя в квадрат обусловленность матрицы?

Спасибо!

[Евгений Машеров]

0. Насколько действительно серьёзна эта проблема? Всё-таки "число обусловленности" это скорее к решению уравнений относится. Ну и закладываться на наихудший случай можно, но слишком накладно. Опять же - 24 бита это мантисса типа single? Может, просто использовать немудрёный double? Или это точность входных данных, определяемая, скажем, применяемым АЦП? В общем, нужно ли сразу паниковать, или "и так сойдёт"?
1. Можно ли повторно получить матрицу А ("считаем на лету" или "данные на магнитной ленте, прочитали до конца и можно перемотать")? В последнем случае можно попробовать итеративно уточнять. Скажем, имея оценку V и D, каждую i-тую строку $A_i$ домножать слева $B_i=A_iD^{-1}V^T$, и затем использовать для расчёта $B^TB$, разложение которой даст поправки к V и D.
2. Если можно рассматривать первые M строк, как "типовые", можно по ним рассчитать разложение и его использовать для коррекции последующих строк.
3. Для "не возведя в квадрат", боюсь, понадобится произвольный доступ ко всем элементам А. Скажем, описанный в древнем ("...на языке Algol-60") справочнике Уилкинсона и Райнша алгоритм через переход к двухдиагональной матрице и затем QR-разложение её. Я не вижу способа реализовать при указанном способе доступа (но не присягну, что его вовсе нет, хотя очень сильно это подозреваю).

[ilghiz]

Спасибо большое Евгений Машеров, за ответы!


0. ресурсов не хватает, поэтому стала серьезной.
1. Данные приходят с АЦП, имеют 24 бита (в скором времени будут 32 бита иметь), и, перемотать их нельзя.

Цитата:

2. Если можно рассматривать первые M строк, как "типовые"

пока делал именно так, или уменьшал мантиссу у сумматора.

Данных может быть много, вплоть до $2^{30}$, в этом случае, и если мантисса 32 бита, надо копить на 94 бита.

Накапливание происходит на FPGA (у нее очень мало памяти), считать можно уже в обычном компьютере, где ограничения по памяти почти нет, но собирать данные и сразу передавать в компьютер не получается, трафик большой. Грубо говоря, у меня за секунду приходит $N=2\times10^8$ и $M=200$ по 24 бита, то есть десятки гигабайт, то есть адекватным способом забрать их в компьютер нельзя.

В то же время, $A^T A$ для $M=200$ по 94 бит - это 240КБайт - что есть очень обозримо.

При внимательном осмотре данных, вернее переводя собственные числа $A^T A$ в сингулярные $A$, в спектре сингулярных чисел $A$ у меня имеется разброс больше мантиссы исходных чисел, но меньше мантиссы + число отсчетов. То есть пусть мантисса 24 бита, и 1 миллион отсчетов. В спектре сингулярных чисел $A$ у меня имеется важные для меня сингулярные числа, которые меньше главного примерно в $2^{28}$ раз.

Это все приводит к тому, что 60 бит - это минимум, с которым я долен копить $A^T A$, хотя надо реально точнее, то есть фактически по полной $20+24 \times 2=68$ и даже на обычном процессоре, когда я считаю задачу с собственными значениями для $A^T A$, надо как-то переходить к 4-ной точности, и так просто с лапака это не сделать. Приходится разлагать $A^T A$ Холецким с ведущим элементом, вручную его запрограммировав в quad-double, ожидая, что такое разложение уменьшит обусловленность не в корень, но хоть во сколько-то, а потом у каждого фактора считать сингулярное разложение.

Если у меня будет 94 бита, то тут уже может и точности не хватить. Поэтому есть непреодолимое желание не возводить обусловленность в квадрат еще на процессе получения данных, ибо даже $N=2^{30}$ и 32 битная точность исходных данных - это всего-то 62 бита, и похоже у этой матрицы важные для меня сингулярные числа будут только в $2^{35}$-$2^{40}$ раз меньше главного.