From 1c59678cd28ef330844959e1781a462d274bc8f4 Mon Sep 17 00:00:00 2001
From: seonghobae <8172694+seonghobae@users.noreply.github.com>
Date: Sat, 27 Jun 2026 03:47:09 +0000
Subject: [PATCH] Remove redundant O(N^3) matrix inversions in vuongtest

Pre-compute and pass the non-inverted matrix to avoid re-calculating the inverse of the inverse.
---
 .jules/bolt.md |  3 +++
 R/vuongtest.R  | 16 +++++++++++-----
 2 files changed, 14 insertions(+), 5 deletions(-)
 create mode 100644 .jules/bolt.md

diff --git a/.jules/bolt.md b/.jules/bolt.md
new file mode 100644
index 0000000..d73341a
--- /dev/null
+++ b/.jules/bolt.md
@@ -0,0 +1,3 @@
+## 2024-05-24 - Redundant Matrix Inversions in Statistical Operations
+**Learning:** Inverting a matrix and then inverting it again is a redundant O(N^3) operation that causes a massive performance bottleneck and potential precision loss. In `nonnest2`, the code computed `A <- chol2inv(chol(tmpvc))` and then later computed `chol2inv(chol(A))` to get back the original matrix `tmpvc`.
+**Action:** When reviewing statistical/mathematical code, trace the lifecycle of expensive matrices to find opportunities to pass pre-computed forms rather than recalculating inverses or Cholesky decompositions.
diff --git a/R/vuongtest.R b/R/vuongtest.R
index 57c5f6b..c3d9cee 100644
--- a/R/vuongtest.R
+++ b/R/vuongtest.R
@@ -231,6 +231,10 @@ calcAB <- function(object, n, scfun, vc){
     ## in case mirt vcov was not estimated
     if(nrow(tmpvc) == 1 & is.na(tmpvc[1,1])) stop("Please re-estimate the mirt model with SE=TRUE")
   }
+
+  # Optimization: We also pre-compute Ainv = tmpvc to save O(N^3) time later
+  # in calcLambda by avoiding chol2inv(chol(A)).
+  Ainv <- tmpvc
   A <- chol2inv(chol(tmpvc))
 
   ## Eq (2.2)
@@ -253,7 +257,7 @@ calcAB <- function(object, n, scfun, vc){
   sc.cp <- crossprod(sc)/n
   B <- matrix(sc.cp, nrow(A), nrow(A))
 
-  list(A=A, B=B, sc=sc)
+  list(A=A, Ainv=Ainv, B=B, sc=sc)
 }
 
 ## a function to get the cross-product from Eq (2.7)
@@ -271,10 +275,12 @@ calcLambda <- function(object1, object2, n, score1, score2, vc1, vc2) {
   AB2 <- calcAB(object2, n, score2, vc2)
   Bc <- calcBcross(AB1$sc, AB2$sc, n)
 
-  W <- cbind(rbind(-AB1$B %*% chol2inv(chol(AB1$A)),
-                   t(Bc) %*% chol2inv(chol(AB1$A))),
-             rbind(-Bc %*% chol2inv(chol(AB2$A)),
-                   AB2$B %*% chol2inv(chol(AB2$A))))
+  # Optimization: AB1$Ainv and AB2$Ainv are already A^-1, so we use them directly
+  # instead of computing chol2inv(chol(A)).
+  W <- cbind(rbind(-AB1$B %*% AB1$Ainv,
+                   t(Bc) %*% AB1$Ainv),
+             rbind(-Bc %*% AB2$Ainv,
+                   AB2$B %*% AB2$Ainv))
 
   lamstar <- eigen(W, only.values=TRUE)$values
   ## Discard imaginary part, as it only occurs for tiny eigenvalues?