Какова семантика rowSums(X %*% C * X)
Я пытаюсь понять функцию stats::mahalanobis. Вот его источник, но, пожалуйста, сосредоточьтесь на последней строке или, точнее, на части rowSums(x %*% cov * x).
> mahalanobis
function (x, center, cov, inverted = FALSE, ...)
{
x <- if (is.vector(x))
matrix(x, ncol = length(x))
else as.matrix(x)
if (!isFALSE(center))
x <- sweep(x, 2L, center)
if (!inverted)
cov <- solve(cov, ...)
setNames(rowSums(x %*% cov * x), rownames(x))
}
Здесь x — матрица размера n на p, тогда как cov — матрица размера p на p. Их содержание не имеет значения для целей этого вопроса.
Согласно документу, mahalanobis вычисляет квадрат расстояния Махаланобиса для всех строк в x. Я воспринял это как подсказку и нашел аналог rowSums(X %*% C * X) с apply. (это совершенно нормально, если вы понятия не имеете, о чем я говорю; этот абзац просто служит объяснением того, как я придумал форму apply)
> X = matrix(rnorm(1:6), nrow = 3)
> C = matrix(rnorm(1:4), nrow = 2)
> rowSums(X %*% C * X)
[1] -0.03377298 0.49306538 -0.16615078
> apply(X, 1, function(row) {
+ t(row) %*% C %*% row
+ })
[1] -0.03377298 0.49306538 -0.16615078
Теперь возникает вопрос, почему они эквивалентны? Я предполагаю, что нужно сделать какой-то умный матричный раздел, чтобы понять причину эквивалентности, но я недостаточно просвещен, чтобы увидеть это.
@RuiBarradas Спасибо за ваше предложение! Это вроде как работает, но теперь мой мозг забит всеми этими арифметическими деталями и не может уловить основную идею этой формы. Я бы предпочел ответ с чуть более высокого уровня.
Ответы 2
Так же, как вместо
sapply(1:5, `*`, 2)
# [1] 2 4 6 8 10
или петля, которую мы предпочитаем
1:5 * 2
# [1] 2 4 6 8 10
поскольку это векторизованное решение, делающее то же самое - поэлементное умножение,
rowSums(X %*% C * X)
# [1] 0.2484329 0.5583787 0.2303054
можно увидеть лучше, чем
apply(X, 1, function(row) t(row) %*% C %*% row)
# [1] 0.2484329 0.5583787 0.2303054
причем оба они снова делают то же самое, только первый из них более лаконичен.
В частности, в моем первом примере мы перешли от скаляров к векторам, а теперь идем от векторов к матрицам. Первый,
X %*% C
# [,1] [,2]
# [1,] 0.7611212 0.6519212
# [2,] -0.4293461 0.6905117
# [3,] 1.2917590 -1.2970376
соответствует
apply(X, 1, function(row) t(row) %*% C)
# [,1] [,2] [,3]
# [1,] 0.7611212 -0.4293461 1.291759
# [2,] 0.6519212 0.6905117 -1.297038
Теперь второе произведение в t(row) %*% C %*% row делает две вещи: 1) поэлементное умножение t(row) %*% C и row, 2) суммирование. Точно так же * в X %*% C * X выполняет 1) и rowSums выполняет суммирование, 2).
Так вот, в данном случае нет существенных ухищрений по изменению порядка операций, разбиения или чего-то еще; это просто использование существующих матричных операций, которые повторяют одни и те же действия с каждой строкой/столбцом для нас.
Дополнительный:
library(microbenchmark)
microbenchmark(rowSums = rowSums(X %*% C * X),
apply = apply(X, 1, function(row) t(row) %*% C %*% row),
times = 100000)
# Unit: microseconds
# expr min lq mean median uq max neval cld
# rowSums 3.565 4.488 5.995129 5.117 5.589 4940.691 1e+05 a
# apply 24.126 26.402 32.539559 27.191 28.615 129234.613 1e+05 b
Подождите, значит, семья apply не может воспользоваться преимуществами векторизации? Спасибо за устранение недоразумения для меня!
@nalzok, матричные продукты в t(row) %*% C %*% row можно рассматривать как векторизованную версию использования ручных циклов, поэтому в некотором смысле ваше решение apply имеет некоторую векторизацию. Однако, кроме этого, это не совсем векторизация, если мы думаем об этом с точки зрения скорости, см., например, stackoverflow.com/q/28983292/1320535. Я также добавил сравнение производительности.
Спасибо за разъяснения! Было бы еще лучше, если бы вы могли показать мне, как делается бенчмарк, потому что сокращения lq и cld не имеют для меня особого смысла?
@nalzok, добавил. Я полагаю, что lq и uq обозначают нижний и верхний квантили.
Если A и B — любые две совместимые матрицы, а a и b — любые два вектора одинаковой длины, мы воспользуемся этими фактами. Первый говорит, что первая строка A * B равна первой строке A, умноженной на первую строку B. Вторая говорит, что первая строка A %*% B равна первой строке A все время B. Третий говорит матричное умножение двух векторов может быть выражено как сумма их поэлементного умножения.
(A * B)[i, ] = A[i, ] * B[i, ] by the defintion of elementwise multiplication [1]
(A %*% B)[i, ] = A[i, ] %*% B as taking ith row is same as premultplying by ei [2]
a %*% b = sum(a * b) by definition of %*% [3]
Таким образом мы получаем:
rowSums(X %*% C * X)[i]
= sum((X %*% C * X)[i, ])
= sum((X %*% C)[i, ] * X[i, ]) by [1]
= (X %*% C)[i, ] %*% X[i, ] by [3]
= X[i, ] %*% C %*% X[i, ] by [2]
= apply(X, 1, function(row) t(row) %*% C %*% row)[i]
Я предлагаю создать матрицу меньшего размера, скажем, 2x2
Xи выполнить вычисления вручную. Должно стать очевидным, почему они одинаковы.