rowSums(X %*% C * X) 的语义是什么答案

【问题标题】：What's the semantics of rowSums(X %*% C * X)rowSums(X %*% C * X) 的语义是什么
【发布时间】：2019-08-07 23:21:12
【问题描述】：

我正在尝试理解函数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

现在问题变成了为什么它们是等价的？我想人们需要做一些巧妙的矩阵划分来理解等价背后的基本原理，但我还不够开明。

【问题讨论】：

我的建议是创建一个更小的 2x2 矩阵 X 并手动进行计算。应该很明显为什么它们是相同的。
@RuiBarradas 感谢您的建议！这有点工作，但现在我的大脑充满了所有这些算术细节，无法理解这种形式背后的重要概念。我更希望得到更高层次的答案。

标签： r math matrix vectorization matrix-multiplication

【解决方案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 * 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 代表上下分位数，我相信。

【解决方案2】：

如果 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]

【讨论】：