Какова семантика rowSums (X% *% C * X) - PullRequest
Новая
       41

Какова семантика rowSums (X% *% C * X)

2 голосов
/ 17 марта 2019

Я пытаюсь понять функцию 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

Теперь возникает вопрос, почему они эквивалентны? Я полагаю, что нужно сделать какой-нибудь умный матричный раздел, чтобы понять причину эквивалентности, но я не достаточно образован, чтобы это увидеть.

Ответы [ 2 ]

5 голосов
/ 17 марта 2019

Так же, как вместо 

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
3 голосов
/ 17 марта 2019

Если 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]
...