diff --git a/DESCRIPTION b/DESCRIPTION
index c2bd031fa..117d5f0d1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -47,7 +47,6 @@ Imports:
   ape (>= 5.0),
   cli (>= 3.0),
   colorspace,
-  memoise,
   Rdpack (>= 0.7),
   shiny,
   shinyjs,
diff --git a/NAMESPACE b/NAMESPACE
index 3fcd911aa..4657fc93f 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -178,7 +178,6 @@ importFrom(graphics,plot)
 importFrom(graphics,points)
 importFrom(graphics,rect)
 importFrom(graphics,segments)
-importFrom(memoise,memoise)
 importFrom(parallel,makeCluster)
 importFrom(parallel,parCapply)
 importFrom(parallel,parLapply)
diff --git a/R/Information.R b/R/Information.R
index ad5be0774..82be9ac9b 100644
--- a/R/Information.R
+++ b/R/Information.R
@@ -169,39 +169,54 @@ MeilaMutualInformation <- function(split1, split2) {
 #' @references \insertAllCited{}
 #' 
 #' @encoding UTF-8
-#' @importFrom memoise memoise
 #' @export
-AllSplitPairings <- memoise(function(n) {
+AllSplitPairings <- local({
+  # Simple memoization cache
+  cache <- new.env(parent = emptyenv())
   
-  if (n < 4L) stop("No informative splits with < 4 taxa")
+  function(n) {
+    # Check cache first
+    key <- as.character(n)
+    if (exists(key, envir = cache)) {
+      return(get(key, envir = cache))
+    }
+    
+    # Compute result
+    if (n < 4L) stop("No informative splits with < 4 taxa")
+    
+    # smallHalves <- 1L + seq_len(ceiling(n / 2) - 2L)
+    dataRows <- 2L
   
-  # smallHalves <- 1L + seq_len(ceiling(n / 2) - 2L)
-  dataRows <- 2L
-
-  unevenPairs <- matrix(
-    # For i in 2:largestSmallSplit
-    #TODO: Make faster by not calculating bottom triangle
-    unlist(lapply(1L + seq_len(n - 3L), function(inA) {
-      # For j in 2:(n - 2)
-      nCa <- choose(n, inA)
-      outA <- n - inA
-      hA <- Entropy(c(inA, outA) / n)
-      unlist(lapply(1L + seq_len(n - 3L), function(inB) {
-        outB <- n - inB
-        hB <- Entropy(c(inB, outB) / n)
-        vapply(max(0, inA + inB - n):min(inA, inB), function(inAB) {
-          association <- c(inAB, inA - inAB, inB - inAB, n + inAB - inA - inB)
-          jointEntropies <- Entropy(association / n)
-          
-          c(#inA, inB, inAB, 
-            #npairs = NPartitionPairs(association), nis = choose(n, i),
-            nTotal = nCa * choose(inA, inAB) * choose(outA, inB - inAB),
-            VoI = jointEntropies + jointEntropies - hA - hB)
-        }, double(dataRows))
-      }))
-    })), dataRows, dimnames=list(c("nTotal", "VoI"), NULL))
-  
-  tapply(unevenPairs["nTotal", ], unevenPairs["VoI", ], sum)
+    unevenPairs <- matrix(
+      # For i in 2:largestSmallSplit
+      #TODO: Make faster by not calculating bottom triangle
+      unlist(lapply(1L + seq_len(n - 3L), function(inA) {
+        # For j in 2:(n - 2)
+        nCa <- choose(n, inA)
+        outA <- n - inA
+        hA <- Entropy(c(inA, outA) / n)
+        unlist(lapply(1L + seq_len(n - 3L), function(inB) {
+          outB <- n - inB
+          hB <- Entropy(c(inB, outB) / n)
+          vapply(max(0, inA + inB - n):min(inA, inB), function(inAB) {
+            association <- c(inAB, inA - inAB, inB - inAB, n + inAB - inA - inB)
+            jointEntropies <- Entropy(association / n)
+            
+            c(#inA, inB, inAB, 
+              #npairs = NPartitionPairs(association), nis = choose(n, i),
+              nTotal = nCa * choose(inA, inAB) * choose(outA, inB - inAB),
+              VoI = jointEntropies + jointEntropies - hA - hB)
+          }, double(dataRows))
+        }))
+      })), dataRows, dimnames=list(c("nTotal", "VoI"), NULL))
+    
+    result <- tapply(unevenPairs["nTotal", ], unevenPairs["VoI", ], sum)
+    
+    # Cache the result
+    assign(key, result, envir = cache)
+    
+    return(result)
+  }
 })
 
 #' Entropy of two splits