rowSums(X %*% C * X) 的语义是什么
What's the semantics of 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
现在问题变成了为什么它们是等价的？我想需要做一些巧妙的矩阵划分才能理解等价背后的基本原理，但我还不够开明。
就像而不是
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).
因此，在这种情况下，没有改变操作顺序、分区或其他任何重要技巧；它只是利用现有的矩阵运算，为我们每个 row/column 重复相同的动作。
额外：
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
如果 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]
我正在尝试理解函数 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
现在问题变成了为什么它们是等价的？我想需要做一些巧妙的矩阵划分才能理解等价背后的基本原理，但我还不够开明。
就像而不是
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).
因此，在这种情况下，没有改变操作顺序、分区或其他任何重要技巧；它只是利用现有的矩阵运算，为我们每个 row/column 重复相同的动作。
额外：
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
如果 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]