Sum of correlation matrix convergence

Question

Consider a correlation matrix P but with diagonal elements set to zero. I want to determine the minimum order n for which the partial sum diag(3) + P + P %^% 2 + P %^% 3 + ... + P %^% n converges in L1 norm within some tolerance. I looked into this related question, but it does not apply in my case, since it doesn't keep the powers of P or sum them up.

I can do this in a really lengthy and lousy way with for loops, but I don't want to, since I have a big data frame with many time windows. I'm looking for something efficient.

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

To sum the powers of P, I did:

library(expm)
matrix(mapply(sum, diag(3), P, P %^% 2, P %^% 3, MoreArgs = list(na.rm = TRUE)), ncol = 3)

Note that the %^% operator is from package expm.

Answer 1

x %^% n computes the nth power of x efficiently, but it is inefficient to compute x %^% i for all i from 0 to n, because each x %^% i requires O(log(i)) matrix multiplications.

In general, the most efficient way to compute all of the powers of x up to the nth is recursive multiplication by x, possibly taking advantage of the diagonalizability of x.

The difference is nontrivial for large n: whereas

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

requires O(log(n!)) = O(n * log(n)) matrix multiplications,

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

requires just O(n).

Here is a function that recursively computes the powers of a matrix x and their sum until it encounters a power whose 1-norm is less than tol. It begins by checking that the spectral radius of x is less than 1, which is a necessary and sufficient condition for convergence of the norm of x %^% n to 0 and thus a necessary condition for convergence of the power series. It does not attempt to diagonalize x, which would simplify computation of the power series but complicate computation of norms.

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)
}

Your matrix P has spectral radius greater than 1, hence:

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

We can always construct a matrix P whose spectral radius is less than 1, for the purpose of testing 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

Hence convergence is attained at n = 60. You can check that the reported sum is correct by comparing against the directly (but inefficiently) calculated value:

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

And just for fun:

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