R choldrop 删除并排名第一 Cholesky 因子更新

说明

给定矩阵 A 的 Cholesky 因子 R ， choldrop 查找 A[-k,-k] 的 Cholesky 因子，其中 k 是整数。 cholup 查找A + uu^T(更新)或A - uu^T(降级)的因子。

用法

choldrop(R,k)
cholup(R,u,up)

参数

R	矩阵的 Cholesky 因子，A 。
k	要删除的 A 的行和列。
u	定义排名一更新的向量。
up	如果TRUE计算更新，否则downdate。

细节

首先考虑choldrop。如果 R 是上三角，则 t(R[,-k])%*%R[,-k] == A[-k,-k] ，但 R[,-k] 在第一个 sub-diagonal 上有元素，从第 k 列开始。为了得到 A[-k,-k] 的三角 Cholesky 因子，我们可以从左侧应用一系列吉文斯旋转来消除 sub-diagonal 元素。例程就是这样做的。如果 R 是下三角因子，则需要从右侧进行吉文斯旋转来删除多余的元素。如果n 是R 的维度，则更新具有O(n^2) 计算成本。

cholup(假设 R 是上三角)根据 R^TR + uu^T = [u,R^T][u,R^T]^T = [u,R^T]Q^TQ[u,R^T]^T 的观察进行更新，因此我们可以构造 Q 以便 Q[u,R^T]^T=[0,R_1^T]^T ，其中 R_1 是修改后的因子。 Q 由一系列吉文斯旋转构造而成，以便将 u 的元素归零。降日期是类似的，只是必须使用双曲旋转来代替吉文斯旋转 - 有关详细信息，请参阅 Golub 和 van Loan(2013 年，第 6.5.4 节)。仅当 A - uu^T 为正定时，降期才有效。同样，计算成本为 O(n^2) 。

请注意，更新是面向矢量的，因此不易通过使用优化的 BLAS 来加速。更新被设置为相对缓存友好，因为在上三角情况下，连续的吉文斯旋转被存储以用于按列顺序应用，而不是在计算后立即按行应用。即便如此，上三角更新还是比下三角更新稍慢。

例子

  require(mgcv)
  set.seed(0)
  n <- 6
  A <- crossprod(matrix(runif(n*n),n,n))
  R0 <- chol(A)
  k <- 3
  Rd <- choldrop(R0,k)
  range(Rd-chol(A[-k,-k]))
  Rd;chol(A[-k,-k])
  
  ## same but using lower triangular factor A = LL'
  L <- t(R0)
  Ld <- choldrop(L,k)
  range(Ld-t(chol(A[-k,-k])))
  Ld;t(chol(A[-k,-k]))

  ## Rank one update example
  u <- runif(n)
  R <- cholup(R0,u,TRUE)
  Ru <- chol(A+u %*% t(u)) ## direct for comparison
  R;Ru
  range(R-Ru)

  ## Downdate - just going back from R to R0
  Rd <-  cholup(R,u,FALSE)
  R0;Rd
  range(R0-Rd)
  

作者

Simon N. Wood simon.wood@r-project.org

参考

Golub GH and CF Van Loan (2013) Matrix Computations (4th edition) Johns Hopkins