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

Question

我想了解功能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-by-p矩阵，而cov是p-by-p矩阵。它们的内容对于这个问题的目的无关紧要。

根据该文件，mahalanobis计算x中所有行的Mahalanobis平方距离。我把它作为一个提示，并找到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

现在问题变成为什么它们相同？我想人们需要做一些聪明的矩阵分区才能理解等价背后的基本原理，但我还没有开明到足以看到它。

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

Answer 2

如果A和B是任何两个适形矩阵，并且a和b是任何两个相同长度的矢量，我们将使用这些事实。第一行说A * B的第一行等于B的第一行A的第一行。第二行说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]