Solver Guide#

Use this page when the equations are already on the page and the remaining question is numerical workflow: solve(), boundary_update(), boundary_search(), or sensitivity_analysis().

Read it after Modeling Guide for direct model construction, or after Getting Started once the repository examples are running. If the workflow already fails and you mainly need recovery steps, go straight to Troubleshooting.

Workflow Decision Table#

Use this workflow	When you should use it	What you get back
`solve()`	boundaries are already fixed	solved state + per-iteration error history
`boundary_update()`	the model can update boundaries from the solved grid	boundary-update state + per-iteration history
`boundary_search()`	one or more boundary values must satisfy a numerical contact/value condition	boundary-search state
`sensitivity_analysis()`	you want a path of solutions across parameter values	summary table + saved grids

`Solver(...)`: Construction Rules#

You can initialize the solver in one of two ways:

from boundary + model,
from an existing solved grid.

The most common constructor is:

solver = fjb.Solver(
    boundary=boundary,
    model=model,
    policy_guess=True,
    number=500,
    config=fjb.Config(pi_method="scan", derivative_method="central"),
)

Important options:

policy_guess: if True, trust the policy initializer; if False, force an early improvement step.
number: grid size. Larger grids are more accurate but more expensive.
config: all iteration, tolerance, and derivative settings.

How FinHJB Solves a One-Dimensional HJB#

The table above tells you which workflow to call. This section fills in the missing layer: what those workflows are actually doing to a one-dimensional HJB.

At the most abstract level, the repository solves problems of the form

\[ 0 = \sup_{\pi \in \Pi} \mathcal{H}\bigl(V(s), V_s(s), V_{ss}(s), s; \pi\bigr), \]

where the control can be scalar or vector-valued. If the model has a jump term, the implementation passes it into Model.hjb_residual(...) separately as jump.

Step 1: Discretize The Continuous Problem On The Interior Grid#

Grid.reset() first discretizes the state interval as

\[ s_0 = s_{\min} < s_1 < \cdots < s_{N-2} < s_{N-1} = s_{\max}, \]

then fixes the boundary values

\[ v_0 = v_{\text{left}}, \qquad v_{N-1} = v_{\text{right}}, \]

so the actual unknown solved by iteration is the interior vector

\[ v_{\text{inter}} = (v_1, \dots, v_{N-2}). \]

At each interior point \(s_i\), the code builds a discrete residual

\[ F_i(v_{i-1}, v_i, v_{i+1}; \pi_i) = 0. \]

For the first derivative, the current implementation uses the scheme selected by Config.dv_func, with central difference as the default:

\[ D_h v_i = \frac{v_{i+1} - v_{i-1}}{2h}. \]

For the second derivative, it uses the centered stencil

\[ D_{hh} v_i = \frac{v_{i+1} - 2v_i + v_{i-1}}{h^2}. \]

When Grid.update_with_v_inter() reconstructs the full v, dv, and d2v, it also computes boundary derivatives with second-order one-sided stencils, so the reported boundary diagnostics and the interior solve remain consistent:

v = [v_left, *v_inter, v_right]
dv = [left one-sided, dv_func(interior), right one-sided]
d2v = [left one-sided, centered interior, right one-sided]

That is why the main numerical unknown in FinHJB is not the full v array but v_inter with fixed endpoint values.

Step 2: Hold Policy Fixed And Run Policy Evaluation#

Given a current policy \(\pi\), PolicyEvaluation solves the discrete fixed-policy HJB system

\[ F(v; \pi) = 0. \]

The implementation does not treat this as a naive pointwise replacement rule. It takes Newton-style steps:

\[ J(v^{(k)}; \pi)\,\Delta v^{(k)} = -F(v^{(k)}; \pi), \qquad v^{(k+1)} = v^{(k)} + \Delta v^{(k)}. \]

The Jacobian is tridiagonal because each \(F_i\) depends only on the neighboring triple \((v_{i-1}, v_i, v_{i+1})\). In code, this maps to:

residual_pointwise(...): evaluates Model.hjb_residual(...) at one interior point;
jax.jacrev(..., argnums=(0, 1, 2)): obtains \(\partial F_i / \partial v_{i-1}\), \(\partial F_i / \partial v_i\), and \(\partial F_i / \partial v_{i+1}\);
jax.vmap(...): applies that local system across all interior points;
jax.lax.linalg.tridiagonal_solve(...): solves the Newton step.

Mechanically, the update is:

residuals, dl, d, du = vmapped_pointwise_system(grid)
dv_update = tridiagonal_solve(dl, d, du, -residuals)
grid = grid.replace(v_inter=grid.v_inter + dv_update)

EvaluationState stores the numerical status of that inner loop, including:

hjb_residuals: pointwise residuals on the interior grid;
last_update_step: the norm of the latest Newton step;
best_error and patience_counter: whether the iteration is still improving;
converged: whether the loop met pe_tol.

Stopping is defined at this layer: the loop stops when the update norm falls below pe_tol, or when improvement stalls for pe_patience rounds.

Step 3: Update Policy And Run Policy Improvement#

The outer PolicyIteration loop is:

hold the current policy fixed and run policy evaluation;
update the policy using the new v, dv, and d2v;
check whether the policy change is already small enough.

Mathematically, this is

\[ v^{k+1} = \operatorname{Eval}(\pi^k), \qquad \pi^{k+1} = \operatorname{Improve}(v^{k+1}), \]

until

\[ \max_j \lVert \pi^{k+1}_j - \pi^k_j \rVert \]

falls below pi_tol.

Implementation-wise, AbstractPolicy.update() executes two kinds of ordered policy steps:

@explicit_policy: closed-form updates that write directly into grid.policy[...];
@implicit_policy: FOCs written as local nonlinear root problems at each grid point.

For the latter, the current repository supports pointwise GaussNewton, Broyden, LevenbergMarquardt, and the custom NewtonRaphson. So policy improvement is not a hard-coded formula. It is whatever inverse mapping from value derivatives to controls the Policy class declares.

PolicyIteration currently exposes two backends:

scan: an explicit evaluation \(\rightarrow\) improvement loop with a full history of outer errors;
anderson: the same fixed-point map, but accelerated as a grid -> next_grid fixed-point problem.

In both cases, the default outer-loop error is the maximum norm across policy-array changes.

Step 4: If The Boundary Is Unknown, Wrap It As Boundary Search#

When the boundary itself is unknown, FinHJB does not merge boundary updates into policy iteration. Instead it builds an outer residual map

\[ G(b) = C\bigl(\mathrm{Solve}(b)\bigr), \]

where \(b\) is the candidate boundary vector, Solve(b) means “solve the interior HJB under that boundary,” and \(C\) comes from BoundaryConditionTarget.condition_func.

The actual _create_objective_func(...) pipeline inside boundary_search() is:

overwrite the searched boundary fields with the candidate b;
call reset() to rebuild the grid, value guess, and policy start;
run the inner HJB solver;
re-read boundary_condition() on the solved grid;
evaluate each target residual and return both the residual vector and the solved grid.

The code skeleton is:

def residual_func(boundary_params):
    boundary = initial_grid.boundary.update_boundaries(...)
    temp_grid = initial_grid.replace(boundary=boundary).reset()
    pi_state, _ = inner_func(temp_grid)
    solved_grid = pi_state.grid
    residuals = jnp.array([
        target.condition_func(solved_grid) for target in final_targets
    ])
    return residuals, solved_grid

So boundary_search() is solving an outer zero-residual boundary problem, not modifying the HJB algebra itself.

Step 5: What The Boundary Search Methods Are Actually Doing#

The current methods fall into three algorithmic roles:

bisection: scalar bracket search. With multiple targets, the implementation performs nested recursion in the order returned by boundary_condition(), so that list order is also the outer-to-inner search order. It uses each target’s own low, high, tol, and max_iter.
hybr, broyden, broyden1, krylov: vector root solvers for \(G(b)=0\), all controlled by Config.bs_tol and Config.bs_max_iter.
lm, gauss_newton, lbfgs: least-squares-style methods. The first two use the residual map in least-squares form, while lbfgs minimizes \(\sum_k G_k(b)^2\), so it is not a direct root solver in the strict sense.

So switching methods inside boundary_search() means switching how the outer boundary residual is solved, not switching how the interior HJB is discretized.

How `boundary_update()` Differs From `boundary_search()`#

Both workflows let boundaries move, but they are not solving the same outer problem.

boundary_search() solves

\[ G(b) = 0, \]

which means “find a boundary that makes a contact or smooth-pasting condition hold.”

boundary_update() is different. It does not run a root search on a residual map. It expects the model to return

boundary_dict, boundary_error = model.update_boundary(grid)

so its outer logic is closer to:

solve under the current boundary;
read off a new boundary directly from the solved grid;
use boundary_error to decide whether to continue.

If the solved grid directly implies what the next boundary should be, boundary_update() is the natural workflow. If you only know that some boundary condition must equal zero, use boundary_search().

`solve()`: Fixed-Boundary Policy Iteration#

Use:

state, history = solver.solve()

If you want the discretization and inner/outer iteration logic behind this call first, read the section above: How FinHJB Solves a One-Dimensional HJB.

Best when:

the problem has no endogenous boundary search,
you want to debug the core HJB residual before introducing another moving part,
you want the simplest possible success/failure signal.

Typical return:

state: often a PolicyIterationState,
history: a vector of iteration errors.

Useful first checks:

print(type(state).__name__)
print(history.shape)
print(state.df.head())

In a typical liquidation fixed-boundary run from this repository:

the state type is PolicyIterationState,
the history length is around a few dozen iterations,
the DataFrame columns include s, v, dv, d2v, and policy columns.

`boundary_update()`: Re-Solve While Updating Boundaries#

Use:

state, history = solver.boundary_update()

Precondition:

your model implements update_boundary(grid).

This workflow is appropriate when:

some boundary value is implied by the solved interior grid,
the boundary can be updated directly from the current solution,
you want an outer loop over “solve -> update boundary -> solve again.”

Example use in the hedging script:

locate a refinancing target m from p'(m) = 1 + gamma,
update the left boundary value from value-matching.

Useful checks:

print(type(state).__name__)
print(history.shape)
print(state.grid.boundary)

`boundary_search()`: Search for a Boundary That Satisfies a Condition#

Use:

search_state = solver.boundary_search(method="bisection", verbose=False)

If you want to see how a candidate boundary becomes a residual map for an outer search solver, read the section above: How FinHJB Solves a One-Dimensional HJB.

This is the core BCW workflow. Use it when:

one boundary value is not known in advance,
your model provides one or more BoundaryConditionTarget objects,
you want the solver to search for a value where a contact condition holds.

Supported methods:

bisection
hybr
lm
broyden
gauss_newton
lbfgs
krylov
broyden1

How The Methods Differ#

bisection is the only method that uses BoundaryConditionTarget.low, high, tol, and max_iter.
With bisection, every searched target must provide low and high.
With multi-boundary bisection, the order returned by model.boundary_condition() becomes the nested outer-to-inner search order.
hybr, lm, broyden, gauss_newton, krylov, and broyden1 treat the problem as a vector root-search problem and use Config.bs_tol and Config.bs_max_iter.
lbfgs does not solve the root problem directly. It minimizes the sum of squared residuals and is best treated as a least-squares fallback.

Practical Starting Rules#

if you have one scalar boundary target and a reliable bracket, start with bisection.
if you have two boundary targets and reliable brackets, bisection is still a sensible first default.
if you have three or more boundary targets and want a robust default, start with hybr.
if the residual map is smooth and naturally least-squares-like, try lm or gauss_newton.
if you want a quasi-Newton alternative, try broyden or broyden1.
if you only want an approximate residual minimizer, try lbfgs last.

These are implementation-level rules of thumb for the current FinHJB search backends. They are good defaults, not universal mathematical guarantees.

For finhjb-model-coder, keep one more rule in mind: if the one- or two-target bisection default fails the post-generation solve loop, the final generated code should explicitly promote the search method to hybr or another supported fallback and record that repair.

What To Inspect After Boundary Search#

state = solver.boundary_search(method="bisection", verbose=False)
grid = state.grid

print(grid.boundary)
print(grid.dv[-1], grid.d2v[-1])

For the BCW liquidation example, the high-value diagnostics are:

solved s_max rather than the initial guess,
grid.dv[-1] close to 1,
grid.d2v[-1] close to 0.

`sensitivity_analysis()`: Follow a Parameter Path#

Use:

result = solver.sensitivity_analysis(
    method="hybr",
    param_name="sigma",
    param_values=jnp.linspace(0.05, 0.20, 10),
)

This workflow is for comparative statics and continuation-style sweeps.

It returns a SensitivityResult with:

result.df: summary table across parameter values,
result.grids: a Grids container storing the solved grids.

In a small example from this repository, result.df.columns include:

sigma,
boundary_error,
converged,
s_min,
s_max,
v_left,
v_right.

That means you can analyze both:

whether the continuation succeeded numerically,
how the boundary itself moved as the parameter changed.

Configuration Tuning#

Config controls both numerical stability and runtime.

Good Default Starting Point#

For a new model, start simple:

config = fjb.Config(
    derivative_method="central",
    pi_method="scan",
    pi_max_iter=50,
    pi_tol=1e-6,
)

Why:

central is usually the safest first derivative scheme,
scan is a straightforward first policy-iteration method,
moderate tolerances tell you whether the formulation is sane before you spend more time tuning.

When Not To Use `central`#

For theory-to-code work with finhjb-model-coder, do not treat central as a universal default.

if the diffusion term becomes very small near the left boundary, prefer forward
if the diffusion term becomes very small near the right boundary, prefer backward
if diffusion stays materially positive at both edges, central remains the natural first choice

The point is not stylistic purity. The derivative scheme should reflect where the HJB becomes numerically delicate near the boundaries.

What To Tune First#

If the solve is unstable, adjust in this order:

verify the model equations and boundaries,
reduce model complexity or use a simpler initial guess,
increase number only after the base formulation is stable,
then tighten tolerances.

If boundary search is unstable:

verify the boundary target itself,
check the bracket for bisection,
inspect grid.dv[-1] and grid.d2v[-1],
only then try a different root-search method.

Common Failure Modes#

`solve()` runs but the result looks economically strange#

Do not immediately blame the solver. First inspect:

whether Policy.initialize is reasonable,
whether hjb_residual signs are correct,
whether s_min, s_max, v_left, and v_right are internally consistent.

`boundary_search()` does not settle#

Most often, one of these is wrong:

the target function is not the one you actually want,
the bracket does not contain a sign change,
the fixed-boundary solve is already unstable before search starts.

`sensitivity_analysis()` is slow#

That is normal when each point requires a full nonlinear solve. Start with a short parameter grid and expand only after you trust the path.

Solver Guide#

Workflow Decision Table#

Solver(...): Construction Rules#

How FinHJB Solves a One-Dimensional HJB#

Step 1: Discretize The Continuous Problem On The Interior Grid#

Step 2: Hold Policy Fixed And Run Policy Evaluation#

Step 3: Update Policy And Run Policy Improvement#

Step 4: If The Boundary Is Unknown, Wrap It As Boundary Search#

Step 5: What The Boundary Search Methods Are Actually Doing#

How boundary_update() Differs From boundary_search()#

solve(): Fixed-Boundary Policy Iteration#

boundary_update(): Re-Solve While Updating Boundaries#

boundary_search(): Search for a Boundary That Satisfies a Condition#

How The Methods Differ#

Practical Starting Rules#

What To Inspect After Boundary Search#

sensitivity_analysis(): Follow a Parameter Path#

Configuration Tuning#

Good Default Starting Point#

When Not To Use central#

What To Tune First#

Common Failure Modes#

solve() runs but the result looks economically strange#

boundary_search() does not settle#

sensitivity_analysis() is slow#

Related Pages#

This Page

`Solver(...)`: Construction Rules#

How `boundary_update()` Differs From `boundary_search()`#

`solve()`: Fixed-Boundary Policy Iteration#

`boundary_update()`: Re-Solve While Updating Boundaries#

`boundary_search()`: Search for a Boundary That Satisfies a Condition#

`sensitivity_analysis()`: Follow a Parameter Path#

When Not To Use `central`#

`solve()` runs but the result looks economically strange#

`boundary_search()` does not settle#

`sensitivity_analysis()` is slow#