diff --git a/.jules/bolt.md b/.jules/bolt.md
new file mode 100644
index 0000000..79d1d86
--- /dev/null
+++ b/.jules/bolt.md
@@ -0,0 +1,3 @@
+## 2024-05-24 - Avoid Double Matrix Inversion
+**Learning:** Found a performance bottleneck where a covariance matrix is inverted, passed to another function, and then inverted again (`chol2inv(chol(A))`). This is not only an unnecessary $O(p^3)$ performance overhead, but also risks precision loss due to repeated floating point calculations.
+**Action:** When a matrix and its inverse are needed, consider whether passing the original un-inverted matrix is sufficient, or pass both. In `vuongtest.R`, we can just pass the original variance-covariance matrix (`tmpvc`) rather than inverting it in `calcAB` and inverting it back in `calcLambda`.
diff --git a/R/vuongtest.R b/R/vuongtest.R
index 57c5f6b..8368062 100644
--- a/R/vuongtest.R
+++ b/R/vuongtest.R
@@ -232,6 +232,10 @@ calcAB <- function(object, n, scfun, vc){
     if(nrow(tmpvc) == 1 & is.na(tmpvc[1,1])) stop("Please re-estimate the mirt model with SE=TRUE")
   }
   A <- chol2inv(chol(tmpvc))
+  ## Bolt performance optimization: return un-inverted matrix `tmpvc` as `Ainv`
+  ## to avoid the $O(p^3)$ performance overhead and precision loss of
+  ## double-inverting in `calcLambda`.
+  Ainv <- tmpvc
 
   ## Eq (2.2)
   if(!is.null(scfun)){
@@ -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,10 @@ 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))))
+  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?