diff --git a/.github/workflows/modsem-benchmarks-lms.yml b/.github/workflows/modsem-benchmarks-lms.yml
new file mode 100644
index 00000000..08bfd53b
--- /dev/null
+++ b/.github/workflows/modsem-benchmarks-lms.yml
@@ -0,0 +1,112 @@
+name: modsem benchmarks (LMS)
+
+on:
+  pull_request:
+  workflow_dispatch:
+
+jobs:
+  benchmark_lms:
+    if: github.event_name != 'pull_request' || github.head_ref == 'lms-optimizations'
+    runs-on: ubuntu-latest
+    env:
+      RSPM: https://packagemanager.posit.co/cran/__linux__/jammy/latest
+      RESULTS_DIR: benchmark-results
+      MODSEM_BENCH_ITERS: 25
+      MODSEM_BENCH_MODELS: oneInt_lms,TPB_lms
+    steps:
+      - name: Checkout PR branch
+        uses: actions/checkout@v4
+        with:
+          fetch-depth: 0
+
+      - name: Checkout main branch
+        uses: actions/checkout@v4
+        with:
+          ref: main
+          path: main-branch
+          fetch-depth: 0
+
+      - name: Setup R
+        uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - name: Install remotes helper
+        run: Rscript -e 'install.packages("remotes", repos = Sys.getenv("RSPM"))'
+
+      - name: Install CRAN modsem
+        env:
+          TARGET_LIB: ${{ runner.temp }}/cran-lib
+        run: |
+          mkdir -p "$TARGET_LIB"
+          Rscript -e 'install.packages("modsem", lib=Sys.getenv("TARGET_LIB"), repos = Sys.getenv("RSPM"))'
+
+      - name: Benchmark CRAN modsem
+        env:
+          R_LIBS_USER: ${{ runner.temp }}/cran-lib
+        run: |
+          mkdir -p "$RESULTS_DIR"
+          Rscript inst/benchmarks/benchmark_models.R --output="$RESULTS_DIR/cran.csv"
+
+      - name: Install main modsem
+        env:
+          TARGET_LIB: ${{ runner.temp }}/main-lib
+        run: |
+          mkdir -p "$TARGET_LIB"
+          Rscript -e 'remotes::install_local("main-branch", lib=Sys.getenv("TARGET_LIB"), dependencies = TRUE, upgrade = "never")'
+
+      - name: Benchmark main modsem
+        env:
+          R_LIBS_USER: ${{ runner.temp }}/main-lib
+        run: |
+          Rscript inst/benchmarks/benchmark_models.R --output="$RESULTS_DIR/main.csv"
+
+      - name: Install PR modsem
+        env:
+          TARGET_LIB: ${{ runner.temp }}/pr-lib
+        run: |
+          mkdir -p "$TARGET_LIB"
+          Rscript -e 'remotes::install_local(".", lib=Sys.getenv("TARGET_LIB"), dependencies = TRUE, upgrade = "never")'
+
+      - name: Benchmark PR modsem
+        env:
+          R_LIBS_USER: ${{ runner.temp }}/pr-lib
+        run: |
+          Rscript inst/benchmarks/benchmark_models.R --output="$RESULTS_DIR/pr.csv"
+
+      - name: Summarize benchmark results
+        run: |
+          Rscript - <<'RS'
+          files <- list(
+            cran = file.path(Sys.getenv("RESULTS_DIR"), "cran.csv"),
+            main = file.path(Sys.getenv("RESULTS_DIR"), "main.csv"),
+            pr   = file.path(Sys.getenv("RESULTS_DIR"), "pr.csv")
+          )
+          stopifnot(all(file.exists(unlist(files))))
+          read_bench <- function(path, label) {
+            df <- utils::read.csv(path)
+            df$version <- label
+            df
+          }
+          df_list <- Map(read_bench, files, names(files))
+          df <- do.call(rbind, df_list)
+          summary <- aggregate(elapsed ~ model + version, df, mean)
+          utils::write.csv(df, file.path(Sys.getenv("RESULTS_DIR"), "raw_timings.csv"), row.names = FALSE)
+          utils::write.csv(summary, file.path(Sys.getenv("RESULTS_DIR"), "summary.csv"), row.names = FALSE)
+          print(summary)
+          wide <- reshape(summary, idvar = "model", timevar = "version", direction = "wide")
+          required <- c("elapsed.cran", "elapsed.main", "elapsed.pr")
+          missing <- setdiff(required, names(wide))
+          if (!length(missing)) {
+            wide$ratio_pr_vs_cran <- wide$elapsed.cran / wide$elapsed.pr
+            wide$ratio_pr_vs_main <- wide$elapsed.main / wide$elapsed.pr
+          }
+          cat("\nWide summary (seconds):\n")
+          print(wide)
+          RS
+
+      - name: Upload benchmark artifacts
+        uses: actions/upload-artifact@v4
+        with:
+          name: modsem-benchmark-results-lms
+          path: ${{ env.RESULTS_DIR }}
diff --git a/.github/workflows/modsem-benchmarks-qml.yml b/.github/workflows/modsem-benchmarks-qml.yml
new file mode 100644
index 00000000..4505755e
--- /dev/null
+++ b/.github/workflows/modsem-benchmarks-qml.yml
@@ -0,0 +1,112 @@
+name: modsem benchmarks (QML)
+
+on:
+  pull_request:
+  workflow_dispatch:
+
+jobs:
+  benchmark_qml:
+    if: github.event_name != 'pull_request' || github.head_ref == 'lms-optimizations'
+    runs-on: ubuntu-latest
+    env:
+      RSPM: https://packagemanager.posit.co/cran/__linux__/jammy/latest
+      RESULTS_DIR: benchmark-results
+      MODSEM_BENCH_ITERS: 25
+      MODSEM_BENCH_MODELS: oneInt_qml,TPB_qml
+    steps:
+      - name: Checkout PR branch
+        uses: actions/checkout@v4
+        with:
+          fetch-depth: 0
+
+      - name: Checkout main branch
+        uses: actions/checkout@v4
+        with:
+          ref: main
+          path: main-branch
+          fetch-depth: 0
+
+      - name: Setup R
+        uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - name: Install remotes helper
+        run: Rscript -e 'install.packages("remotes", repos = Sys.getenv("RSPM"))'
+
+      - name: Install CRAN modsem
+        env:
+          TARGET_LIB: ${{ runner.temp }}/cran-lib
+        run: |
+          mkdir -p "$TARGET_LIB"
+          Rscript -e 'install.packages("modsem", lib=Sys.getenv("TARGET_LIB"), repos = Sys.getenv("RSPM"))'
+
+      - name: Benchmark CRAN modsem
+        env:
+          R_LIBS_USER: ${{ runner.temp }}/cran-lib
+        run: |
+          mkdir -p "$RESULTS_DIR"
+          Rscript inst/benchmarks/benchmark_models.R --output="$RESULTS_DIR/cran.csv"
+
+      - name: Install main modsem
+        env:
+          TARGET_LIB: ${{ runner.temp }}/main-lib
+        run: |
+          mkdir -p "$TARGET_LIB"
+          Rscript -e 'remotes::install_local("main-branch", lib=Sys.getenv("TARGET_LIB"), dependencies = TRUE, upgrade = "never")'
+
+      - name: Benchmark main modsem
+        env:
+          R_LIBS_USER: ${{ runner.temp }}/main-lib
+        run: |
+          Rscript inst/benchmarks/benchmark_models.R --output="$RESULTS_DIR/main.csv"
+
+      - name: Install PR modsem
+        env:
+          TARGET_LIB: ${{ runner.temp }}/pr-lib
+        run: |
+          mkdir -p "$TARGET_LIB"
+          Rscript -e 'remotes::install_local(".", lib=Sys.getenv("TARGET_LIB"), dependencies = TRUE, upgrade = "never")'
+
+      - name: Benchmark PR modsem
+        env:
+          R_LIBS_USER: ${{ runner.temp }}/pr-lib
+        run: |
+          Rscript inst/benchmarks/benchmark_models.R --output="$RESULTS_DIR/pr.csv"
+
+      - name: Summarize benchmark results
+        run: |
+          Rscript - <<'RS'
+          files <- list(
+            cran = file.path(Sys.getenv("RESULTS_DIR"), "cran.csv"),
+            main = file.path(Sys.getenv("RESULTS_DIR"), "main.csv"),
+            pr   = file.path(Sys.getenv("RESULTS_DIR"), "pr.csv")
+          )
+          stopifnot(all(file.exists(unlist(files))))
+          read_bench <- function(path, label) {
+            df <- utils::read.csv(path)
+            df$version <- label
+            df
+          }
+          df_list <- Map(read_bench, files, names(files))
+          df <- do.call(rbind, df_list)
+          summary <- aggregate(elapsed ~ model + version, df, mean)
+          utils::write.csv(df, file.path(Sys.getenv("RESULTS_DIR"), "raw_timings.csv"), row.names = FALSE)
+          utils::write.csv(summary, file.path(Sys.getenv("RESULTS_DIR"), "summary.csv"), row.names = FALSE)
+          print(summary)
+          wide <- reshape(summary, idvar = "model", timevar = "version", direction = "wide")
+          required <- c("elapsed.cran", "elapsed.main", "elapsed.pr")
+          missing <- setdiff(required, names(wide))
+          if (!length(missing)) {
+            wide$ratio_pr_vs_cran <- wide$elapsed.cran / wide$elapsed.pr
+            wide$ratio_pr_vs_main <- wide$elapsed.main / wide$elapsed.pr
+          }
+          cat("\nWide summary (seconds):\n")
+          print(wide)
+          RS
+
+      - name: Upload benchmark artifacts
+        uses: actions/upload-artifact@v4
+        with:
+          name: modsem-benchmark-results-qml
+          path: ${{ env.RESULTS_DIR }}
diff --git a/inst/benchmarks/README.md b/inst/benchmarks/README.md
new file mode 100644
index 00000000..cfd44818
--- /dev/null
+++ b/inst/benchmarks/README.md
@@ -0,0 +1,31 @@
+Performance benchmarks
+======================
+
+`inst/benchmarks/benchmark_models.R` fits four representative models:
+the `oneInt` and `TPB` examples under both LMS and QML estimation. Elapsed
+times are written to a CSV so different builds (CRAN, `main`, PR) can be
+compared. Run it with any installed version of `modsem`:
+
+```bash
+Rscript inst/benchmarks/benchmark_models.R \
+  --output=benchmark-results/pr.csv \
+  --iters=3
+```
+
+- Arguments:
+  - `--output=<file>`: CSV destination (default `benchmark-results.csv`)
+  - `--iters=<n>`: number of repetitions per model (default 25, override with
+    the flag or `MODSEM_BENCH_ITERS`)
+  - `--models=name1,name2`: restrict the run to specific targets (same names
+    as printed in the output). You can also set `MODSEM_BENCH_MODELS`.
+
+The script prints per-run timing as well as summary statistics, and the CSV
+contains one row per model/iteration/method. It respects `R_LIBS_USER`, so you
+can point it at different installed versions (CRAN release, `main`, PR build,
+…) to compare performance.
+
+GitHub PRs automatically run `.github/workflows/modsem-benchmarks-lms.yml` and
+`.github/workflows/modsem-benchmarks-qml.yml`, which install the CRAN release,
+the current `main` branch, and the PR commit into separate libraries, run this
+script for each subset of models, and upload the resulting CSVs and summaries
+as workflow artifacts.
diff --git a/inst/benchmarks/benchmark_models.R b/inst/benchmarks/benchmark_models.R
new file mode 100644
index 00000000..7cd4b23f
--- /dev/null
+++ b/inst/benchmarks/benchmark_models.R
@@ -0,0 +1,169 @@
+#!/usr/bin/env Rscript
+
+parse_args <- function(args) {
+  default_iters <- as.integer(Sys.getenv("MODSEM_BENCH_ITERS", "25"))
+  default_models <- Sys.getenv("MODSEM_BENCH_MODELS", "")
+  opts <- list(
+    output = "benchmark-results.csv",
+    iterations = default_iters,
+    models = default_models
+  )
+
+  for (a in args) {
+    if (grepl("^--output=", a)) {
+      opts$output <- sub("^--output=", "", a)
+    } else if (grepl("^--iters=", a)) {
+      opts$iterations <- suppressWarnings(as.integer(sub("^--iters=", "", a)))
+    } else if (grepl("^--models=", a)) {
+      opts$models <- sub("^--models=", "", a)
+    } else {
+      stop("Unknown argument: ", a)
+    }
+  }
+
+  if (!is.finite(opts$iterations) || opts$iterations < 1L)
+    opts$iterations <- 1L
+
+  opts
+}
+
+cfg <- parse_args(commandArgs(trailingOnly = TRUE))
+
+suppressPackageStartupMessages({
+  if (!requireNamespace("modsem", quietly = TRUE))
+    stop("Package 'modsem' is not installed.")
+  library(modsem)
+})
+
+data("oneInt", package = "modsem")
+data("TPB", package = "modsem")
+
+oneInt_model <- "
+# Outer Model
+  X =~ x1 + lx2 * x2 + lx3 * x3
+  Z =~ z1 + z2 + z3
+  Y =~ y1 + y2 + y3
+
+# Inner Model
+  Y ~ X + Z
+
+# Constraints
+  x1 ~~ vx1*x1
+  x2 ~~ vx2*x2
+  x3 ~~ vx3*x3
+
+  proj_vx1 := 0.8
+  proj_vx2 := (lx2 ^ 2) * 0.8
+  proj_vx3 := (lx3 ^ 2) * 0.8
+
+  vx1 == 1 - proj_vx1
+  vx2 == 1 - proj_vx2
+  vx3 == 1 - proj_vx3
+"
+
+tpb_model <- "
+# Outer Model (Based on Hagger et al., 2007)
+  ATT =~ att1 + att2 + att3 + att4 + att5
+  SN  =~ sn1 + sn2
+  PBC =~ pbc1 + pbc2 + pbc3
+  INT =~ int1 + int2 + int3
+  BEH =~ b1 + b2
+
+# Inner Model (Based on Steinmetz et al., 2011)
+  INT ~ ATT + SN + PBC
+  BEH ~ INT + PBC
+  BEH ~ PBC:INT
+"
+
+benchmarks <- list(
+  list(
+    name = "oneInt_lms",
+    run = function() modsem(
+      model.syntax     = oneInt_model,
+      data             = oneInt,
+      method           = "lms",
+      standardize.data = TRUE,
+      mean.observed    = FALSE,
+      calc.se          = FALSE,
+      verbose          = FALSE
+    )
+  ),
+  list(
+    name = "oneInt_qml",
+    run = function() modsem(
+      model.syntax     = oneInt_model,
+      data             = oneInt,
+      method           = "qml",
+      standardize.data = TRUE,
+      mean.observed    = FALSE,
+      calc.se          = FALSE,
+      verbose          = FALSE
+    )
+  ),
+  list(
+    name = "TPB_lms",
+    run = function() modsem(
+      model.syntax = tpb_model,
+      data         = TPB,
+      method       = "lms",
+      calc.se      = FALSE,
+      verbose      = FALSE,
+      nodes        = 32
+    )
+  ),
+  list(
+    name = "TPB_qml",
+    run = function() modsem(
+      model.syntax = tpb_model,
+      data         = TPB,
+      method       = "qml",
+      calc.se      = FALSE,
+      verbose      = FALSE
+    )
+  )
+)
+
+if (nzchar(cfg$models)) {
+  select <- trimws(strsplit(cfg$models, ",", fixed = TRUE)[[1]])
+  select <- select[nzchar(select)]
+  if (length(select)) {
+    avail <- vapply(benchmarks, `[[`, character(1), "name")
+    missing <- setdiff(select, avail)
+    if (length(missing))
+      stop("Unknown benchmark(s): ", paste(missing, collapse = ", "))
+    benchmarks <- Filter(function(b) b$name %in% select, benchmarks)
+  }
+}
+
+stopifnot(length(benchmarks) > 0)
+
+all_timings <- lapply(benchmarks, function(entry) {
+  vals <- numeric(cfg$iterations)
+  for (i in seq_len(cfg$iterations)) {
+    cat(sprintf("Running %s (iteration %d/%d)...\n",
+                entry$name, i, cfg$iterations))
+    vals[i] <- system.time({
+      fit <- entry$run()
+      rm(fit)
+    })[["elapsed"]]
+    invisible(gc())
+  }
+  data.frame(
+    model = entry$name,
+    iteration = seq_len(cfg$iterations),
+    elapsed = vals
+  )
+})
+
+timings <- do.call(rbind, all_timings)
+
+dir.create(dirname(cfg$output), recursive = TRUE, showWarnings = FALSE)
+utils::write.csv(timings, cfg$output, row.names = FALSE)
+
+print(timings)
+
+summary_df <- aggregate(elapsed ~ model, data = timings,
+  FUN = function(x) c(mean = mean(x), median = stats::median(x), sd = stats::sd(x)))
+
+cat("\nSummary (seconds):\n")
+print(summary_df)
diff --git a/src/equations_lms.cpp b/src/equations_lms.cpp
index e12bb5f9..1b979d21 100644
--- a/src/equations_lms.cpp
+++ b/src/equations_lms.cpp
@@ -104,7 +104,7 @@ inline arma::vec make_zvec(unsigned k, unsigned numXis, const arma::vec& z) {
 
 struct LMSModel {
   arma::mat A, Oxx, Oex, Ie, lY, lX, tY, tX, Gx, Ge,
-    a, beta0, Psi, d, e;
+    a, beta0, Psi, d, e, Oi;
   unsigned  k       = 0;
   unsigned  numXis  = 0;
 
@@ -133,31 +133,35 @@ struct LMSModel {
     Psi    = Rcpp::as<arma::mat>(matrices["psi"]);
     d      = Rcpp::as<arma::mat>(matrices["thetaDelta"]);
     e      = Rcpp::as<arma::mat>(matrices["thetaEpsilon"]);
+    Oi     = make_Oi(k, numXis);
   }
 
   arma::vec mu(const arma::vec& z) const {
     const arma::vec zVec = make_zvec(k, numXis, z);
-    const arma::mat kronZ = arma::kron(Ie, beta0 + A * zVec);
-    const arma::mat Binv = arma::inv(Ie - Ge - kronZ.t() * Oex);
-
-    const arma::vec muXi  = beta0 + A * zVec;
-    const arma::vec muEta = Binv * (a +
-        Gx * (beta0 + A * zVec) +
-        kronZ.t() * Oxx * (beta0 + A * zVec));
+    const arma::vec muXi = beta0 + A * zVec;
+    const arma::mat kronZ = arma::kron(Ie, muXi);
+    const arma::mat B = Ie - Ge - kronZ.t() * Oex;
+    const arma::vec rhs = a + Gx * muXi + kronZ.t() * Oxx * muXi;
+    const arma::vec muEta = arma::solve(B, rhs, arma::solve_opts::fast);
 
     return tX + lX * arma::join_cols(muXi, muEta);
   }
 
   arma::mat Sigma(const arma::vec& z) const {
     const arma::vec zVec  = make_zvec(k, numXis, z);
-    const arma::mat kronZ = arma::kron(Ie, beta0 + A * zVec);
-    const arma::mat Binv = arma::inv(Ie - Ge - kronZ.t() * Oex);
-
-    const arma::mat Oi = make_Oi(k, numXis);
-    const arma::mat Eta = Binv * (Gx * A + kronZ.t() * Oxx * A);
-    const arma::mat varXi    = A * Oi * A.t();
-    const arma::mat varEta   = Eta * Oi * Eta.t() + Binv * Psi * Binv.t();
-    const arma::mat covXiEta = A * Oi * Eta.t();
+    const arma::vec muXi  = beta0 + A * zVec;
+    const arma::mat kronZ = arma::kron(Ie, muXi);
+    const arma::mat B = Ie - Ge - kronZ.t() * Oex;
+    const arma::mat rhs    = Gx * A + kronZ.t() * Oxx * A;
+    const arma::mat Eta    = arma::solve(B, rhs, arma::solve_opts::fast);
+    const arma::mat AOi    = A * Oi;
+    const arma::mat varXi  = AOi * A.t();
+    const arma::mat BinvT  = arma::solve(arma::trans(B),
+                                         arma::eye<arma::mat>(B.n_rows, B.n_cols),
+                                         arma::solve_opts::fast);
+    const arma::mat resvarEta = arma::solve(B, Psi * BinvT, arma::solve_opts::fast);
+    const arma::mat varEta = Eta * Oi * Eta.t() + resvarEta;
+    const arma::mat covXiEta = AOi * Eta.t();
 
     const arma::mat vcovXiEta = arma::join_cols(
         arma::join_rows(varXi, covXiEta),
@@ -185,6 +189,7 @@ struct LMSModel {
     c.Psi   = arma::mat(Psi);
     c.d     = arma::mat(d);
     c.e     = arma::mat(e);
+    c.Oi    = arma::mat(Oi);
 
     return c;
   }
diff --git a/src/equations_qml.cpp b/src/equations_qml.cpp
index 6de5d698..284d14f3 100644
--- a/src/equations_qml.cpp
+++ b/src/equations_qml.cpp
@@ -38,13 +38,11 @@ arma::mat muQmlCpp(Rcpp::List m, int t, int ncores = 1) {
 
       const arma::mat kronXi_t = kronXi.submat(firstRow, 0, lastRow, lastColKOxx);
       const arma::mat Binv_t   = Binv  .submat(firstRow, 0, lastRow, lastColBinv);
+      const arma::vec muXi     = beta0 + L1 * X.row(i).t();
+      const arma::vec kronTerm = kronXi_t * omegaXiXi * muXi;
+      const arma::vec rhs = trOmegaSigma + alpha + gammaXi * muXi + kronTerm;
 
-      Ey.row(i) = ( Binv_t *
-        ( trOmegaSigma + alpha
-          + gammaXi * (beta0 + L1 * X.row(i).t())
-          + kronXi_t * omegaXiXi * (beta0 + L1 * X.row(i).t()) )
-        + L2 * U.row(i).t()
-      ).t();
+      Ey.row(i) = ( Binv_t * rhs + L2 * U.row(i).t() ).t();
     }
 
   } else {
@@ -55,13 +53,11 @@ arma::mat muQmlCpp(Rcpp::List m, int t, int ncores = 1) {
     for (int i = 0; i < t; ++i) {
       const int firstRow = i * numEta, lastRow = (i + 1) * numEta - 1;
       const arma::mat kronXi_t = kronXi.submat(firstRow, 0, lastRow, lastColKOxx);
+      const arma::vec muXi     = beta0 + L1 * X.row(i).t();
+      const arma::vec kronTerm = kronXi_t * omegaXiXi * muXi;
+      const arma::vec rhs = trOmegaSigma + alpha + gammaXi * muXi + kronTerm;
 
-      Ey.row(i) = ( Binv *
-        ( trOmegaSigma + alpha
-          + gammaXi * (beta0 + L1 * X.row(i).t())
-          + kronXi_t * omegaXiXi * (beta0 + L1 * X.row(i).t()) )
-        + L2 * U.row(i).t()
-      ).t();
+      Ey.row(i) = ( Binv * rhs + L2 * U.row(i).t() ).t();
     }
   }
 
@@ -94,6 +90,7 @@ arma::mat sigmaQmlCpp(Rcpp::List m, int t, int ncores = 1) {
             lastColKOxx  = numXi * numEta - 1;
 
   const arma::mat omegaXiXi2T = omegaXiXi + transposeOmega(omegaXiXi, numEta); // omega + omega'
+  const arma::mat kronXiOmega = kronXi * omegaXiXi2T;
 
   if (Binv.n_rows > static_cast<unsigned>(numEta)) {
     #ifdef _OPENMP
@@ -102,12 +99,12 @@ arma::mat sigmaQmlCpp(Rcpp::List m, int t, int ncores = 1) {
 
     for (int i = 0; i < t; ++i) {
       const int firstRow = i * numEta, lastRow = (i + 1) * numEta - 1;
-      const arma::mat kronXi_t = kronXi.submat(firstRow, 0, lastRow, lastColKOxx);
+      const arma::mat kronOmega_t = kronXiOmega.submat(firstRow, 0, lastRow, lastColSigma);
       const arma::mat Binv_t   = Binv  .submat(firstRow, 0, lastRow, lastColSigma);
       const arma::mat Sigma2   = Binv_t * psi * Binv_t.t() + Sigma2Theta;
 
       const arma::mat BinvGammaXi2Omega =
-        Binv_t * gammaXi + kronXi_t * omegaXiXi2T;
+        Binv_t * gammaXi + kronOmega_t;
 
       sigmaE.submat(firstRow, 0, lastRow, lastColSigma) =
         BinvGammaXi2Omega * Sigma1 * BinvGammaXi2Omega.t() +
@@ -124,9 +121,9 @@ arma::mat sigmaQmlCpp(Rcpp::List m, int t, int ncores = 1) {
 
     for (int i = 0; i < t; ++i) {
       const int firstRow = i * numEta, lastRow = (i + 1) * numEta - 1;
-      const arma::mat kronXi_t = kronXi.submat(firstRow, 0, lastRow, lastColKOxx);
+      const arma::mat kronOmega_t = kronXiOmega.submat(firstRow, 0, lastRow, lastColSigma);
       const arma::mat BinvGammaXi2Omega =
-        Binv * gammaXi + kronXi_t * omegaXiXi2T;
+        Binv * gammaXi + kronOmega_t;
 
       sigmaE.submat(firstRow, 0, lastRow, lastColSigma) =
         BinvGammaXi2Omega * Sigma1 * BinvGammaXi2Omega.t() +