Add the endogenous-grid solver family: EGM, DC-EGM, and NEGM

[hmgaudecker]

Summary

This PR adds the endogenous-grid solver family to pylcm and makes solver choice an
explicit part of a regime's economic structure:

• EGM for the classic one-dimensional consumption-savings problem;
• DCEGM for a liquid Euler margin combined with discrete choices or other
  non-Euler axes;
• NEGM for a liquid inner DC-EGM problem nested inside an outer search over a
  durable or illiquid margin;
• typed upper-envelope backends, including a certified exact backend and four
  approximate alternatives.

The implementation also introduces a constraint-routing contract. A solver must
account for every applicable constraint on every candidate-production route: it
must evaluate it, prove it by construction, compile its boundary, or reject the
model during construction.

Motivation

Grid search spends most of its work evaluating a dense continuous-action grid.
EGM reverses the Euler equation instead: it starts from an exogenous
post-decision-state grid and constructs the current-state grid on which the Euler
condition holds. This makes finely resolved consumption-savings problems much
cheaper, but only when the model exposes the economic roles and dependency
structure that the algorithm needs.

DC-EGM extends that idea to models with discrete choices. Those choices generate
overlapping, potentially non-concave value branches, so the solver must take an
upper envelope. NEGM adds a second continuous margin whose kinks or adjustment
costs make a second Euler inversion inappropriate; it retains DC-EGM for the
liquid margin and searches the outer margin directly.

Public interface

Solver and envelope names live under lcm.solvers. The specialized economic
declarations live under lcm.consumption_savings_regime:

from lcm.consumption_savings_regime import (
    ConsumptionSavingsRegime,
    LiquidMargin,
)
from lcm.solvers import DCEGM, ExactEnvelope


LIQUID = LiquidMargin(
    state="wealth",
    action="consumption",
    resources="wealth",
    post_decision_state="savings",
)

working = ConsumptionSavingsRegime(
    liquid=LIQUID,
    solver=DCEGM(
        savings_grid=savings_grid,
        envelope=ExactEnvelope(),
    ),
    states={"wealth": wealth_grid},
    actions={"consumption": consumption_grid},
    functions={
        "utility": utility,
        "savings": savings,
    },
    state_transitions={"wealth": next_wealth},
    transition=next_regime,
)

When resources are simply the liquid state, the state can fill that role
directly; an identity resources(wealth) -> wealth DAG node is unnecessary.
LiquidMargin.resources can instead name a resources function or a
NetOfAdjustmentCost declaration when the budget has another investment
decision with adjustment cost (e.g., housing).

The regime declares economic roles; the solver object contains numerical choices.
This keeps names such as wealth, consumption, and savings out of solver
configuration and lets validation explain structural mismatches at Model(...)
construction.

Regime/solver	Intended structure
Regime / GridSearch	General state-action product and arbitrary feasible set
ConsumptionSavingsRegime / EGM	Classic one-dimensional cash-on-hand problem, without discrete/process axes or an envelope
ConsumptionSavingsRegime / DCEGM	One liquid Euler margin plus discrete, passive, or stochastic axes; overlapping branches are resolved by an upper envelope
NestedConsumptionSavingsRegime / NEGM	Inner DC-EGM over the liquid margin and an outer search over a durable/illiquid post-decision grid

ConsumptionSavingsRegime and NestedConsumptionSavingsRegime also accept
GridSearch, so the same economic declaration can be used as a reference solution.
A plain Regime cannot select an endogenous-grid solver because it does not declare
the required margin roles.

Constraints and Condition

Ordinary constraint callables remain fully supported:

def can_pay_premium(cash, premium):
    return cash >= premium

For constraints whose comparison structure matters to a solver, the same Boolean
can be written as a Condition:

from lcm import implies, ref


constraints = {
    "borrowing_limit": ref("savings") >= 0.0,
    "means_test": ref("assets") < ref("asset_limit"),
    "insured_can_pay": implies(
        premise=ref("insured") == 1,
        consequent=ref("cash") >= ref("premium"),
    ),
}

Condition supports named/literal and named/named comparisons, intersections
(&), unions (|), complements (~), and implication. A reference can name a
state, action, DAG output, parameter, or declared margin role.

The callable and structured forms have the same meaning and produce the same
feasibility mask. Use Condition when:

• a solver needs to retain the named comparison in order to prove the constraint
  from its grid or compile a boundary; or
• the structured form is clearer as a declarative authoring choice.

It is not otherwise required. Condition.from_callable(...) preserves an arbitrary
callable but is intentionally opaque: it does not give a solver comparison structure
that was absent from the callable.

This distinction matters for endogenous-grid solvers because they do not construct
the same candidates as grid search. For each solver route and phase, every constraint
receives exactly one disposition:

• evaluate it where all of its inputs exist;
• prove it because the construction already enforces it;
• compile a structured boundary into candidate construction; or
• reject the model with an explanation of the unavailable inputs or unsupported
  shape.

The borrowing limit is the main example. The lower node of an EGM-family savings
grid already enforces a lower bound on the post-decision liquid state. Authors can
make that contract explicit with either a structured comparison or the margin-aware
helper:

from lcm.consumption_savings_regime import post_decision_lower_bound


constraints = {
    "borrowing_limit": post_decision_lower_bound(
        margin=LIQUID,
        lower=0.0,
    )
}

Under EGM and DCEGM, pylcm checks that the declared value equals the savings
grid's lower node and then proves the constraint by construction. Under
GridSearch, the returned ordinary predicate is evaluated normally. The declaration
is optional if the savings grid is intended to be the sole source of the borrowing
limit. An opaque callable that happens to spell the same inequality cannot be silently
treated as proven, because the solver cannot inspect what it says.

Solver contracts

The endogenous-grid solvers differentiate the model's declared utility and law of
motion rather than reconstructing a hard-coded savings equation. An analytical
inverse marginal utility function is optional; when absent, pylcm differentiates
utility and numerically inverts marginal utility.

Model construction validates the contract before solving. Among other checks, it
verifies:

• that declared margin roles resolve to continuous states/actions and assembled DAG
  nodes;
• that transition and utility dependencies are available at the stage where a solver
  needs them;
• that every constraint has a legal disposition on every candidate route;
• that the savings grid and any declared post-decision lower bound agree;
• that a nested regime has exactly the liquid and outer structure NEGM supports; and
• that backend-specific numerical configuration belongs to the selected typed
  backend.

Unsupported structures therefore fail during Model(...) construction instead of
being accepted and producing a plausible but invalid policy.

Upper envelopes

DCEGM.envelope is a typed configuration:

• ExactEnvelope certifies ownership among the finite candidates supplied to it;
• FUESEnvelope provides the Fast Upper-Envelope Scan;
• RFCEnvelope provides Rooftop Cut;
• LTMEnvelope provides a quadratic local-upper-bound reference;
• MSSEnvelope provides a HARK-style left-to-right segment sweep.

ExactEnvelope is the default. It relies on pylcm's own compiled exact-affine kernel,
built and shipped as part of the package, to make candidate-ownership decisions in
exact integer arithmetic over the stored floating-point operand bits. If that kernel
is unavailable, selecting ExactEnvelope raises a capability error rather than
falling back to floating arithmetic that cannot provide the same guarantee. Builds
that deliberately do not need the certified backend can set
LCM_SKIP_EXACT_AFFINE=1 and select one of the approximate envelope configurations.

The certification is deliberately scoped: it identifies the envelope of the finite
candidate set presented to the kernel. The resolution of a sampled constrained branch
remains controlled by DCEGM.n_constrained_points.

What this changes internally

• Adds one- and two-margin regime declarations and keeps their invariants out of the
  general Regime module.
• Adds canonical EGM continuation carries, Euler inversion, endogenous-grid sorting,
  interpolation, and policy publication.
• Carries discrete states/actions, stochastic processes, passive states, fixed
  parameters, and target-specific laws through DC-EGM's candidate construction.
• Applies the intrinsic liquid-budget mask in simulation as well as solution.
• Adds route-aware constraint planning shared by grid search and the EGM family.
• Adds the native exact-affine envelope kernel with CPU support and optional CUDA
  compilation.
• Adds benchmark models and focused analytical, brute-force-comparison, batching,
  constraint-routing, envelope, installation, and import-surface tests.

Documentation

The documentation starts from model and constraint shape, then explains which solver
fits that shape. It includes:

• the EGM / DCEGM / NEGM decision boundary;
• ordinary callable and structured Condition authoring;
• named/literal, named/named, intersection, union, complement, and implication syntax;
• the precise relationship between a post-decision lower-bound declaration and the
  savings grid;
• typed envelope selection and the scope of exact certification; and
• complete one-margin and nested consumption-savings examples.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[read-the-docs-community]

Documentation build overview

📚 pylcm | 🛠️ Build #34165719 | 📁 Comparing 0ea035c against latest (dc59124)

  🔍 Preview build  

41 files changed · + 4 added · ± 37 modified

+ Added

• adding-a-solver/index.html
• age-specialized/index.html
• credits/index.html
• epstein-zin/index.html

± Modified

• index.html
• approximating-continuous-shocks/index.html
• architecture/index.html
• benchmarking/index.html
• benchmarking-1/index.html
• beta-delta/index.html
• continuous-stochastic-processes/index.html
• conventions/index.html
• debugging/index.html
• defining-models/index.html
• and 27 more...

[github-actions]

Benchmark comparison (main → HEAD)

Comparing 9062f58e (main) → 0ea035c7 (HEAD)

Benchmark	Statistic	before	after	Ratio	Alert
aca-baseline	execution time	9.696 s	10.058 s	1.04
	peak GPU mem	603 MB	604 MB	1.00
	compilation time	439.00 s	440.78 s	1.00
	peak CPU mem	8.93 GB	10.34 GB	1.16	❌
aca-baseline-debug	execution time	39.252 s	36.291 s	0.92
	peak GPU mem	603 MB	602 MB	1.00
	compilation time	491.38 s	490.12 s	1.00
	peak CPU mem	6.94 GB	6.53 GB	0.94
Mahler-Yum	execution time	4.634 s	4.922 s	1.06
	peak GPU mem	788 MB	788 MB	1.00
	compilation time	15.91 s	15.67 s	0.99
	peak CPU mem	1.51 GB	1.35 GB	0.89
Precautionary Savings - Solve	execution time	25.7 ms	26.2 ms	1.02
	peak GPU mem	34 MB	34 MB	1.00
	compilation time	2.43 s	2.36 s	0.97
	peak CPU mem	1.09 GB	1.07 GB	0.98
Precautionary Savings - Simulate	execution time	70.1 ms	62.2 ms	0.89
	peak GPU mem	328 MB	328 MB	1.00
	compilation time	4.71 s	5.14 s	1.09
	peak CPU mem	1.27 GB	1.23 GB	0.97
Precautionary Savings - Solve & Simulate	execution time	101.7 ms	100.5 ms	0.99
	peak GPU mem	1.68 GB	1.68 GB	1.00
	compilation time	7.08 s	7.13 s	1.01
	peak CPU mem	1.30 GB	1.25 GB	0.97
Precautionary Savings - Solve & Simulate (irreg)	execution time	236.9 ms	233.9 ms	0.99
	peak GPU mem	2.20 GB	2.20 GB	1.00
	compilation time	7.28 s	7.46 s	1.03
	peak CPU mem	1.37 GB	1.32 GB	0.96
IskhakovEtAl2017DCEGMSimulate	execution time		211.0 ms
	compilation time		7.43 s
	peak CPU mem		1.46 GB
IskhakovEtAl2017DCEGMSolve	execution time		1.891 s
	compilation time		7.24 s
	peak CPU mem		1.37 GB
IskhakovEtAl2017Simulate	execution time	227.5 ms	230.0 ms	1.01
	compilation time	7.85 s	7.97 s	1.02
	peak CPU mem	1.28 GB	1.22 GB	0.95
IskhakovEtAl2017Solve	execution time	50.7 ms	54.5 ms	1.08
	compilation time	2.69 s	2.65 s	0.99
	peak CPU mem	1.12 GB	1.11 GB	0.99
IskhakovEtAl2017DCEGMSimulateGpuPeakMem	peak GPU mem		440 MB
IskhakovEtAl2017DCEGMSolveGpuPeakMem	peak GPU mem		0 MB
IskhakovEtAl2017SimulateGpuPeakMem	peak GPU mem	440 MB	440 MB	1.00
IskhakovEtAl2017SolveGpuPeakMem	peak GPU mem	67 MB	67 MB	1.00

[mj023]

Great additions, I was obviusly not able to check everything super thoroughly, but the algorithms look interesting and most of it seems to be well tested againt multiple benchmarks. Just some general thoughts:

General Structure:

I think the structure of the EGM folder is a bit cluttered, if possible, I think it would be better if the solvers could be split up so that each one has their own submodule.

Solvers:

As I understand it, the solvers are all specialized for slightly different models, which is fine for now. I feel like there is some overlap between the solvers, for example the one_asset_egm_step could reasonably be solved by the standard DC-EGM. I think we should overall try to find the most general solver that works well on GPU and focus on that for further development, maybe even remove some of the others later.

Upper Envelope Kernels:

I think here the same applies as for the solvers, we should try to at some point determine which methods are good and only support these. A good method would ideally return a constant number of gridpoints and be easily parallelizable to work on the GPU. The best fitting here will probably be the interpolation based algorithms like the query-algorithm or ltm-algorithm (which are nearly the same on the first look?). I am now quite sure, that both can also be implemented in a way, such that they are on average close to $\Theta(n)$ and can still be parallelized (we would need to write the CUDA kernel ourselves).

[codecov]

Codecov Report

❌ Patch coverage is 90.28444% with 304 lines in your changes missing coverage. Please review.
✅ Project coverage is 91.42%. Comparing base (9062f58) to head (0ea035c).

Files with missing lines	Patch %	Lines
src/_lcm/egm/continuation.py	83.65%	68 Missing ⚠️
src/_lcm/egm/upper_envelope/certified_sign.py	63.93%	44 Missing ⚠️
hatch_build.py	72.80%	31 Missing ⚠️
benchmarks/asv/bench_aca_baseline.py	0.00%	27 Missing ⚠️
benchmarks/asv/bench_iskhakov_et_al_2017.py	67.14%	23 Missing ⚠️
src/_lcm/constraints/ir.py	90.38%	20 Missing ⚠️
src/_lcm/egm/upper_envelope/_exact_affine/ffi.py	85.27%	19 Missing ⚠️
src/_lcm/egm/kernel_scope.py	84.21%	15 Missing ⚠️
benchmarks/publish.py	0.00%	12 Missing ⚠️
src/_lcm/egm/upper_envelope/__init__.py	92.63%	7 Missing ⚠️
... and 14 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #390      +/-   ##
==========================================
+ Coverage   91.25%   91.42%   +0.16%     
==========================================
  Files         177      249      +72     
  Lines       16956    24630    +7674     
==========================================
+ Hits        15473    22517    +7044     
- Misses       1483     2113     +630     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[timmens]

Very nice!

1. Plain EGM bypasses its lower-bound proof
   (src/_lcm/solution/egm.py): validate_model rejects every constraint before
   route planning, including a matching post_decision_lower_bound. Let
   proof-carrying constraints reach the planner and reject only unmet ones.

2. Condition is broader than current use (src/lcm/condition.py,
   src/_lcm/constraints/): the only shipped structural consumer recognizes one
   numeric lower bound; no boundary compiler exists. A typed LowerBound or
   InequalityConstraint would cover this with much less complexity. What concrete
   second consumer justifies the Boolean DSL?

3. name_in_dag leaks internals
   (src/lcm/consumption_savings_regime.py): rename it to output or name, or
   infer it.

4. Regime.solver overstates support (src/lcm/regime.py): it lists EGM,
   DCEGM, and NEGM, although these require specialized regime classes. Document the
   pairing explicitly.

[hmgaudecker]

Thanks! 1, 3, 4 are implmented.

2. **`Condition` is broader than current use** (`src/lcm/condition.py`,
   `src/_lcm/constraints/`): the only shipped structural consumer recognizes one
   numeric lower bound; no boundary compiler exists. A typed `LowerBound` or
   `InequalityConstraint` would cover this with much less complexity. What concrete
   second consumer justifies the Boolean DSL?

The genesis here is that we needed it in #400, but then introduced it here so we don't introduce a transient interface. #400 already contains the substantive boundary-processing algorithm: it maps declared case and piecewise-affine metadata into asset-axis breakpoints and interval-specific EGM problems. But your comment was constructive in the sense that some elements have not been wired in yet, will do so now!