如何找到矩阵何时与循环收敛

Question

我得到了一个矩阵：

P <- matrix(c(0, 0, 0, 0.5, 0, 0.5, 0.1, 0.1, 0, 0.4, 0, 0.4, 0, 0.2, 0.2, 0.3, 0, 0.3, 0, 0, 0.3, 0.5, 0, 0.2, 0, 0, 0, 0.4, 0.6, 0, 0, 0, 0, 0, 0.4, 0.6), nrow = 6, ncol = 6, byrow = TRUE)

使用函数mpow、rows_equal、matrices_equal。我想知道P^n 何时收敛，换句话说，n 是什么，当矩阵中的所有行都相等时以及何时P^n = P^(n+1)。

仅通过查看函数i 已经设法推断出在n=19-21 附近矩阵将收敛。

虽然，我想使用循环找到正确的 n。下面是函数mpow、rows_equal 和matrices_equal。我知道它们可以用不同的方式写，但请保持原样。

mpow <- function(P, n, d=4) {
  if (n == 0) diag(nrow(P)))
  else if (n== 1) P
  else P %*% mpow(P, n - 1))
  }

rows_equal <- function(P, d = 4) {
  P_new <- trunc(P * 10^d)
  for (k in 2:nrow(P_new)) {
    if (!all(P_new[1, ] == P_new[k, ])) {
      return(FALSE)}
      }
    return(TRUE)
    }

matrices_equal <- function(A, B, d = 4) {
  A_new <- trunc(A * 10^d)
  B_new <-trunc(B * 10^d)
  if (all(A_new == B_new)) TRUE else FALSE
  }

现在，要编写循环，我们应该按照以下方式进行：

首先像这样创建一个函数：

when_converged <- function(P) {...}

和

for (n in 1:50)

在 t.ex n = 50 时尝试。

虽然我不知道如何正确编写代码来做到这一点，但有人可以帮我吗？

感谢您阅读我的问题。

Answer 1

实际上，一个更好的方法是这样做：

## transition probability matrix
P <- matrix(c(0, 0, 0, 0.5, 0, 0.5, 0.1, 0.1, 0, 0.4, 0, 0.4, 0, 0.2, 0.2, 0.3, 0, 0.3, 0, 0, 0.3, 0.5, 0, 0.2, 0, 0, 0, 0.4, 0.6, 0, 0, 0, 0, 0, 0.4, 0.6), nrow = 6, ncol = 6, byrow = TRUE)

## a function to find stationary distribution
stydis <- function(P, tol = 1e-16) {
  n <- 1; e <- 1
  P0 <- P  ## transition matrix P0
  while(e > tol) {
    P <- P %*% P0  ## resulting matrix P
    e <- max(abs(sweep(P, 2, colMeans(P))))
    n <- n + 1
    }
  cat(paste("convergence after",n,"steps\n"))
  P[1, ]
  }

那么当你调用函数时：

stydis(P)
# convergence after 71 steps
# [1] 0.002590674 0.025906736 0.116580311 0.310880829 0.272020725 0.272020725

函数stydis，基本上是连续做的：

P <- P %*% P0

直到达到P 的收敛。收敛性在数值上由差异矩阵的 L1 范数确定：

sweep(P, 2, colMeans(P))

L1 范数是所有矩阵元素的最大绝对值。当 L1 范数低于1e-16 时，就会发生收敛。

如您所见，收敛需要 71 个步骤。现在，我们可以通过控制tol（容差）来获得更快的“收敛”：

stydis(P, tol = 1e-4)
# convergence after 17 steps
# [1] 0.002589361 0.025898057 0.116564506 0.310881819 0.272068444 0.271997814

但如果你检查：

mpow(P, 17)
#             [,1]       [,2]      [,3]      [,4]      [,5]      [,6]
# [1,] 0.002589361 0.02589806 0.1165645 0.3108818 0.2720684 0.2719978
# [2,] 0.002589415 0.02589722 0.1165599 0.3108747 0.2720749 0.2720039
# [3,] 0.002589738 0.02589714 0.1165539 0.3108615 0.2720788 0.2720189
# [4,] 0.002590797 0.02590083 0.1165520 0.3108412 0.2720638 0.2720515
# [5,] 0.002592925 0.02592074 0.1166035 0.3108739 0.2719451 0.2720638
# [6,] 0.002588814 0.02590459 0.1166029 0.3109419 0.2720166 0.2719451

只有前 4 位数字与您输入的 tol = 1e-4 相同。

浮点数最多有16位，所以我建议你使用tol = 1e-16进行可靠的收敛测试。