Быстрое определение, является ли плотная матрица низкого ранга

13

В программном проекте, над которым я работаю, некоторые вычисления намного проще для плотных матриц низкого ранга. В некоторых проблемных случаях используются плотные матрицы низкого ранга, но они даны мне полностью, а не как факторы, поэтому мне придется проверять ранг и фактор матрицы, если я хочу воспользоваться преимуществами структуры низкого ранга ,

Матрицы, о которых идет речь, обычно являются полностью или почти полностью плотными, причем n колеблется от ста до нескольких тысяч. Если матрица имеет низкий ранг (скажем, менее 5-10), то вычисление SVD и использование его для формирования факторинга низкого ранга стоит усилий. Однако, если матрица не низкого ранга, то усилия будут потрачены впустую.

Таким образом, я хотел бы найти быстрый и достаточно надежный способ определения, является ли ранг низким, прежде чем вкладывать усилия в полную факторизацию SVD. Если в какой-то момент становится ясно, что ранг выше порога, процесс может быть немедленно остановлен. Если процедура ошибочно объявляет матрицу низкого ранга, когда это не так, это не является большой проблемой, так как я все еще буду делать полный SVD для подтверждения низкого ранга и нахождения факторизации низкого ранга.

Варианты, которые я рассмотрел, включают ранг, раскрывающий LU или QR-факторизацию, за которым следует полный SVD в качестве проверки. Есть ли другие подходы, которые я должен рассмотреть?

Брайан Борхерс
источник

Ответы:

8

Из этой статьи я недавно узнал об одном уловке . Вы начинаете делать QR-код, раскрывающий ранг, и останавливаетесь после первых отражений Домохозяина, когда у вас есть матрица вида [ R 1 R 12 0 R 22 ] , где R 1 треугольник размера k × k , а R 22 обычно не треугольный (так как мы остановились после первых k итераций нашего основного цикла). На данный момент, вы проверяете , если | | R 22 | | & le ; & epsi ; : если оно выполнено, тоk

[R1R120R22],
R1k×kR22kR22εAнаходится на расстоянии не более от матрицы ранга k ; в противном случае это не должно быть (исключая ошибки с цифрами).εk

Эта процедура стоит для плотной матрицы n × n .O(n2k)n×n

Федерико Полони
источник
По сути, это тот подход, который я описал в этом вопросе. Я думаю, что предложенный ответ Вольфганга Бангерта мог бы быть лучше, чем . O(n2k)
Брайан Борчерс
7

Проблема, конечно, заключается в том, что вычисление истинного ранга (например, посредством разложения QR) на самом деле не дешевле, чем вычисление низкосортного представления матрицы.

Лучшее, что вы можете сделать, это использовать рандомизированный алгоритм для поиска аппроксимаций низкого ранга. Они могут, по крайней мере, теоретически, быть значительно быстрее, чем работать со всей матрицей, потому что, по сути, они только вычисляют разложения для проекций матрицы на случайные подпространства.

Хороший вопрос, стоит ли это для матрицы размером , но если ваши проблемы действительно станут большими, я подозреваю, что это окупится.100×100

Вольфганг Бангерт
источник
Из того, что я знаю об этих алгоритмах, они производят матрицу низкого ранга, которая достаточно близка по норме к данной матрице. Мне нужно знать, существует ли (например) матрица ранга 10 или меньше, которая очень близка к данной матрице (скажем, относительная ошибка 1,0e-10 или лучше.)
Брайан Борхерс
Да, но вы также можете выполнить QR-декомпозицию спроецированной (низкоразмерной) матрицы, и если эта декомпозиция обнаружит отсутствие полного ранга, то у вас также будет оригинальная матрица с дефицитом ранга. Разве это не тот критерий, который нужен для разложения QR-кода на исходную матрицу?
Вольфганг Бангерт
kkk1kAkknO(k2n)AO(kn2)
k=nkkn2n3
1

Другой подход, который стоит попробовать, - использовать адаптивную перекрестную аппроксимацию (ACA). Это довольно популярный алгоритм, который имеет много реализаций, доступных онлайн. Для справки, вы можете увидеть оригинал статьи:

ACA и ее вариации (скажем, ACA +, гибридная перекрестная аппроксимация HCA) могут использоваться в различных сценариях. Вы уже рассчитали всю плотную матрицу, что является одним из преимуществ, так как вы сможете точно рассчитать остатки, если это необходимо.

O(Nr)Nr(ϵ)rϵO(N2r)

Антон Меньшов
источник
0

Для простого случая, когда матрица A является симметричным положительно определенным, вычислите его, скажем, 20 самых больших собственных значений, и посмотрите, если они 0, or compare norms. ARPACK is fast for this; more important, it needs only a function ИксAИкс, Так что для общегоAпосмотрите на собственные значения ATA (как LinOp, не создавая его.)

scipy.sparse.linalg.svds делает это: LinOp(ATA) Arpack, для A любого размера:

from scipy.sparse.linalg import svds
sing = svds( A, k=20, tol=1e-4, return_singular_vectors=False )  # v0=random
# runtimes on random-normal n x n:
# n = 100, 1k, 2k
#       5, 130, 770 ms
Денис
источник