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