如何简单高效地编写二次型代码

Question

我正在尝试编写二次型 Z'(S)^{-1} Z

代码如下

z <- matrix(rnorm(200 * 100), 200, 100)
S <- cov(z)
quad.naive <- function(z, S) {
        Sinv <- solve(S)
        rowSums((z %*% Sinv) * z) 
}

但是，我不确定我是否完全理解函数的最后一行

rowSums((z %*% Sinv) * z)

因为天真，我们应该输入与数学公式完全相同的内容

t(Z) %*% Sinv %*% Z

所以，任何人都可以解释为什么行总和形式与朴素的数学形式相同，尤其是。为什么在两个度量（z 和 Sin）相乘后，它使用逐元素相乘符号 * 乘以 Z，而不是使用 %*%。

(z %*% Sinv) * z

Answer 1

下面的评论有点太长了。

"I'm trying to code a quadratic form Z'(S)^{-1} Z"我觉得二次型不对

假设 Z 是一个 m x n 矩阵。那么：

S = cov(Z) 是一个 n x n 矩阵
S^-1 是一个 n x n 矩阵
t(Z) 是一个 n x m 矩阵

因此 Z' S^-1 Z（在 R 中：t(Z) %*% solve(S) %*% Z）将意味着将矩阵与以下维度相乘

(n x m) (n x m) (m x n)

这显然行不通。

也许你的意思是 Z %*% solve(S) %*% t(Z) 其中 returns 一个 m x m 矩阵，其对角线与 rowSums(Z %*% Sinv * Z).

相同

更根本的是：quadratic form 不应该是标量吗？或者你在谈论不同的二次形式？

Answer 2

好的，根据我们在评论中的交流和 link 你给本书相关部分的内容 Advanced Statistical Computing 我认为我明白了问题是。

我post这是一个单独的（真实的）答案，以避免混淆未来可能想要阅读评论中的思路的读者。

让我们 return 使用您 post 中给出的代码（从 1.3.3 多元正态分布部分复制）

set.seed(2017-07-13)
z <- matrix(rnorm(200 * 100), 200, 100)
S <- cov(z)
quad.naive <- function(z, S) {
        Sinv <- solve(S)
        rowSums((z %*% Sinv) * z)
}

考虑到二次型定义为随机p×1列向量的标量z' Sigma^-1 z（或者R语言t(z) %*% solve(Sigma) %*% z），可能会出现两个问题：

为什么 z 作为 matrix 给出（而不是书中所述的 p 维列向量），并且
在quad.naive中使用rowSums的原因是什么？

首先，请记住，二次型是单个随机多变量样本的标量。 quad.naive 实际上 returning 是多元样本（复数！）中二次型的分布。 z 这里包含来自 p = 100 维法线的 200 个样本。

则S是100×100的协方差矩阵，solve(S)return是S的逆矩阵。数量 z %*% Sinv * z（由于 R 的运算符优先级，附加括号不是必需的）returns z 的每个样本的 t(z) %*% solve(Sigma) %*% z 的对角线元素作为矩阵中的行向量.然后采用 rowSums 与采用迹线相同（即具有二次形式 return 每个样本的标量）。另请注意，使用 diag(z %*% solve(Sigma) %*% t(z)) 会得到相同的结果，但在 quad.naive 中，我们避免了双矩阵乘法和附加转置。

还有一个更基本的问题：为什么要看二次型的分布？可以证明，标准正态变量中某些二次型的分布服从卡方分布（参见 Mathai and Provost, Quadratic Forms in Random Variables: Theory and Applications and Normal distribution - Quadratic forms）

具体来说，我们可以证明 p × 1 列向量的二次型 (x - μ)' Σ^-1 (x - μ) 是自由度为 p 的卡方分布。

为了说明这一点，让我们从双变量标准正态中抽取 100 个样本，并计算每个样本的二次形式。

set.seed(2020)
nSamples <- 100
z <- matrix(rnorm(nSamples * 2), nSamples, 2)
S <- cov(z)
Sinv <- solve(S)
dquadform <- rowSums(z %*% Sinv * z)

我们可以将分布可视化为直方图，并叠加 2 个自由度的理论卡方密度。

library(ggplot2)
bw = 0.2
ggplot(data.frame(x = dquadform), aes(x)) +
    geom_histogram(binwidth = bw) +
    stat_function(fun = function(x) dchisq(x, df = 2) * nSamples * bw)

最后，Kolmogorov-Smirnov 检验将二次型分布与具有 2 个自由度的累积卡方分布进行比较的结果导致我们未能拒绝零假设（两个分布相等） .

ks.test(dquadform, pchisq, df = 2)
#
#   One-sample Kolmogorov-Smirnov test
#
#data:  dquadform
#D = 0.063395, p-value = 0.8164
#alternative hypothesis: two-sided

如何简单高效地编写二次型代码

How to code quadratic form both naively and efficiently

r

quadratic