From e117f7a563ec9bee4f5bcd780099a8e5c140f1ac Mon Sep 17 00:00:00 2001 From: opencode agent Date: Fri, 17 Jul 2026 21:40:55 +0200 Subject: [PATCH] Reducing the number of times the gradients and fit are computed to improve the speed --- include/lesstimate/glmnet_class.h | 104 +++++++++++++++++++----------- include/lesstimate/ista_class.h | 32 ++++----- include/lesstimate/model.h | 26 ++++---- 3 files changed, 98 insertions(+), 64 deletions(-) diff --git a/include/lesstimate/glmnet_class.h b/include/lesstimate/glmnet_class.h index fda8614..107a306 100644 --- a/include/lesstimate/glmnet_class.h +++ b/include/lesstimate/glmnet_class.h @@ -238,7 +238,15 @@ namespace lessSEM const double sigma, const double gamma, const int maxIterLine, - const int verbose) + const int verbose, + // The line search needs to compute the fit and the gradients. Because + // the next optimization step also needs those, we can save them here and + // make them available to the caller by reference. + // The additional flag acceptedOut tells the caller if the line search + // created a non-acceptable step. + double &fitAcceptedOut, + arma::rowvec &modelGradientsAcceptedOut, + bool &acceptedOut) { static_cast(verbose); // currently not used; for later use @@ -253,6 +261,10 @@ namespace lessSEM double p_k; // new penalty value double f_k; // new combined fit + // signal to the caller that no accepted step is available yet + acceptedOut = false; + fitAcceptedOut = arma::datum::nan; + // get penalized M2LL for step size 0: double pen_0 = penalty_.getValue(parameters_kMinus1, @@ -348,6 +360,10 @@ namespace lessSEM // go to next iteration and test smaller step size continue; } + // Update the fit / gradient and acceptOut flag for the caller + fitAcceptedOut = fit_k; + modelGradientsAcceptedOut = gradients_k; + acceptedOut = true; // else break; } @@ -403,17 +419,15 @@ namespace lessSEM arma::rowvec direction(startingValues.n_elem); // prepare fit elements - // fit of the smooth part of the fit function - double fit_k = model_.fit(parameters_k, - parameterLabels) + - smoothPenalty_.getValue(parameters_k, - parameterLabels, - tuningParameters); - double fit_kMinus1 = model_.fit(parameters_kMinus1, - parameterLabels) + - smoothPenalty_.getValue(parameters_kMinus1, - parameterLabels, - tuningParameters); + // fit of the smooth part of the fit function. At this part + // fit_k and fit_kMinus1 are identical. + double fit_starting = model_.fit(startingValues, + parameterLabels) + + smoothPenalty_.getValue(startingValues, + parameterLabels, + tuningParameters); + double fit_k = fit_starting; + double fit_kMinus1 = fit_starting; // add non-differentiable part double penalizedFit_k = fit_k + penalty_.getValue(parameters_k, @@ -432,17 +446,15 @@ namespace lessSEM // prepare gradient elements // NOTE: We combine the gradients of the smooth functions (the log-Likelihood) - // of the model and the smooth penalty function (e.g., ridge) + // of the model and the smooth penalty function (e.g., ridge). As above, + // parameters_k and parameters_kMinus1 are both equal to startingValues, so the + // model gradient and smooth penalty gradient only have to be evaluated once. arma::rowvec gradients_k = model_.gradients(parameters_k, parameterLabels) + smoothPenalty_.getGradients(parameters_k, parameterLabels, tuningParameters); // ridge part - arma::rowvec gradients_kMinus1 = model_.gradients(parameters_kMinus1, - parameterLabels) + - smoothPenalty_.getGradients(parameters_kMinus1, - parameterLabels, - tuningParameters); // ridge part + arma::rowvec gradients_kMinus1 = gradients_k; // prepare Hessian elements arma::mat Hessian_k(startingValues.n_elem, startingValues.n_elem, arma::fill::zeros), @@ -472,12 +484,6 @@ namespace lessSEM Rcpp::checkUserInterrupt(); #endif - // the gradients will be used by the inner iteration to compute the new - // parameters - gradients_kMinus1 = model_.gradients(parameters_kMinus1, parameterLabels) + - smoothPenalty_.getGradients(parameters_kMinus1, parameterLabels, tuningParameters); // ridge part - - // find step direction direction = glmnetInner(parameters_kMinus1, gradients_kMinus1, Hessian_kMinus1, @@ -487,7 +493,13 @@ namespace lessSEM control_.breakInner, control_.verbose); - // find length of step in direction + // find length of step in direction. The line search already evaluates the + // smooth fit and the (model-only) gradients at the accepted step; we reuse + // them below instead of re-evaluating. + double fit_ls; + arma::rowvec modelGradients_ls(gradients_kMinus1.n_elem); + modelGradients_ls.fill(arma::datum::nan); + bool ls_accepted = false; parameters_k = glmnetLineSearch(model_, penalty_, smoothPenalty_, @@ -504,20 +516,36 @@ namespace lessSEM control_.sigma, control_.gamma, control_.maxIterLine, - control_.verbose); + control_.verbose, + fit_ls, + modelGradients_ls, + ls_accepted); - // get gradients of differentiable part - gradients_k = model_.gradients(parameters_k, - parameterLabels) + - smoothPenalty_.getGradients(parameters_k, - parameterLabels, - tuningParameters); - // fit of the smooth part of the fit function - fit_k = model_.fit(parameters_k, - parameterLabels) + - smoothPenalty_.getValue(parameters_k, - parameterLabels, - tuningParameters); + if (ls_accepted) + { + // reuse the fit and the model-only gradients from the accepted line + // search step (the smooth penalty parts are added here) + gradients_k = modelGradients_ls + + smoothPenalty_.getGradients(parameters_k, + parameterLabels, + tuningParameters); + fit_k = fit_ls; + } + else + { + // line search did not accept a step: fall back to a fresh evaluation at + // the returned parameters_k + gradients_k = model_.gradients(parameters_k, + parameterLabels) + + smoothPenalty_.getGradients(parameters_k, + parameterLabels, + tuningParameters); + fit_k = model_.fit(parameters_k, + parameterLabels) + + smoothPenalty_.getValue(parameters_k, + parameterLabels, + tuningParameters); + } // add non-differentiable part penalizedFit_k = fit_k + penalty_.getValue(parameters_k, diff --git a/include/lesstimate/ista_class.h b/include/lesstimate/ista_class.h index c3489fc..9d6fe25 100644 --- a/include/lesstimate/ista_class.h +++ b/include/lesstimate/ista_class.h @@ -213,12 +213,11 @@ namespace lessSEM arma::mat quadr, parchTimeGrad; numericVector randomNumber; // for stochastic Barzilai Borwein - // prepare fit elements +// prepare fit elements double fit_k = (1.0 / control_.sampleSize) * model_.fit(startingValues, parameterLabels) + - smoothPenalty_.getValue(parameters_k, parameterLabels, smoothTuningParameters), // ridge penalty part - fit_kMinus1 = (1.0 / control_.sampleSize) * model_.fit(startingValues, parameterLabels) + - smoothPenalty_.getValue(parameters_kMinus1, parameterLabels, smoothTuningParameters), // ridge penalty part, - penalty_k = 0.0; + smoothPenalty_.getValue(parameters_k, parameterLabels, smoothTuningParameters), + fit_kMinus1 = fit_k, + penalty_k = 0.0; double penalizedFit_k, penalizedFit_kMinus1; arma::rowvec gradients_k, gradients_kMinus1, gradient_y_k; @@ -235,14 +234,12 @@ namespace lessSEM // prepare gradient elements // NOTE: We combine the gradients of the smooth functions (the log-Likelihood) - // of the model and the smooth penalty function (e.g., ridge) + // of the model and the smooth penalty function (e.g., ridge). gradients_k = (1.0 / control_.sampleSize) * model_.gradients(parameters_k, parameterLabels) + smoothPenalty_.getGradients(parameters_k, parameterLabels, smoothTuningParameters); // ridge part - gradients_kMinus1 = (1.0 / control_.sampleSize) * model_.gradients(parameters_kMinus1, parameterLabels) + - smoothPenalty_.getGradients(parameters_kMinus1, parameterLabels, smoothTuningParameters); // ridge part + gradients_kMinus1 = gradients_k; // for acceleration: - gradient_y_k = (1.0 / control_.sampleSize) * model_.gradients(parameters_kMinus1, parameterLabels) + - smoothPenalty_.getGradients(parameters_kMinus1, parameterLabels, smoothTuningParameters); // ridge part + gradient_y_k = gradients_k; // breaking flags bool breakInner = false, // if true, the inner iteration is exited @@ -389,11 +386,16 @@ namespace lessSEM continue; } - gradients_k = (1.0 / control_.sampleSize) * model_.gradients(parameters_k, - parameterLabels) + - smoothPenalty_.getGradients(parameters_k, - parameterLabels, - smoothTuningParameters); // ridge part + // gradients_k is already computed at parameters_k when the inner + // loop converged. We recompute only if the inner loop did not converge. + if (!breakInner) + { + gradients_k = (1.0 / control_.sampleSize) * model_.gradients(parameters_k, + parameterLabels) + + smoothPenalty_.getGradients(parameters_k, + parameterLabels, + smoothTuningParameters); // ridge part + } fits(outer_iteration + 1) = penalizedFit_k; diff --git a/include/lesstimate/model.h b/include/lesstimate/model.h index 50d7696..46cbdd0 100644 --- a/include/lesstimate/model.h +++ b/include/lesstimate/model.h @@ -24,8 +24,8 @@ namespace lessSEM * @param parameterLabels stringVector with parameterLabels * @return double */ - virtual double fit(arma::rowvec parameterValues, - stringVector parameterLabels) = 0; + virtual double fit(const arma::rowvec& parameterValues, + const stringVector& parameterLabels) = 0; /** * @brief gradients method with arguments parameterValues(arma::rowvec) and parameterLabels(stringVector; see common_headers.h) * specifying the parameter values and the labels of the paramters. The function should return the gradients(arma::rowvec). @@ -35,28 +35,32 @@ namespace lessSEM * @param parameterLabels stringVector with parameterLabels * @return arma::rowvec gradients */ - virtual arma::rowvec gradients(arma::rowvec parameterValues, - stringVector parameterLabels) + virtual arma::rowvec gradients(const arma::rowvec& parameterValues, + const stringVector& parameterLabels) { - arma::rowvec gradients(parameterValues.n_elem); + // The default central-difference implementation mutates the parameter + // vector to apply +/- stepSize. With a const reference interface, take a + // local copy here; the caller's value is left untouched. + arma::rowvec perturbedParameters = parameterValues; + arma::rowvec gradients(perturbedParameters.n_elem); gradients.fill(arma::fill::zeros); // define stepSize used in numerically approximated gradients: double stepSize = 1e-5; - for (unsigned int i = 0; i < parameterValues.n_elem; i++) + for (unsigned int i = 0; i < perturbedParameters.n_elem; i++) { // step forward - parameterValues(i) += stepSize; - gradients(i) = fit(parameterValues, + perturbedParameters(i) += stepSize; + gradients(i) = fit(perturbedParameters, parameterLabels); // step backward - parameterValues(i) -= 2.0 * stepSize; - gradients(i) -= fit(parameterValues, + perturbedParameters(i) -= 2.0 * stepSize; + gradients(i) -= fit(perturbedParameters, parameterLabels); // reset - parameterValues(i) += stepSize; + perturbedParameters(i) += stepSize; // compute gradient gradients(i) /= 2.0 * stepSize;