6

这个问题在参考是来自代码高尔夫挑战的观察。

提交的 R 解决方案是一个可行的解决方案，但我们中的一些人（也许只有我）似乎对为什么X=m需要进行初始重新分配感到目瞪口呆。

@Giuseppe 对代码进行了一些修改，因此我将为读者写一些评论。

function(m){
            X=m
            # Re-assign input m as X

            while(any(X-(X=X%*%m))) 0
            # Instead of doing the meat of the calculation in the code block after `while`
            # OP exploited its infinite looping properties to perform the
            # calculations within the condition check.
            # `-` here is an abuse of inequality check and relies on `any` to coerce
            # the numeric to logical. See `as.logical(.Machine$double.xmin)`
            # The code basically multiplies the matrix `X` with the starting matrix `m`            
            # Until the condition is met: X == X%*%m

            X
            # Return result
           }

据我所知。乘法X%*%m等价于X%*%X因为X是 的一个迭代自乘版本m。矩阵收敛后，将其他副本相乘m或X不改变其值。请参阅线性代数教科书或v(m)%*%v(m)%*%v(m)%*%v(m)%*%v(m)%*%m%*%m将上述函数定义为v。好玩吧？

那么问题来了，为什么@CodesInChaos 对这个想法的实现行不通呢？

function(m){while(any(m!=(m=m%*%m)))0 m}

这是由浮点精度问题引起的吗？或者这是由代码中的 a 函数引起的，例如不等式检查或 .Primitive("any")？我不相信这是as.logical由于 R 似乎强制错误小于.Machine$double.xmin0 造成的。

这是上面的演示。m我们只是循环并获取和之间的差异m%*%m。当我们尝试收敛随机矩阵时，该误差变为 0。它似乎收敛然后吹到 0/INF 最终取决于输入。

mat = matrix(c(7/10, 4/10, 3/10, 6/10), 2, 2, byrow = T)

m = mat
for (i in 1:25) {
  m = m%*%m
  cat("Mean Error:", mean(m-(m=m%*%m)), 
      "\n Float to Logical:", as.logical(m-(m=m%*%m)),
      "\n iter", i, "\n")
}

关于为什么这是一个浮点数学问题的一些额外想法

1）循环表明这可能不是any任何逻辑检查/转换步骤的问题，而是与浮点矩阵数学有关。

2）@user202729 在原始线程中的评论说这个问题在果冻中仍然存在，代码高尔夫语言更加相信这可能是一个浮点问题。

4

2 回答 2

5

不同的方法迭代不同的函数，都从种子值开始m。如果该固定点是稳定的并且种子在该固定点的吸引力范围内，则函数迭代只会收敛到给定的固定点。

在原始代码中，您正在迭代函数

f <- function(X) X %*% m

极限矩阵是一个稳定的不动点，假设（在 Code Gulf 问题中说明）存在明确定义的极限。由于函数定义取决于m，因此不动点是 的函数也就不足为奇了m。

另一方面，建议的变体使用m = m %*% m是通过迭代函数获得的

g <- function(X) X %*% X

请注意，所有幂等矩阵都是此函数的不动点，但显然它们不能都是稳定的不动点。显然，原来固定函数中的极限矩阵不是一个稳定的不动点g（即使它是一个不动点）。

要真正确定这一点，您需要深入了解函数迭代下的矩阵不动点理论，以说明为什么在这种情况下的不动点g是不稳定的。

于 2018-04-04T15:12:12.250 回答
-1

这确实是一个浮点数学问题。要查看它，请查看此函数的结果：

test2 <- function(m) {
  c <- 0
  res <- list()
  while (any(m!=(m=m%*%m))) {
    c <- c + 1
    res[[c]] <- m
  }
  print(c)
  res
}

要测试具有一定容差的相等性，您可以使用：

test3 <- function(m) {
  while (!isTRUE(all.equal(m, m <- m %*% m))) 0
  m
}
于 2018-04-04T14:40:57.597 回答