Когда "crossprod" предпочтительнее "%*%", а когда нет?

Question

Когда именно crossprod(X,Y) предпочтительнее t(X) %*% Y когда X и Y обе матрицы? В документации сказано

Данные матрицы x а также y в качестве аргументов вернуть матричный перекрестный продукт. Это формально эквивалентно (но обычно немного быстрее) t(x) %*% y (crossprod) или же x %*% t(y) (tcrossprod).

Так когда же это не быстрее? При поиске в Интернете я нашел несколько источников, которые утверждали, что crossprod обычно является предпочтительным и должен использоваться по умолчанию (например, здесь), или что это зависит, и результаты сравнений неубедительны. Например, Дуглас Бейтс (2007) говорит, что:

Обратите внимание, что в более новых версиях R и библиотеки BLAS (по состоянию на лето 2007 года) R's %*% способен обнаружить множество нулей в mm и сократить многие операции, и, следовательно, гораздо быстрее для такой разреженной матрицы, чем crossprod который в настоящее время не использует такие оптимизации [...]

Так, когда я должен использовать %*% и когда crossprod?

Answer 1

Я кратко осветил эту проблему в нескольких моих ответах на вопросы статистических вычислений. Последующая часть моего этого ответа является относительно полной. Обратите внимание, что мое обсуждение там, а также следующее действительно предполагают, что вы знаете, что такое BLAS, особенно что такое подпрограммы BLAS уровня 3 dgemm а также dsyrk являются.

Мой ответ здесь предоставит некоторые доказательства и сравнительный анализ, чтобы заверить вас в моем аргументе там. Кроме того, будет дано объяснение комментарию Дугласа Бейтса.

Что crossprod а также "%*%" действительно?

Давайте проверим исходный код для обеих подпрограмм. Исходный код уровня C для базового пакета R в основном находится в

R-release/src/main

В частности, матричные операции определены в

R-release/src/main/array.c

Сейчас,

• "%*%" соответствует рутине C matprod, который вызывает dgemm с transa = "N" а также transb = "N";
• crossprod видно, что это составная функция, соответствующая symcrossprod, crossprod, symtcrossprod а также tcrossprod (аналоги для сложных матриц, как ccrossprod, не перечислены здесь).

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

Профилирование уровня R может выявить эти издержки. Рассмотрим следующий пример:

A <- matrix(ruinf(1000 * 1000), 1000)
B <- matrix(ruinf(1000 * 1000), 1000)

Rprof("brutal.out")
result <- t.default(A) %*% B
Rprof(NULL)
summaryRprof("brutal.out")

Rprof("smart.out")
result <- crossprod(A, B)
Rprof(NULL)
summaryRprof("smart.out")

Обратите внимание, как t.default вступает в результат профилирования для первого случая.

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

Когда значительны затраты на транспозицию?

Позволять A быть k * m матрица и B k * n матрица, затем матричное умножение A'B (' означает транспонирование) имеет счетчик FLOP (количество сложений и умножений с плавающей запятой) 2mnk, Если вы делаете t(A) %*% B, транспонирование накладных расходов mk (количество элементов в A), поэтому соотношение

useful computation : overhead = 2n : 1

Если не n велико, накладные расходы не могут быть амортизированы при фактическом умножении матриц.

Самый крайний случай n = 1 у вас есть матрица из одного столбца (или вектор) для B, Рассмотрим сравнительный анализ:

library(microbenchmark)
A <- matrix(runif(2000 * 2000), 2000)
x <- runif(2000)
microbenchmark(t.default(A) %*% x, crossprod(A, x))

Вы увидите, что crossprod в несколько раз быстрее!

Когда "%*%" структурно уступает?

Как упоминалось в моем связанном ответе (а также в заметках Бейтса с результатами сравнительного анализа), если вы это сделаете A'A, crossprod обязательно будет в 2 раза быстрее.

Что делать, если у вас есть разреженные матрицы?

Наличие разреженных матриц не меняет основного вывода, приведенного выше. R пакет Matrix для настройки разреженных матриц также имеет разреженные методы вычисления "%*%" а также crossprod, Так что вы все еще должны ожидать crossprod быть немного быстрее.

Так что же означает комментарий Бейтса о BLAS "по состоянию на лето 2007 года"?

Это не имеет ничего общего с разреженными матрицами. BLAS строго предназначен для плотной числовой линейной алгебры.

То, к чему это относится, является различием в вариантах цикла, используемых в справочнике Fli Netlib dgemm, Существует два варианта цикла для планирования умножения матрицы на матрицу. op(A) * op(B) (* здесь обозначает матричное умножение, а не поэлементное произведение!), а вариант полностью определяется настройкой транспонирования op(A):

• за op(A) = A' версия "внутренний продукт" используется в самом внутреннем цикле, и в этом случае обнаружение нуля невозможно;
• за op(A) = A, Используется версия "AXPY", и любой ноль в op(B) может быть исключен из расчета.

Теперь подумайте, как R звонит dgemm, Первый случай соответствует crossprod в то время как второй случай "%*%" (так же как tcrossprod).

В этом аспекте, если ваш B матрица имеет много нулей, в то время как она все еще находится в формате плотной матрицы, то t(A) %*% B будет быстрее чем crossprod(A, B)просто потому, что вариант цикла для первого более эффективен.

Самый яркий пример - это когда B является диагональной матрицей или полосовой матрицей. Давайте рассмотрим полосчатую матрицу (на самом деле симметричная трехдиагональная матрица здесь):

B_diag <- diag(runif(1000))
B_subdiag <- rbind(0, cbind(diag(runif(999)), 0))
B <- B_diag + B_subdiag + t(B_subdiag)

Теперь позвольте A все еще быть полной плотной матрицей

A <- matrix(runif(1000 * 1000), 1000)

Тогда сравните

library(microbenchmark)
microbenchmark(t.default(A) %*% B, crossprod(A, B))

Вы увидите, что "%*%" намного быстрее!

Полагаю, лучшим примером является матрица B треугольная матрица Это довольно часто встречается на практике. Треугольная матрица не рассматривается как разреженная матрица (и также не должна храниться как разреженная матрица).

Внимание: если вы используете R с оптимизированной библиотекой BLAS, такой как OpenBLAS или Intel MKL, вы все равно увидите, что crossprod быстрее. Однако это действительно потому, что стратегия блокирования и кэширования в любых оптимизированных библиотеках BLAS разрушает шаблон планирования вариантов цикла, как в эталонном BLAS Netlib F77. В результате, любое "обнаружение нуля" невозможно. Итак, что вы заметите, так это то, что для этого конкретного примера эталонная BLAS F77 даже быстрее, чем оптимизированная BLAS.