diff --git a/.jules/bolt.md b/.jules/bolt.md
new file mode 100644
index 0000000..6e74844
--- /dev/null
+++ b/.jules/bolt.md
@@ -0,0 +1,3 @@
+## 2024-05-18 - Caching Matrix Inversions in vuongtest
+**Learning:** In statistical R packages involving likelihood calculations, matrix inversion is a common performance bottleneck. Repeating the exact same inversion inside matrix bindings multiplies this cost unnecessarily.
+**Action:** Always inspect matrix construction blocks for duplicated expensive operations and extract them into cached local variables before assembling the final block matrix.
diff --git a/R/vuongtest.R b/R/vuongtest.R
index 57c5f6b..3864774 100644
--- a/R/vuongtest.R
+++ b/R/vuongtest.R
@@ -271,10 +271,13 @@ 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))))
+  # Bolt: cache matrix inversions to avoid recalculating the O(n^3) inverse twice for each model
+  invA1 <- chol2inv(chol(AB1$A))
+  invA2 <- chol2inv(chol(AB2$A))
+  W <- cbind(rbind(-AB1$B %*% invA1,
+                   t(Bc) %*% invA1),
+             rbind(-Bc %*% invA2,
+                   AB2$B %*% invA2))
 
   lamstar <- eigen(W, only.values=TRUE)$values
   ## Discard imaginary part, as it only occurs for tiny eigenvalues?