计算二次形式的有效方法：避免 for 循环？答案

【问题标题】：Efficient way of calculating quadratic forms: avoid for loops?计算二次形式的有效方法：避免 for 循环？
【发布时间】：2025-12-24 19:10:12
【问题描述】：

我想计算 N（N 很大）quadratic forms。我正在使用 R 包“模拟器”中的命令“quad.form”。如何在不使用 for 循环的情况下实现这一点？

到目前为止，我正在使用

library(emulator)

A = matrix(1,ncol=5,nrow=5) # A matrix

x = matrix(1:25,ncol=5,nrow=5) # The vectors of interest in the QF

# for loop
QF = vector()
for(i in 1:5){
QF[i] = quad.form(A,x[,i])
}

有没有更直接有效的方法来计算这些二次型？

有趣的是

quad.form(A,x)

比 for 循环快（10 倍），但我只需要这个结果的对角线。因此，它仍然是计算 N 个感兴趣的二次形式的低效方法。

【问题讨论】：

标签： r quadratic

【解决方案1】：

怎么样

colSums(x * (A %*% x))

？至少得到这个例子的正确答案......而且应该更快！

library("rbenchmark")
A <- matrix(1, ncol=500, nrow=500)
x <- matrix(1:25, ncol=500, nrow=500)

library("emulator")
aa <- function(A,x) apply(x, 2, function (y) quad.form(A,y))
cs <- function(A,x) colSums(x * (A %*% x))
dq <- function(A,x) diag(quad.form(A,x))
all.equal(cs(A,x),dq(A,x))  ## TRUE
all.equal(cs(A,x),aa(A,x))  ## TRUE
benchmark(aa(A,x),
          cs(A,x),
          dq(A,x))
##       test replications elapsed relative user.self sys.self
## 1 aa(A, x)          100  13.121    1.346    13.085    0.024
## 2 cs(A, x)          100   9.746    1.000     9.521    0.224
## 3 dq(A, x)          100  26.369    2.706    25.773    0.592

【讨论】：

它并没有比apply快多少... ;)
同意。让我有点吃惊。
这可能已经足够快了，但如果它太慢，你可以（1）尝试Rcpp解决方案或（2）尝试将问题并行化成块......
谢谢。第一个命令在我的示例中运行良好（A 是 5x5，x 是 5x2000）并且确实快得多。不过，存在一些（可忽略不计的）数值差异。

【解决方案2】：

使用apply函数：

apply(x, 2, function (y) quad.form(A,y))

如果您使矩阵更大 (500x500)，很明显使用 apply 的速度大约是使用 quad.form(A,x) 的两倍：

A <- matrix(1, ncol=500, nrow=500)
x <- matrix(1:25, ncol=500, nrow=500)

system.time(apply(x, 2, function (y) quad.form(A,y)))
# user  system elapsed 
# 0.183   0.000   0.183 

system.time(quad.form(A,x))
# user  system elapsed 
# 0.314   0.000   0.314

编辑

而@Ben Bolker 的回答比apply 快1/3 左右：

system.time(colSums(x * (A %*% x)))
# user  system elapsed 
# 0.123   0.000   0.123

【讨论】：

非常感谢您的回答。这是 Bolker 答案的一个有趣的替代方案。在我的真实场景中，您的方法的效率介于 Ben Bolker 的方法和我在问题末尾描述的方法之间。