diff --git a/.jules/bolt.md b/.jules/bolt.md
new file mode 100644
index 0000000..2e5e980
--- /dev/null
+++ b/.jules/bolt.md
@@ -0,0 +1,3 @@
+## 2024-05-24 - Expensive operations in R
+**Learning:** Common performance bottlenecks in math-heavy R libraries often come from duplicated operations in algebraic functions, such as Cholesky decompositions or matrix inversions. Due to CRAN checks, any arbitrary test scripts we write to evaluate these improvements shouldn't be left in the root package directory.
+**Action:** When working on math-heavy code, always search for repeated expensive computations, cache their results, and verify the improvements. Make sure to delete any local test scripts prior to submission to avoid interfering with CRAN submission checks.
diff --git a/R/vuongtest.R b/R/vuongtest.R
index 57c5f6b..cc82f7f 100644
--- a/R/vuongtest.R
+++ b/R/vuongtest.R
@@ -271,10 +271,14 @@ 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 Optimization: Cache matrix inversions to prevent redundant $O(p^3)$ computation
+  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?