рандомизированные сингулярные значения SVD

Question

рандомизированное SVD разбивает матрицу, извлекая первые k сингулярных значений / векторов, используя k + p случайных проекций.это работает на удивление хорошо для больших матриц.

Мой вопрос касается сингулярных значений, которые выводятся из алгоритма.почему значения не равны первым k-единичным значениям, если вы делаете полный SVD?

Ниже у меня есть простая реализация в R. Любые предложения по улучшению производительности будут оценены.

rsvd = function(A, k=10, p=5) {
       n = nrow(A)
       y = A %*% matrix(rnorm(n * (k+p)), nrow=n)
       q = qr.Q(qr(y))
       b = t(q) %*% A
       svd = svd(b)
       list(u=q %*% svd$u, d=svd$d, v=svd$v)
       }

> set.seed(10)

> A <- matrix(rnorm(500*500),500,500)

> svd(A)$d[1:15]
 [1] 44.94307 44.48235 43.78984 43.44626 43.27146 43.15066 42.79720 42.54440 42.27439 42.21873 41.79763 41.51349 41.48338 41.35024 41.18068

> rsvd.o(A,10,5)$d
 [1] 34.83741 33.83411 33.09522 32.65761 32.34326 31.80868 31.38253 30.96395 30.79063 30.34387 30.04538 29.56061 29.24128 29.12612 27.61804

Answer 1

Расчет

Я считаю, что ваш алгоритм является модификацией алгоритма Martinsson et al. .Если я правильно понял, это особенно важно для приближений для матриц низкого ранга.Я могу ошибаться, хотя.

Разницу легко объяснить огромной разницей между фактическим рангом A (500) и значениями k (10) и p (5).Кроме того, Мартинссон и др. Отмечают, что значение для p на самом деле должно быть больше, чем выбранное значение для k.

Поэтому, если мы применим ваше решение с учетом их соображений, используя:

set.seed(10)
A <- matrix(rnorm(500*500),500,500) # rank 500
B <- matrix(rnorm(500*50),500,500)  # rank 50

Мы считаем, что использование большего значения p по-прежнему приводит к огромному ускорению по сравнению с исходным алгоритмом SVD.

> system.time(t1 <- svd(A)$d[1:5])
   user  system elapsed 
    0.8     0.0     0.8 

> system.time(t2 <- rsvd(A,10,5)$d[1:5])
   user  system elapsed 
   0.01    0.00    0.02 

> system.time(t3 <- rsvd(A,10,30)$d[1:5])
   user  system elapsed 
   0.04    0.00    0.03 

> system.time(t4 <- svd(B)$d[1:5]       )
   user  system elapsed 
   0.55    0.00    0.55 

> system.time(t5 <-rsvd(B,10,5)$d[1:5]  )
   user  system elapsed 
   0.02    0.00    0.02 

> system.time(t6 <-rsvd(B,10,30)$d[1:5]  )
   user  system elapsed 
   0.05    0.00    0.05 

> system.time(t7 <-rsvd(B,25,30)$d[1:5]  )
   user  system elapsed 
   0.06    0.00    0.06 

Но мы видим, что использование более высокого p для более низкого рангаматрица действительно дает лучшее приближение.Если мы позволим k также приблизиться к рангу немного ближе, разница между реальным решением и приближением станет ок.0, в то время как прирост скорости все еще существенен.

> round(mean(t2/t1),2)
[1] 0.77

> round(mean(t3/t1),2)
[1] 0.82

> round(mean(t5/t4),2)
[1] 0.92

> round(mean(t6/t4),2)
[1] 0.97

> round(mean(t7/t4),2)
[1] 1

Итак, в общем, я считаю, что можно сделать вывод, что:

• p следует выбирать так, чтобы p> k (Мартинссон вызываетэто l если я прав)
• k не должно сильно отличаться от ранга (A)
• Для матриц низкого ранга результат обычно лучше.

---

Оптимизация:

Насколько я понимаю, это отличный способ сделать это.Я действительно не мог найти более оптимальный путь.Единственное, что я могу сказать, это то, что конструкция t(q) %*% A не рекомендуется.Для этого нужно использовать crossprod(q,A), что должно быть немного быстрее.Но в вашем примере разницы не было.

Answer 2

В работе Halko, Martinsson и Tropp также рекомендуется выполнить несколько мощных итераций, прежде чем вычислять QR. По умолчанию мы выполняем 3 мощных итерации в реализации в scikit-learn , и мы обнаружили, что на практике очень хорошо работает .