R中的奇怪错误

Question

由于某种原因，每当我运行我的函数时，我的代码中的 if 语句都会在 R 中产生一些奇怪的错误：

Error in if (P > 0) { : missing value where TRUE/FALSE needed

如果 B 是单位矩阵，我的函数应该返回对称的半正定矩阵 A 的最小特征值。我使用

测试了我的功能

set.seed(12345)
a <- crossprod(matrix (rnorm(40), 10 ,4)) / 10
b <- diag(1, 4, 4)

它确实产生了正确的答案（我使用 eigen 函数进行了检查。）

无论如何，我需要将 theta 设置为最小化某个函数的值。对于我使用的示例，该值是 (-Q - sqrt(under.sqrt)) / (2 * P)，但通常它取决于 P 的符号。为什么我的 if 语句会给我这样的错误？我已经被困在这里一段时间了。任何帮助，将不胜感激。

myFunction <- function(A, B, eps = 1e-6, itmax = 100, verbose = FALSE) {

  n <- nrow(A)
  m <- ncol(A)
  x <- rep (1, n)
  u <- A %*% x
  v <- B %*% x
  r <- t(u) %*% x
  s <- t(v) %*% x
  y <- r / s

  itel <- 1

  repeat {
    tmax <- 0

    for(j in 1:m) {
      P <- (u[j] * B[j, j]) - (v[j] * A[j, j])
      Q <- (r * B[j, j]) - (s * A[j, j])
      R <- (r * v[j]) - (s * u[j])

      under.sqrt <- Q^2 - 4 * P * R

      # Error right here
      if (P > 0) {
        theta <- (-Q + sqrt(under.sqrt)) / (2 * P)
      }

      else {
        theta <- (-Q - sqrt(under.sqrt)) / (2 * P)
      }


      x[j] <- x[j] + theta
      r <- r + 2 * theta * u[j] + A[j, j] * theta * theta
      s <- s + 2 * theta * v[j] + B[j, j] * theta * theta
      u <- u + theta * A[ ,j]
      v <- v + theta * B[, j]

      y <- r / s

      tmax <- max(tmax, abs(theta))
    }

      if ((tmax < eps) || (itel == itmax)) {
        break()
      }

      itel <- itel + 1  

  }

  return(y)
}

Answer 1

你的算法不稳定。使用options(error=recover)可以浏览出错的函数：

Browse[1]> P
[1] NaN
Browse[1]> theta
     [,1]
[1,]  Inf

在这里您可以看到P 是NaN，这是由theta 爆炸到无限值引起的。这会传递给 x、u、v 并最终传递给 P，并导致比较失败。

Answer 2

这是我的赌注：P 是 NA

> if(NA > 0){}
Error in if (NA > 0) { : missing value where TRUE/FALSE needed

P 渲染 NaN：

> if(NaN > 0){}
Error in if (NaN > 0) { : missing value where TRUE/FALSE needed

关于P是向量的猜测：

> if(1:2 > 0){}
NULL
Warning message:
In if (1:2 > 0) { :
  the condition has length > 1 and only the first element will be used

Answer 3

你确定P 是一个单一的值吗？即，不是长度> 1的向量？

Answer 4

试试ifelse:

theta <- ifelse(P > 0,
                (-Q + sqrt(under.sqrt)) / (2 * P),
                (-Q - sqrt(under.sqrt)) / (2 * P))

ifelse 可以处理 NA 比较：

ifelse(NA>0,print("x"),print("y"))
[1] NA

if(NA>0) print("x") else print("y")
Error in if (NA > 0) print("x") else print("y") : 
  missing value where TRUE/FALSE needed