diff --git a/.jules/bolt.md b/.jules/bolt.md
index f658475..1656099 100644
--- a/.jules/bolt.md
+++ b/.jules/bolt.md
@@ -9,9 +9,15 @@
 ## 2024-05-25 - Avoid O(N^2) memory reallocation in R loops
 **Learning:** Using `do.call(cbind, ...)` to grow an N-row object across K submodels causes $O(NK^2)$ cumulative copying and $O(NK)$ peak storage.
 **Action:** Accumulate sums directly to keep $O(N)$ accumulator storage and $O(NK)$ total accumulation work.
+
 ## 2026-07-14 - Matrix Cross Product Optimization
 **Learning:** In R, matrix multiplication of the form `t(X) %*% Y` explicitly allocates memory for the transposed matrix. Using the optimized base function `crossprod(X, Y)` avoids this allocation.
 **Action:** Always replace `t(X) %*% Y` with `crossprod(X, Y)` for faster and more memory-efficient cross-product calculations.
+
 ## 2024-05-15 - [R Performance: ifelse Overhead]
 **Learning:** In R, ifelse evaluates both true and false branches entirely before subsetting, which is very inefficient for vector operations.
 **Action:** Optimize this by preallocating with res <- Y * 0 to preserve attributes and using vectorized subsetting like if any cond res subset <- ...
+
+## 2026-08-11 - Matrix indexing vs allocation and multiplication
+**Learning:** In R codebases, allocating a zero matrix, updating specific elements using `matrix(c(rows, cols), ncol=2)`, and multiplying before `rowSums` is extremely inefficient (O(N*K) space and time) compared to directly subsetting the values using the index matrix `[cbind(rows, cols)]`.
+**Action:** When extracting one element per row from a matrix based on a vector of column indices, never build an indicator matrix to multiply. Always use matrix subsetting `mat[cbind(1:nrow(mat), col_indices)]` which is magnitudes faster.
diff --git a/R/llcont.R b/R/llcont.R
index d8e496a..0bad6a0 100644
--- a/R/llcont.R
+++ b/R/llcont.R
@@ -363,12 +363,14 @@ llcont.nls <- function (x, ...) {
 llcont.polr <- function(x, ...) {
   m <- x$model
   y <- unclass(model.response(m))
-  wherey <- matrix(c(as.numeric(names(y)), y), ncol=2)
-  idx <- matrix(0, nrow=length(y), ncol=length(x$lev))
-  idx[wherey] <- 1
-
-  ## Bolt: replaced apply(..., 1, sum) with optimized rowSums() for performance
-  model.weights(m) * log(rowSums(idx * x$fitted.values))
+  wherey <- cbind(seq_along(y), as.numeric(y))
+
+  ## Bolt: replaced matrix creation and rowSums with direct matrix subsetting for performance
+  w <- model.weights(m)
+  res <- if (is.null(w)) log(x$fitted.values[wherey]) else w * log(x$fitted.values[wherey])
+  names(res) <- rownames(x$fitted.values)
+  if (is.null(names(res))) names(res) <- names(y)
+  res
 }
 
 ################################################################
diff --git a/benchmark_hurdle_ifelse.R b/benchmark_hurdle_ifelse.R
deleted file mode 100644
index eda1283..0000000
--- a/benchmark_hurdle_ifelse.R
+++ /dev/null
@@ -1,48 +0,0 @@
-# Reproducible benchmark harness
-# Install the optional benchmark dependency with:
-# install.packages("microbenchmark")
-library(microbenchmark)
-
-run_zeroPoisson_orig <- function(Z, parms, offsetz, weights, Y0, Y1) {
-  mu <- as.vector(exp(Z %*% parms + offsetz))
-  loglik0 <- -mu
-  Y0 * weights * loglik0 + ifelse(Y1, weights * log(1 - exp(loglik0)), 0)
-}
-
-run_zeroPoisson_opt <- function(Z, parms, offsetz, weights, Y0, Y1) {
-  mu <- as.vector(exp(Z %*% parms + offsetz))
-  loglik0 <- -mu
-  res_Y1 <- Y1 * 0
-  cond <- Y1; cond[is.na(cond)] <- FALSE
-  if (any(cond)) {
-    w_c <- if (length(weights) == 1) rep_len(weights, sum(cond)) else weights[cond]
-    res_Y1[cond] <- w_c * log(1 - exp(loglik0[cond]))
-  }
-  Y0 * weights * loglik0 + res_Y1
-}
-
-# Generate mostly zeros (so Y1 is mostly FALSE)
-set.seed(20260811)
-n <- 1000000
-Z <- matrix(rnorm(n*2), n, 2)
-parms <- c(0.5, -0.5)
-offsetz <- rep(0, n)
-Y <- rbinom(n, 1, 0.1)
-Y0 <- Y <= 0
-Y1 <- Y > 0
-weights <- 1
-
-bm <- microbenchmark(
-  original = run_zeroPoisson_orig(Z, parms, offsetz, weights, Y0, Y1),
-  optimized = run_zeroPoisson_opt(Z, parms, offsetz, weights, Y0, Y1),
-  times = 100,
-  control = list(warmup = 10)
-)
-print(bm)
-
-med_orig <- median(bm$time[bm$expr == "original"])
-med_opt <- median(bm$time[bm$expr == "optimized"])
-improvement <- (med_orig - med_opt) / med_orig
-
-cat(sprintf("Performance improvement: %.2f%%\n", improvement * 100))
-cat("Timing is descriptive; compare thresholds only in a controlled environment.\n")