Add Wolfe line search to Laplace approximation by SteveBronder · Pull Request #3250 · stan-dev/math · GitHub

Summary

This PR makes the following changes for the laplace approximation:

Adds a wolfe line search to the Newton solver used in the laplace approximation to improve convergence.

The example code provided from Laplace Bug when passing Eigen::Map in tuple of functor arguments #3205 fails on develop. The issue arose that the initial value of 0 for theta started the model in the tail of the distribution. The quick line search we did which only tested half of a newton step was not robust enough for this model to reach convergance. This PR adds a full wolfe line search to the Newton solver used in the laplace approximation to improve convergence in such cases.
The graphic below shows the difference in estimates of the log likelihood for laplace relative to integrate_1d on the roach test data plotted along the mu and sigma estimates. There is still a bias relative to integrate_1d as mu becomes negative and sigma becomes larger, but it is much nicer than before.

The main loop for laplace_marginal_density_est is expensive as it requires calculating either a diagonal hessian or block diagonal hessian with 2nd order autodiff. The wolfe line search only requires the gradients of the likelihood with respect to theta. So with that in mind the wolfe line search tries pretty aggressively get the best step size. If our initial step size is successful, we try to keep doubling until we hit a step size where the strong wolfe conditions fail and then return the information for the step right before that failure. If our initial step size does not satisfy strong wolfe then we do a bracketed zoom with cubic interpolation until till we find a step size that satisfies the strong wolfe conditions.
Tests for the wolfe line search are added to test/unit/math/laplace/wolfe_line_search.hpp.

Fixes bugs in the laplace approximation

Fix iteration mismatch between values when line search succeeds
In the last iteration of the laplace approximation we were returning the negative block diagonal hessian and derived matrices from the previous search. This is fine if the line search in that last step failed. But if the line search succeeds then we need to go back and recalculate the negative block diagonal hessian and it's derived quantities.
Breakup diagonal and block hessian functions
Previously we had one block_hessian function that calculated both the block hessian or the diagonal hessian at runtime. But this function is only used in places where we know at compile time whether we want a block or diagonal hessian. So I split out the two functions to avoid unnecessary runtime branching.
barzilai_borwein_step_size
For an initial step size estimate before each line search we use the Barzilai-Borwein method to get an estimate.
Adjoints of ll args only calculated once
Previously we calculated them eargerly in each laplace iteration. But they are not needed within the inner loop so we wait till we finish the inner search then calculate their adjoints once afterwards.
Calculate covariance once at the start and reuse throughout.
We were calculating the covariance matrix from inside of laplace_density_est, but this required us to then return it from that function and imo looked weird. So I pulled it out and now laplace_marginal_density_est is passed the covariance matrix.

Fixes numerical stability in laplace distributions
There were a few places where we could use log_sum_exp etc. so I made those changes.
Fixes "bug" in finite difference step size calculation

Changed from cube root of epsilon to epsilon^(1/7) for 6th order
The finite difference method in Stan was previously using stepsize optimzied a 2nd order method. But the code is a 6th order method. I modified finite_diff_stepsize to use epsilon^(1/7) instead of cbrt(epsilon). With this change all of the laplace tests pass with a much higher tolerance for precision.

Tests

All the AD tests now have a tighter tolerance for the laplace approximation.
There are also tests for the wolfe line search in test/unit/math/laplace/wolfe_line_search.hpp.

./runTests.py test/unit/math/laplace

Release notes

Improve laplace approximation with wolfe line search and bug fixes.

Checklist

Copyright holder: Steve Bronder

The copyright holder is typically you or your assignee, such as a university or company. By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing
- unit tests pass (to run, use: ./runTests.py test/unit)
- header checks pass, (make test-headers)
- dependencies checks pass, (make test-math-dependencies)
- docs build, (make doxygen)
- code passes the built in C++ standards checks (make cpplint)
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested

…ch method

…line-search

… the lambdas

…lues for W, B, etc. are used

…order stepsize

@SteveBronder the latest couple commits here look like mistakes meant for #3271

@SteveBronder f993bdb didn't fix the issue in stan-dev/stanc3#1586. You need to add a decay to the result of tuple_element_t in is_inner_tuple_type_v to get the intended behavior when the tuple of options contains e.g. const double&

Thanks for resolving all the above. I think there was one missed, and a few new things from the changes since then

The idea is this could be exposed to the user, right? We probably want to workshop the name if so

I'm open to new names!

I also ran the tests locally and noticed that nearly as many were being marked as SKIPPED as actually ran -- is that expected?

I also ran the tests locally and noticed that nearly as many were being marked as SKIPPED as actually ran -- is that expected?

Yes. Only a few tests run without the wolfe step, but we skip tests where having turning off the wolfe line search causes a failure. This is expected as the wolfe step stops overshooting from a bad initial position.

fyi jenkins is failing just because of the api change

--- Compiling C++ code ---
clang++-6.0 -std=c++17 -D_REENTRANT -Wno-sign-compare -Wno-ignored-attributes      -I stan/lib/stan_math/lib/tbb_2020.3/include    -O0 -I src -I stan/src -I stan/lib/rapidjson_1.1.0/ -I lib/CLI11-1.9.1/ -I stan/lib/stan_math/ -I stan/lib/stan_math/lib/eigen_3.4.0 -I stan/lib/stan_math/lib/boost_1.87.0 -I stan/lib/stan_math/lib/sundials_6.1.1/include -I stan/lib/stan_math/lib/sundials_6.1.1/src/sundials    -DBOOST_DISABLE_ASSERTS          -c -include-pch stan/src/stan/model/model_header.hpp.gch/model_header_4_2.hpp.gch -x c++ -o ../stanc3/test/integration/good/code-gen/laplace_nested_tuple3.o ../stanc3/test/integration/good/code-gen/laplace_nested_tuple3.hpp
../stanc3/test/integration/good/code-gen/laplace_neg_binomial_2_log.hpp:404:24: error: no matching function for call to 'laplace_marginal_tol_neg_binomial_2_log_lpmf'
        lp_accum__.add(stan::math::laplace_marginal_tol_neg_binomial_2_log_lpmf<
                       ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
../stanc3/test/integration/good/code-gen/laplace_neg_binomial_2_log.hpp:792:12: note: in instantiation of function template specialization 'laplace_neg_binomial_2_log_model_namespace::laplace_neg_binomial_2_log_model::log_prob_impl<false, false, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, nullptr, nullptr, nullptr>' requested here
    return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
           ^

So I think once we approve this and the stanc3 pr then we can just merge both at the same time

comments from meeting with @charlesm93

Consider logging when jitter is needed -- might be too verbose if happens every time?

I think this would be too verbose to run every time.

Fix string for general usage

Test cutting this out

@charlesm93 So actually there are a few packages that check this and if it the direction is not correct they throw an error

LBFGSpp just hard throws if the descent direction is not correct

https://github.com/yixuan/LBFGSpp/blob/master/include/LBFGSpp/LineSearchBacktracking.h#L62

Ceres conjugate gradient solver flips the direction

https://github.com/ceres-solver/ceres-solver/blob/master/internal/ceres/line_search_direction.cc#L87

So I'm not sure what we should do here. We could just leave it?

I'm going to leave this.

In the wolfe line search we are doing maximization and want the low bracket to have a positive gradient and the high bracket to have a negative gradient so we are searching a space that looks like /\. If we don't flip this then we end up in a starting position that can look like \(/|\) which would goof up the search. So I'm going to leave this

Delete

Delete !wolfe since if wolfe is true we would have completed earlier

Check this is correct because we are looking for a low with directional derivative > 0 and high with directional derivative < 0.

This can be correct if we are trying to exit early here and choose to start with a bisection and hope later in the bisection scheme we either find a good point or a point where we can do the cubic interpolation (`log.g > 0 and high.g < 0)

Check if necessary

Suggested change

      
                // Choose trial alpha: cubic when bracket is good, else bisection.
          
                double alpha_mid = use_cubic ? cubic_or_bisect_max(low, high, opt)
          
                                             : 0.5 * (low.alpha() + high.alpha());
          
                // Choose trial alpha: cubic when bracket is good, else bisection.
          
                double alpha_mid = cubic_or_bisect_max(low, high, opt);

This reverts commit 17002d9.

This reverts commit 495e9c6.

…ntu1" This reverts commit 97846a6.

This reverts commit eab35ec.

I've started a re-run of downstream against the updated stanc3 here: https://jenkins.flatironinstitute.org/blue/organizations/jenkins/Stan%2FStan/detail/downstream_tests/780/pipeline

I don't think this will actually update the build status here (might stay a red x), but it should be sufficient if it passes

So we are going to merge this PR, then stanc. Note that once this is merged the CI will fail for any Stan pull requests for about an hour or so, but we can just rerun that stuff.

SteveBronder added 30 commits September 29, 2025 18:03

update the laplace line search to use a more advanced wolfe line sear…

43c3eef

…ch method

Merge remote-tracking branch 'origin/develop' into fix/laplace-wolfe-…

74a92bc

…line-search

add data for roach data test

2de7a4f

update tests

547288e

move wolfe to its own file

cb5e282

Merge remote-tracking branch 'origin' into fix/laplace-wolfe-line-search

cdaf700

update wolfe

6b22c85

Merge remote-tracking branch 'origin' into fix/laplace-wolfe-line-search

a1d0906

update to use barzilai borwein step size as initial step size estimate

5b6ffff

seperate moto from other lpdf tests

8eff766

update

f542cc5

add WolfeInfo

c845944

use WolfeInfo for extra data

40f1243

put everything for iterations in laplace into structs

6e528d2

update poisson test

d89eeb5

add swap functions

40d889f

cleanup laplace_density_est to reduce repeated code

59b7a2f

update to search for a good initial alpha on a space

c73f5aa

fix code for wolfe line search

773d417

update tests for zoom

b557dad

all tests pass for laplace with new wolfe

b18bf87

use log sum of diagonal of U matrix for solver 3 determinant

98df588

move update_step to be a user passed function

929dd47

cleanup the laplace code to remove some passed by reference values to…

2ebb01a

… the lambdas

cleanup the laplace code to remove some passed by reference values to…

3bbcef3

… the lambdas

update WolfeData with member accessors and use Eval within WolfeData

66ffec9

update docs for wolfe

ff5bee4

update logic in laplace_marginal_desntiy_est so that final updated va…

7a7415a

…lues for W, B, etc. are used

clang format

973144a

change stepsize of finite difference to use 6th order instead of 2nd …

cc5d49a

…order stepsize

WardBrian mentioned this pull request Jan 26, 2026

Update laplace _tol functions for allow_fallthrough argument, tuple of control params stan-dev/stanc3#1586

Merged

2 tasks

use std::decay on tuple_to_laplace_options template argument

f993bdb

SteveBronder and others added 3 commits January 27, 2026 14:20

update with decay_t

102a79a

Merge commit '42063be9833129cdf1d4b2299616ec8aba237d80' into HEAD

0acaf73

[Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1

96e3c2f

WardBrian reviewed Jan 28, 2026

View reviewed changes

SteveBronder and others added 2 commits January 29, 2026 15:17

fixes for @WardBrian review

6288ccd

[Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1

8cefcd4

Remove commented out debugging statements in aki_ex

6789f32

WardBrian reviewed Feb 3, 2026

View reviewed changes

SteveBronder and others added 11 commits February 5, 2026 17:58

update with last review

eab35ec

[Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1

97846a6

fix wolfe table formatting

495e9c6

remove log files

17002d9

Revert "remove log files"

37d2ab4

This reverts commit 17002d9.

Revert "fix wolfe table formatting"

30573a7

This reverts commit 495e9c6.

Revert "[Jenkins] auto-formatting by clang-format version 10.0.0-4ubu…

237f2a4

…ntu1" This reverts commit 97846a6.

Revert "update with last review"

cd363fe

This reverts commit eab35ec.

cleanup tests and address review comments

7a47143

[Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1

fcd1b3c

changed tests compile with O3 and native

ca363fe

SteveBronder merged commit 3ba26b8 into develop Feb 11, 2026
33 of 34 checks passed

WardBrian mentioned this pull request Apr 21, 2026

Laplace tests taking too long in develop Jenkins runs #3309

Open

+               * @param theta_0_size Size of user supplied initial theta
+               * @return tuple representing laplace options exposed to user.
+               */
+              inline auto generate_laplace_options(int theta_0_size) {

+               * @param max_jitter Maximum jitter to try
+               */
+              template <typename LLT, typename B_t>
+              inline void llt_with_jitter(LLT& llt_B, B_t& B, double min_jitter = 1e-10,

+                                             + covariance * state.prev().theta_grad();
+                    state.wolfe_info.init_dir_ = state.prev_g.dot(state.wolfe_info.p_);
+                    // Flip direction if not ascending
+                    state.wolfe_info.flip_direction();

+                    } else if (options.line_search.max_iterations == 0) {
+                      state.curr().update(scratch);
+                      state.wolfe_status.accept_ = true;
+                      finish_update = false;

+                  if (!finite_ok) {  //   f or g is NaN/Inf → shrink
+                    high.alpha() *= 0.5;
+                    high_has_eval = false;
+                    if (high.alpha() < opt.min_alpha) {
+                      break;
+                    }
+                    continue;
+                  }

+                  }
+                  if (llt_B.info() != Eigen::Success) {
+                    throw std::domain_error(
+                        "laplace_marginal_density: Cholesky (Diag) failed");

+                      state.wolfe_status = internal::wolfe_line_search(
+                          state.wolfe_info, update_fun, options.line_search, msgs);
+                    }
+                    if (run_convergence_check) {

+                    best = high;
+                  }
+                  const bool dir_pos = high.dir() > 0;
+                  if (armijo && !wolfe && dir_pos) {  // [2]

+                if (!high_has_eval) {
+                  wolfe_check = update_with_tick(scratch, curr, prev, high, p);
+                  if (wolfe_check.stop_ != WolfeReturn::Continue) {
+                    return wolfe_check;
+                  }
+                  high_has_eval = true;
+                }

		// [3,4,5]
		found_right = true;

Add Wolfe line search to Laplace approximation by SteveBronder · Pull Request #3250 · stan-dev/math · GitHub

Uh oh!

Uh oh!

Add Wolfe line search to Laplace approximation - #3250

Add Wolfe line search to Laplace approximation#3250SteveBronder merged 125 commits intodevelopstan-dev/math:developfrom fix/wolfe-zoom1stan-dev/math:fix/wolfe-zoom1Copy head branch name to clipboard

Conversation

SteveBronder commented Oct 24, 2025

Summary

Tests

Release notes

Checklist

Uh oh!

WardBrian commented Jan 26, 2026

Uh oh!

WardBrian commented Jan 27, 2026

Uh oh!

WardBrian left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

WardBrian commented Jan 28, 2026

Uh oh!

SteveBronder commented Jan 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

SteveBronder commented Jan 30, 2026

Uh oh!

WardBrian left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

WardBrian commented Feb 10, 2026

Uh oh!

SteveBronder commented Feb 11, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Add Wolfe line search to Laplace approximation#3250
SteveBronder merged 125 commits into
developfrom
fix/wolfe-zoom1

SteveBronder commented Jan 29, 2026 •
edited
Loading