Самый быстрый способ вычисления линейных точечных произведений между двумя тощими высокими матрицами в R

Question

Рассмотрим A и B две высокие тощие матрицы размерности 10^8 X 5. т.е.

r=10^8
c=5
A=matrix(runif(r*c,0,1),r,c)
B=matrix(runif(r*c,0,1),r,c)

Я хочу вычислить скалярное произведение каждой строки A с соответствующей строкой B, т.е.

rowSums(A * B)

Но это довольно медленно, и я хотел бы знать, есть ли более быстрый путь.

Answer 1

Это может вас разочаровать, но на уровне R это уже лучшее, что вы можете получить, не написав код на C самостоятельно. Проблема в том, что, делая rowSums(A * B) вы эффективно делаете

C <- A * B
rowSums(C)

Первая строка выполняет полное сканирование трех больших тонких матриц; в то время как вторая строка выполняет полное сканирование 1 большой тонкой матрицы. Таким образом, мы эквивалентно сканируем матрицу высокой тонкости 4 раза (интенсивно использую память).

На самом деле, для такой операции оптимальный алгоритм требует только сканирования n * p Высоко -тонкую матрицу дважды, выполняя скрещивание произведений:

rowsum <- numeric(n)
for j = 1, 2, ... p
  rowsum += A[,i] * B[,i]

Таким образом, мы также избегаем генерации матрицы C, Обратите внимание, что выше это просто поддельный код, а не действительный код R или даже код C. Но идея ясна, и мы хотим запрограммировать это на C.

---

Аналогия с вашей ситуацией - разница в скорости между sum(x * y) а также crossprod(x, y) предполагая x а также y быть большими векторами одинаковой длины.

x <- runif(1e+7)
y <- runif(1e+7)

system.time(sum(x * y))
#   user  system elapsed 
#  0.124   0.032   0.158 

system.time(crossprod(x, y))
#   user  system elapsed 
#  0.036   0.000   0.056 

В первом случае мы сканируем длинный вектор 4 раза, а во втором случае мы сканируем его только дважды.

---

Актуальность в статистических вычислениях

rowSums(A * B) на самом деле эффективная оценка diag(tcrossprod(A, B)) Обычно наблюдается в регрессионных вычислениях, связанных с дисперсией точечного прогнозирования. Например, в обычных линейных квадратах регрессия с тонкими Q матрица из QR-факторизации модельной матрицы, точечная дисперсия подогнанных значений diag(tcrossprod(Q)), который более эффективно вычисляется rowSums(Q ^ 2), Но все же, это по-прежнему не самая быстрая оценка, по причинам, уже объясненным.