From a4a6a2496d554b370e418af5bcaadacb5924fd8f Mon Sep 17 00:00:00 2001 From: kss2k Date: Tue, 13 Jan 2026 10:08:18 +0100 Subject: [PATCH 1/7] Try some optimizations --- R/equations_lms.R | 80 +++++++++++++++++++++++++++++++++---------- src/equations_lms.cpp | 35 ++++++++++--------- 2 files changed, 80 insertions(+), 35 deletions(-) diff --git a/R/equations_lms.R b/R/equations_lms.R index 0858f70d..0729be0b 100644 --- a/R/equations_lms.R +++ b/R/equations_lms.R @@ -106,10 +106,12 @@ estepLmsGroup <- function(submodel, lastQuad = NULL, recalcQuad = FALSE, data.id <- data$data.split[[j]] colidx <- data$colidx[[j]] - pj <- p[offset:end] - wm <- colSums(data.id * pj) / sum(pj) - X <- data.id - matrix(wm, nrow=nrow(data.id), ncol=ncol(data.id), byrow=TRUE) - wcov <- t(X) %*% (X * pj) + pj <- p[offset:end] + pj_sum <- sum(pj) + wm <- drop(crossprod(pj, data.id)) / pj_sum + centered <- sweep(data.id, 2, wm, "-") + weighted <- sweep(centered, 1, sqrt(pj), "*") + wcov <- crossprod(weighted) wMeans[[i]][[j]] <- wm wCovs[[i]][[j]] <- wcov @@ -137,14 +139,32 @@ mstepLms <- function(theta, model, P, optim.method = "L-BFGS-B", epsilon = 1e-6, ...) { + getFilled <- local({ + cache <- new.env(parent = emptyenv()) + cache$theta <- NULL + cache$mod <- NULL + + function(theta) { + if (!is.null(cache$theta) && identical(cache$theta, theta)) { + cache$mod + } else { + cache$theta <- theta + cache$mod <- fillModel(model = model, theta = theta, method = "lms") + cache$mod + } + } + }) + gradient <- function(theta) { + modFilled <- getFilled(theta) gradientCompLogLikLms(theta = theta, model = model, P = P, sign = -1, - epsilon = epsilon) + epsilon = epsilon, modFilled = modFilled) } objective <- function(theta) { + modFilled <- getFilled(theta) compLogLikLms(theta = theta, model = model, P = P, sign = -1, - epsilon = epsilon) + epsilon = epsilon, modFilled = modFilled) } if (optimizer == "nlminb") { @@ -172,9 +192,10 @@ mstepLms <- function(theta, model, P, } -compLogLikLms <- function(theta, model, P, sign = -1, ...) { +compLogLikLms <- function(theta, model, P, sign = -1, modFilled = NULL, ...) { tryCatch({ - modFilled <- fillModel(model = model, theta = theta, method = "lms") + if (is.null(modFilled)) + modFilled <- fillModel(model = model, theta = theta, method = "lms") ll <- 0 for (g in seq_len(model$info$n.groups)) { @@ -211,24 +232,34 @@ compLogLikLmsGroup <- function(submodel, P, sign = -1, ...) { -gradientCompLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6) { +gradientCompLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6, + modFilled = NULL) { gradientAllLogLikLms(theta = theta, model = model, P = P, sign = sign, - epsilon = epsilon, FGRAD = gradLogLikLmsCpp, FOBJECTIVE = compLogLikLmsGroup) + epsilon = epsilon, FGRAD = gradLogLikLmsCpp, + FOBJECTIVE = compLogLikLmsGroup, modFilled = modFilled) } gradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6, - FGRAD, FOBJECTIVE) { + FGRAD, FOBJECTIVE, modFilled = NULL) { hasCovModel <- model$params$gradientStruct$hasCovModel - if (hasCovModel) gradient <- \(...) complicatedGradientAllLogLikLms(..., FOBJECTIVE = FOBJECTIVE) - else gradient <- \(...) simpleGradientAllLogLikLms(..., FGRAD = FGRAD) + if (hasCovModel) { + gradient <- \(...) complicatedGradientAllLogLikLms(..., + FOBJECTIVE = FOBJECTIVE, modFilled = modFilled, theta.base = theta) + } else { + gradient <- \(...) simpleGradientAllLogLikLms(..., + FGRAD = FGRAD, modFilled = modFilled) + } c(gradient(theta = theta, model = model, P = P, sign = sign, epsilon = epsilon)) } -complicatedGradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-4, FOBJECTIVE, sum = TRUE, ...) { +complicatedGradientAllLogLikLms <- function(theta, model, P, sign = -1, + epsilon = 1e-4, FOBJECTIVE, + sum = TRUE, modFilled = NULL, + theta.base = NULL, ...) { params <- model$params SELECT_THETA_LAB <- params$SELECT_THETA_LAB @@ -246,9 +277,17 @@ complicatedGradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon grad <- matrix(0, nrow = n, ncol = k, dimnames = list(NULL, names(theta))) - FOBJECTIVE_GROUP <- function(theta, g) { - modFilled <- fillModel(theta = theta, model = model, method = "lms") - FOBJECTIVE(submodel = modFilled$models[[g]], P = P$P_GROUPS[[g]], sign = sign, ...) + theta.base <- if (is.null(theta.base)) theta else theta.base + use_cache <- !is.null(modFilled) + + FOBJECTIVE_GROUP <- function(theta.eval, g) { + if (use_cache && identical(theta.eval, theta.base)) { + submodel <- modFilled$models[[g]] + } else { + modFilledEval <- fillModel(theta = theta.eval, model = model, method = "lms") + submodel <- modFilledEval$models[[g]] + } + FOBJECTIVE(submodel = submodel, P = P$P_GROUPS[[g]], sign = sign, ...) } for (g in seq_len(model$info$n.groups)) { @@ -276,10 +315,13 @@ complicatedGradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon } -simpleGradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6, FGRAD, N.FGRAD = 1L) { +simpleGradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6, + FGRAD, N.FGRAD = 1L, modFilled = NULL) { # simple gradient which should work if constraints are well-behaved functions # which can be derivated by Deriv::Deriv, and there is no covModel - modelR <- fillModel(model=model, theta=theta, method="lms") + modelR <- if (is.null(modFilled)) + fillModel(model = model, theta = theta, method = "lms") + else modFilled locations <- model$params$gradientStruct$locations Jacobian <- model$params$gradientStruct$Jacobian nlinDerivs <- model$params$gradientStruct$nlinDerivs diff --git a/src/equations_lms.cpp b/src/equations_lms.cpp index e12bb5f9..598459c2 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,33 @@ struct LMSModel { Psi = Rcpp::as(matrices["psi"]); d = Rcpp::as(matrices["thetaDelta"]); e = Rcpp::as(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 tmp = arma::solve(arma::trans(B), Psi, arma::solve_opts::fast); + const arma::mat noise = arma::solve(B, tmp, arma::solve_opts::fast); + const arma::mat varEta = Eta * Oi * Eta.t() + noise; + const arma::mat covXiEta = AOi * Eta.t(); const arma::mat vcovXiEta = arma::join_cols( arma::join_rows(varXi, covXiEta), @@ -185,6 +187,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; } From ee83aa427b474598fa592bc2a9ba0fbe1cf25811 Mon Sep 17 00:00:00 2001 From: kss2k Date: Tue, 13 Jan 2026 10:27:57 +0100 Subject: [PATCH 2/7] Add benchmark workflow --- .github/workflows/modsem-benchmarks.yml | 110 +++++++++++++++++ inst/benchmarks/README.md | 29 +++++ inst/benchmarks/benchmark_models.R | 150 ++++++++++++++++++++++++ 3 files changed, 289 insertions(+) create mode 100644 .github/workflows/modsem-benchmarks.yml create mode 100644 inst/benchmarks/README.md create mode 100644 inst/benchmarks/benchmark_models.R diff --git a/.github/workflows/modsem-benchmarks.yml b/.github/workflows/modsem-benchmarks.yml new file mode 100644 index 00000000..2ac4ee10 --- /dev/null +++ b/.github/workflows/modsem-benchmarks.yml @@ -0,0 +1,110 @@ +name: modsem benchmarks + +on: + pull_request: + workflow_dispatch: + +jobs: + benchmark: + 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 + 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 + path: ${{ env.RESULTS_DIR }} diff --git a/inst/benchmarks/README.md b/inst/benchmarks/README.md new file mode 100644 index 00000000..b7ac9963 --- /dev/null +++ b/inst/benchmarks/README.md @@ -0,0 +1,29 @@ +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=`: CSV destination (default `benchmark-results.csv`) +- `--iters=`: number of repetitions per model. You can also set + `MODSEM_BENCH_ITERS` in the environment. + +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.yml`, which +installs the CRAN release, the current `main` branch, and the PR commit into +separate libraries, runs this script for each, and uploads 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..6a4ce3b9 --- /dev/null +++ b/inst/benchmarks/benchmark_models.R @@ -0,0 +1,150 @@ +#!/usr/bin/env Rscript + +parse_args <- function(args) { + opts <- list( + output = "benchmark-results.csv", + iterations = as.integer(Sys.getenv("MODSEM_BENCH_ITERS", "1")) + ) + + 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 { + stop("Unknown argument: ", a) + } + } + + if (!is.finite(opts$iterations) || opts$iterations < 1L) + opts$iterations <- 1L + + opts +} + +args <- 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 + ) + ) +) + +all_timings <- lapply(benchmarks, function(entry) { + vals <- numeric(args$iterations) + for (i in seq_len(args$iterations)) { + cat(sprintf("Running %s (iteration %d/%d)...\n", + entry$name, i, args$iterations)) + vals[i] <- system.time({ + fit <- entry$run() + rm(fit) + })[["elapsed"]] + invisible(gc()) + } + data.frame( + model = entry$name, + iteration = seq_len(args$iterations), + elapsed = vals + ) +}) + +timings <- do.call(rbind, all_timings) + +dir.create(dirname(args$output), recursive = TRUE, showWarnings = FALSE) +utils::write.csv(timings, args$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) From 15b46ad982eb304c9ef65430161581b3d7e3e748 Mon Sep 17 00:00:00 2001 From: kss2k Date: Tue, 13 Jan 2026 10:31:32 +0100 Subject: [PATCH 3/7] Fix multiplication ordering --- src/equations_lms.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/equations_lms.cpp b/src/equations_lms.cpp index 598459c2..1b979d21 100644 --- a/src/equations_lms.cpp +++ b/src/equations_lms.cpp @@ -156,9 +156,11 @@ struct LMSModel { 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 tmp = arma::solve(arma::trans(B), Psi, arma::solve_opts::fast); - const arma::mat noise = arma::solve(B, tmp, arma::solve_opts::fast); - const arma::mat varEta = Eta * Oi * Eta.t() + noise; + const arma::mat BinvT = arma::solve(arma::trans(B), + arma::eye(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( From be5c115f28f5e086b98c0f92603d25a981b2b317 Mon Sep 17 00:00:00 2001 From: kss2k Date: Tue, 13 Jan 2026 10:37:33 +0100 Subject: [PATCH 4/7] Increase the number of benchmark iterations --- .github/workflows/modsem-benchmarks.yml | 1 + inst/benchmarks/benchmark_models.R | 3 ++- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/.github/workflows/modsem-benchmarks.yml b/.github/workflows/modsem-benchmarks.yml index 2ac4ee10..095b3d24 100644 --- a/.github/workflows/modsem-benchmarks.yml +++ b/.github/workflows/modsem-benchmarks.yml @@ -11,6 +11,7 @@ jobs: env: RSPM: https://packagemanager.posit.co/cran/__linux__/jammy/latest RESULTS_DIR: benchmark-results + MODSEM_BENCH_ITERS: 25 steps: - name: Checkout PR branch uses: actions/checkout@v4 diff --git a/inst/benchmarks/benchmark_models.R b/inst/benchmarks/benchmark_models.R index 6a4ce3b9..10291aaa 100644 --- a/inst/benchmarks/benchmark_models.R +++ b/inst/benchmarks/benchmark_models.R @@ -1,9 +1,10 @@ #!/usr/bin/env Rscript parse_args <- function(args) { + default_iters <- as.integer(Sys.getenv("MODSEM_BENCH_ITERS", "25")) opts <- list( output = "benchmark-results.csv", - iterations = as.integer(Sys.getenv("MODSEM_BENCH_ITERS", "1")) + iterations = default_iters ) for (a in args) { From 7b972429c71083aa15dcd52f42624444ee93625d Mon Sep 17 00:00:00 2001 From: kss2k Date: Tue, 13 Jan 2026 11:52:16 +0100 Subject: [PATCH 5/7] Split benchmark workflow --- ...nchmarks.yml => modsem-benchmarks-lms.yml} | 7 +- .github/workflows/modsem-benchmarks-qml.yml | 112 ++++++++++++++++++ inst/benchmarks/README.md | 20 ++-- inst/benchmarks/benchmark_models.R | 34 ++++-- 4 files changed, 153 insertions(+), 20 deletions(-) rename .github/workflows/{modsem-benchmarks.yml => modsem-benchmarks-lms.yml} (96%) create mode 100644 .github/workflows/modsem-benchmarks-qml.yml diff --git a/.github/workflows/modsem-benchmarks.yml b/.github/workflows/modsem-benchmarks-lms.yml similarity index 96% rename from .github/workflows/modsem-benchmarks.yml rename to .github/workflows/modsem-benchmarks-lms.yml index 095b3d24..08bfd53b 100644 --- a/.github/workflows/modsem-benchmarks.yml +++ b/.github/workflows/modsem-benchmarks-lms.yml @@ -1,17 +1,18 @@ -name: modsem benchmarks +name: modsem benchmarks (LMS) on: pull_request: workflow_dispatch: jobs: - benchmark: + 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 @@ -107,5 +108,5 @@ jobs: - name: Upload benchmark artifacts uses: actions/upload-artifact@v4 with: - name: modsem-benchmark-results + 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 index b7ac9963..cfd44818 100644 --- a/inst/benchmarks/README.md +++ b/inst/benchmarks/README.md @@ -12,18 +12,20 @@ Rscript inst/benchmarks/benchmark_models.R \ --iters=3 ``` -Arguments: - -- `--output=`: CSV destination (default `benchmark-results.csv`) -- `--iters=`: number of repetitions per model. You can also set - `MODSEM_BENCH_ITERS` in the environment. +- Arguments: + - `--output=`: CSV destination (default `benchmark-results.csv`) + - `--iters=`: 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.yml`, which -installs the CRAN release, the current `main` branch, and the PR commit into -separate libraries, runs this script for each, and uploads the resulting CSVs -and summaries as workflow artifacts. +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 index 10291aaa..7cd4b23f 100644 --- a/inst/benchmarks/benchmark_models.R +++ b/inst/benchmarks/benchmark_models.R @@ -2,9 +2,11 @@ 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 + iterations = default_iters, + models = default_models ) for (a in args) { @@ -12,6 +14,8 @@ parse_args <- function(args) { 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) } @@ -23,7 +27,7 @@ parse_args <- function(args) { opts } -args <- parse_args(commandArgs(trailingOnly = TRUE)) +cfg <- parse_args(commandArgs(trailingOnly = TRUE)) suppressPackageStartupMessages({ if (!requireNamespace("modsem", quietly = TRUE)) @@ -119,11 +123,25 @@ benchmarks <- list( ) ) +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(args$iterations) - for (i in seq_len(args$iterations)) { + vals <- numeric(cfg$iterations) + for (i in seq_len(cfg$iterations)) { cat(sprintf("Running %s (iteration %d/%d)...\n", - entry$name, i, args$iterations)) + entry$name, i, cfg$iterations)) vals[i] <- system.time({ fit <- entry$run() rm(fit) @@ -132,15 +150,15 @@ all_timings <- lapply(benchmarks, function(entry) { } data.frame( model = entry$name, - iteration = seq_len(args$iterations), + iteration = seq_len(cfg$iterations), elapsed = vals ) }) timings <- do.call(rbind, all_timings) -dir.create(dirname(args$output), recursive = TRUE, showWarnings = FALSE) -utils::write.csv(timings, args$output, row.names = FALSE) +dir.create(dirname(cfg$output), recursive = TRUE, showWarnings = FALSE) +utils::write.csv(timings, cfg$output, row.names = FALSE) print(timings) From 7751d65d224d755326ad76fad378ac8a27ff1ab3 Mon Sep 17 00:00:00 2001 From: kss2k Date: Tue, 13 Jan 2026 11:52:58 +0100 Subject: [PATCH 6/7] Revert R/equations_lms.R --- R/equations_lms.R | 80 +++++++++++------------------------------------ 1 file changed, 19 insertions(+), 61 deletions(-) diff --git a/R/equations_lms.R b/R/equations_lms.R index 0729be0b..0858f70d 100644 --- a/R/equations_lms.R +++ b/R/equations_lms.R @@ -106,12 +106,10 @@ estepLmsGroup <- function(submodel, lastQuad = NULL, recalcQuad = FALSE, data.id <- data$data.split[[j]] colidx <- data$colidx[[j]] - pj <- p[offset:end] - pj_sum <- sum(pj) - wm <- drop(crossprod(pj, data.id)) / pj_sum - centered <- sweep(data.id, 2, wm, "-") - weighted <- sweep(centered, 1, sqrt(pj), "*") - wcov <- crossprod(weighted) + pj <- p[offset:end] + wm <- colSums(data.id * pj) / sum(pj) + X <- data.id - matrix(wm, nrow=nrow(data.id), ncol=ncol(data.id), byrow=TRUE) + wcov <- t(X) %*% (X * pj) wMeans[[i]][[j]] <- wm wCovs[[i]][[j]] <- wcov @@ -139,32 +137,14 @@ mstepLms <- function(theta, model, P, optim.method = "L-BFGS-B", epsilon = 1e-6, ...) { - getFilled <- local({ - cache <- new.env(parent = emptyenv()) - cache$theta <- NULL - cache$mod <- NULL - - function(theta) { - if (!is.null(cache$theta) && identical(cache$theta, theta)) { - cache$mod - } else { - cache$theta <- theta - cache$mod <- fillModel(model = model, theta = theta, method = "lms") - cache$mod - } - } - }) - gradient <- function(theta) { - modFilled <- getFilled(theta) gradientCompLogLikLms(theta = theta, model = model, P = P, sign = -1, - epsilon = epsilon, modFilled = modFilled) + epsilon = epsilon) } objective <- function(theta) { - modFilled <- getFilled(theta) compLogLikLms(theta = theta, model = model, P = P, sign = -1, - epsilon = epsilon, modFilled = modFilled) + epsilon = epsilon) } if (optimizer == "nlminb") { @@ -192,10 +172,9 @@ mstepLms <- function(theta, model, P, } -compLogLikLms <- function(theta, model, P, sign = -1, modFilled = NULL, ...) { +compLogLikLms <- function(theta, model, P, sign = -1, ...) { tryCatch({ - if (is.null(modFilled)) - modFilled <- fillModel(model = model, theta = theta, method = "lms") + modFilled <- fillModel(model = model, theta = theta, method = "lms") ll <- 0 for (g in seq_len(model$info$n.groups)) { @@ -232,34 +211,24 @@ compLogLikLmsGroup <- function(submodel, P, sign = -1, ...) { -gradientCompLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6, - modFilled = NULL) { +gradientCompLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6) { gradientAllLogLikLms(theta = theta, model = model, P = P, sign = sign, - epsilon = epsilon, FGRAD = gradLogLikLmsCpp, - FOBJECTIVE = compLogLikLmsGroup, modFilled = modFilled) + epsilon = epsilon, FGRAD = gradLogLikLmsCpp, FOBJECTIVE = compLogLikLmsGroup) } gradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6, - FGRAD, FOBJECTIVE, modFilled = NULL) { + FGRAD, FOBJECTIVE) { hasCovModel <- model$params$gradientStruct$hasCovModel - if (hasCovModel) { - gradient <- \(...) complicatedGradientAllLogLikLms(..., - FOBJECTIVE = FOBJECTIVE, modFilled = modFilled, theta.base = theta) - } else { - gradient <- \(...) simpleGradientAllLogLikLms(..., - FGRAD = FGRAD, modFilled = modFilled) - } + if (hasCovModel) gradient <- \(...) complicatedGradientAllLogLikLms(..., FOBJECTIVE = FOBJECTIVE) + else gradient <- \(...) simpleGradientAllLogLikLms(..., FGRAD = FGRAD) c(gradient(theta = theta, model = model, P = P, sign = sign, epsilon = epsilon)) } -complicatedGradientAllLogLikLms <- function(theta, model, P, sign = -1, - epsilon = 1e-4, FOBJECTIVE, - sum = TRUE, modFilled = NULL, - theta.base = NULL, ...) { +complicatedGradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-4, FOBJECTIVE, sum = TRUE, ...) { params <- model$params SELECT_THETA_LAB <- params$SELECT_THETA_LAB @@ -277,17 +246,9 @@ complicatedGradientAllLogLikLms <- function(theta, model, P, sign = -1, grad <- matrix(0, nrow = n, ncol = k, dimnames = list(NULL, names(theta))) - theta.base <- if (is.null(theta.base)) theta else theta.base - use_cache <- !is.null(modFilled) - - FOBJECTIVE_GROUP <- function(theta.eval, g) { - if (use_cache && identical(theta.eval, theta.base)) { - submodel <- modFilled$models[[g]] - } else { - modFilledEval <- fillModel(theta = theta.eval, model = model, method = "lms") - submodel <- modFilledEval$models[[g]] - } - FOBJECTIVE(submodel = submodel, P = P$P_GROUPS[[g]], sign = sign, ...) + FOBJECTIVE_GROUP <- function(theta, g) { + modFilled <- fillModel(theta = theta, model = model, method = "lms") + FOBJECTIVE(submodel = modFilled$models[[g]], P = P$P_GROUPS[[g]], sign = sign, ...) } for (g in seq_len(model$info$n.groups)) { @@ -315,13 +276,10 @@ complicatedGradientAllLogLikLms <- function(theta, model, P, sign = -1, } -simpleGradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6, - FGRAD, N.FGRAD = 1L, modFilled = NULL) { +simpleGradientAllLogLikLms <- function(theta, model, P, sign = -1, epsilon = 1e-6, FGRAD, N.FGRAD = 1L) { # simple gradient which should work if constraints are well-behaved functions # which can be derivated by Deriv::Deriv, and there is no covModel - modelR <- if (is.null(modFilled)) - fillModel(model = model, theta = theta, method = "lms") - else modFilled + modelR <- fillModel(model=model, theta=theta, method="lms") locations <- model$params$gradientStruct$locations Jacobian <- model$params$gradientStruct$Jacobian nlinDerivs <- model$params$gradientStruct$nlinDerivs From 67c0ba9d57ce242cc6944248999fe1d1390b98e6 Mon Sep 17 00:00:00 2001 From: kss2k Date: Tue, 13 Jan 2026 12:15:45 +0100 Subject: [PATCH 7/7] Make som small optimizations in `src/equations_qml.cpp` --- src/equations_qml.cpp | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) 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(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() +