rowSums（X％*％C * X）的语义是什么

Question

I am trying to understand the function stats::mahalanobis. Here is its source, but please just focus on last line, or more specifically, the rowSums(x %*% cov * x) part.

> 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))
}

Here x is a n-by-p matrix, whereas cov is a p-by-p matrix. Their content doesn't matter for the purpose of this question.

According to the document, mahalanobis calculates the squared Mahalanobis distance of all rows in x. I took this as a hint and found a counterpart of rowSums(X %*% C * X) with apply. (it's perfectly fine if you have no idea what I'm talking about; this paragraph merely serves as an explanation of how I came up with the apply form)

> 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

现在的问题变成了为什么它们相等？我想需要做一些聪明的矩阵划分来理解等效性的原理，但是我不足以了解它。

Answer 1

就像代替

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 * X1）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