相关矩阵收敛之和

Question

假设相关矩阵 P 的对角线为零。我想确定阶数 n，其中所有相关矩阵阶数的总和会收敛，即 diag(3)+ P + P%^%2 + P%^%3 + ... + P%^%n 会收敛，这意味着 L1 范数下降到 tol 以下。我查看了但这对我没有帮助，因为它既不保留订单，也不对它们进行总结。我可以使用 for 循环以一种非常冗长和糟糕的方式完成它，但我不想这样做，因为我有一个很大的 df 有很多时间 windows，所以我正在寻找有效的东西。谢谢！

P <- matrix(c(0,0.1,0.8,0.1,0,-0.7,0.8,-0.7,0), nrow = 3, ncol = 3, byrow = TRUE)

一些注意事项：%^% 运算符来自 expm 包。对我使用的矩阵求和 matrix(mapply(sum, diag(3), P, P%^%2, P%^%3, MoreArgs=list(na.rm=T)), ncol=3).

Answer 1

x %^% n 有效地计算了 x 的 n 次方，但是计算 x %^% i 是 低效的 i 从 0 到 n，因为每个 x %^% i 需要 O(log(i)) 次矩阵乘法。

一般来说，计算 x 直到第 n 次方的所有幂的最有效方法是递归乘以 x，可能利用 diagonalizability 共 x。

对于大 n 来说，差异是不小的：而

x2 <- x %^% 2
x3 <- x %^% 3
x4 <- x %^% 4
## and so on

需要O(log(n!)) = O(n * log(n))次矩阵乘法，

x2 <- x  %*% x
x3 <- x2 %*% x
x4 <- x3 %*% x
## and so on

只需要 O(n).

这是一个递归计算矩阵 x 的幂及其总和的函数，直到遇到 1-范数小于 tol 的幂为止。它首先检查 x 的 spectral radius 是否小于 1，这是 x %^% n 范数收敛到 0 的充分必要条件，因此也是收敛的必要条件电源系列。它不会尝试对角化x，这会简化幂级数的计算，但会使范数的计算复杂化。

f <- function(x, tol = 1e-06, nmax = 1e+03) {
    stopifnot(max(abs(eigen(x, only.values = TRUE)$values)) < 1)
    pow <- sum <- diag(nrow(x))
    nrm <- rep.int(NA_real_, nmax + 1)
    i <- 1
    while ((nrm[i] <- norm(pow, "1")) >= tol && i <= nmax) {
        pow <- pow %*% x
        sum <- sum + pow
        i <- i + 1
    }
    list(x = x, tol = tol, nmax = nmax, n = i - 1, sum = sum, 
         norm = nrm[seq_len(i)], converged = nrm[i] < tol)
}

您的矩阵 P 的光谱半径大于 1，因此：

P <- matrix(c(0, 0.1, 0.8, 0.1, 0, -0.7, 0.8, -0.7, 0), 3L, 3L, byrow = TRUE)
f(P)

Error in f(P) : 
  max(abs(eigen(x, only.values = TRUE)$values)) < 1 is not TRUE

我们总可以构造一个谱半径小于1的矩阵P，用于测试f:

set.seed(1L)
m <- 3L
V <- matrix(rnorm(m * m), m, m)
D <- diag(runif(m, -0.9, 0.9))
P <- V %*% D %*% solve(V)
all.equal(sort(eigen(P)$values), sort(diag(D))) # [1] TRUE
(fP <- f(P))

$x
            [,1]      [,2]       [,3]
[1,]  0.26445172 0.5317116 -0.2432849
[2,]  0.04932194 0.6332122  0.1496390
[3,] -0.31174920 0.6847937  0.1682702

$tol
[1] 1e-06

$nmax
[1] 1000

$n
[1] 60

$sum
            [,1]     [,2]        [,3]
[1,]  1.53006915 2.081717 -0.07302465
[2,] -0.04249899 4.047528  0.74063387
[3,] -0.60849191 2.552208  1.83947562

$norm
 [1] 1.000000e+00 1.849717e+00 1.223442e+00 1.008928e+00 7.799426e-01
 [6] 6.131516e-01 4.795602e-01 3.754905e-01 2.938577e-01 2.299751e-01
[11] 1.799651e-01 1.408263e-01 1.101966e-01 8.622768e-02 6.747162e-02
[16] 5.279503e-02 4.131077e-02 3.232455e-02 2.529304e-02 1.979107e-02
[21] 1.548592e-02 1.211727e-02 9.481396e-03 7.418905e-03 5.805067e-03
[26] 4.542288e-03 3.554202e-03 2.781054e-03 2.176090e-03 1.702724e-03
[31] 1.332329e-03 1.042507e-03 8.157298e-04 6.382837e-04 4.994374e-04
[36] 3.907945e-04 3.057848e-04 2.392672e-04 1.872193e-04 1.464934e-04
[41] 1.146266e-04 8.969179e-05 7.018108e-05 5.491455e-05 4.296896e-05
[46] 3.362189e-05 2.630810e-05 2.058529e-05 1.610736e-05 1.260351e-05
[51] 9.861865e-06 7.716607e-06 6.038009e-06 4.724558e-06 3.696822e-06
[56] 2.892650e-06 2.263410e-06 1.771049e-06 1.385792e-06 1.084340e-06
[61] 8.484627e-07

$converged
[1] TRUE

因此在 n = 60 处达到收敛。您可以通过与直接（但效率低下）计算的值进行比较来检查报告的总和是否正确：

library("expm")
all.equal(Reduce(`+`, lapply(0:fP$n, function(i) P %^% i)), fP$sum) # [1] TRUE

并且只是为了好玩：

plot(0:fP$n, fP$norm)

相关矩阵收敛之和

Sum of correlation matrix convergence

r

matrix

sum

convergence