Add the homotopy operator and HomotopyProblem construction#4601
Add the homotopy operator and HomotopyProblem construction#4601BatchZero wants to merge 3 commits into
Conversation
|
Status note: per the review on SciML/SciMLBase.jl#1375 and SciML/NonlinearSolve.jl#962, the homotopy parameter λ is being separated from the parameter vector (the residual becomes |
BatchZero
force-pushed
the
feature/homotopy-lowering
branch
from
6f2690c to
4706e06
Compare
|
Both upstream dependencies are merged and registered: SciML/SciMLBase.jl#1375 → SciMLBase v3.19.0 ( |
| @@ -0,0 +1,167 @@ | |||
| # [Homotopy Initialization](@id homotopy) | |||
There was a problem hiding this comment.
Why even speak to initialization here? It's just a nonlinear solver. The fact that it shows up commonly in initialization is a detail, but there are many other uses for it.
| Problems built through the default path then receive | ||
|
|
||
| ```julia | ||
| initializealg = SciMLBase.OverrideInit(nlsolve = NonlinearSolveBase.HomotopySweep()) |
There was a problem hiding this comment.
Any nonlinear solve of equations which have a homotopy operator should default to defining a homotopy problem.
There was a problem hiding this comment.
Or rather, there should be an AbstractNonlinearProblem constructor that is automatic here.
| - If the initialization would otherwise target an `SCCNonlinearProblem`, a | ||
| homotopy-using system is forced to the monolithic `NonlinearProblem` target, | ||
| since the continuation deforms the coupled system along one `λ`. |
There was a problem hiding this comment.
It should probably split that to SCCs.
| falls back to least squares), a one-time warning is emitted and initialization | ||
| solves the `actual` (`λ = 1`) form instead. Least-squares and per-SCC sweeps | ||
| are future work. | ||
| - **Trivial initialization is silently correct.** If the initialization system |
There was a problem hiding this comment.
Obviously. Why is this even worth mentioning?
| - **`expression = Val{true}` does not auto-inject the `initializealg`.** When | ||
| generating problem-construction expressions, pass | ||
| `initializealg = SciMLBase.OverrideInit(nlsolve = NonlinearSolveBase.HomotopySweep())` | ||
| explicitly. |
| - **Scalar `Real` expressions only**, matching Modelica's restriction; the | ||
| numeric fallback method is `homotopy(::Real, ::Real)`. |
| - **Scalar `Real` expressions only**, matching Modelica's restriction; the | ||
| numeric fallback method is `homotopy(::Real, ::Real)`. | ||
| - **Failure surfaces as a failed initialization.** If the continuation cannot |
There was a problem hiding this comment.
Obviously. Also, not relevant to the MTK side.
| homotopic = has_any_homotopy(isys) | ||
| if homotopic && TProb === SCCNonlinearProblem | ||
| # A continuation sweep deforms the COUPLED system along one λ; per-SCC | ||
| # sweeping is a separate future improvement. Root ModelingToolkit's | ||
| # `get_initialization_problem_type` prefers `SCCNonlinearProblem` for split | ||
| # fully-determined systems — force the monolithic target here. | ||
| TProb = NonlinearProblem | ||
| end | ||
| initprob = TProb{_iip}(isys, op; kwargs..., build_initializeprob = false, is_initializeprob = true) | ||
| if homotopic | ||
| initprob = _maybe_sweep_homotopy_initprob(initprob, isys, _iip, kwargs) | ||
| end | ||
| return initprob |
There was a problem hiding this comment.
Split to the AbstractNonlinearProblem constructor.
| rhss = [_iszero(eq.lhs) ? eq.rhs : eq.rhs - eq.lhs for eq in eqs] | ||
| end | ||
| else | ||
| # NOTE: generate_homotopy_residual in homotopy.jl mirrors the time-independent path here; keep in sync. |
| ) | ||
| ) | ||
| end | ||
| λname = gensym(:__homotopy_λ) |
There was a problem hiding this comment.
gensym will break precompilation. Instead use a hard to hit sentinal with unicodes.
| @@ -1,3 +1,22 @@ | |||
| # Unreleased | |||
There was a problem hiding this comment.
? it'll release right away? Why write that here?
|
Do one PR first with just the HomotopyProblem constructor and direct interface, and the AbstractNonlinearProblem stuff. Then do a follow up with initialization. |
| homotopy(actual, simplified) | ||
| ``` | ||
|
|
||
| During initialization the solver may start from the easy `simplified` form and | ||
| continuously deform it into `actual`; outside initialization the operator | ||
| evaluates to `actual`, as the Modelica specification prescribes. |
There was a problem hiding this comment.
It should be more clear. The simplified equations are a set of equations for which the set of nonlinear equations is easier to get a convergent (Newton) iteration for, and the actual equations are the more complex equations you want to actually solve for. Homotopy solvers start by solving the nonlinear system with the simplified equations, and uses that solution as a starting point to solve the actual problem. There are many ways this can be used. For example if you have equations that have multiple solutions, like a quadratic equation with positive and negative solutions, you can simplify it down to an approximating linear problem which has one problem (the positive solution), and then deform it to the actual in order to stabilize the process of getting that positive solution.
There was a problem hiding this comment.
This also highlights that your implementation of homotopy has a bug since it doesn't start at lambda=0 😅
BatchZero
force-pushed
the
feature/homotopy-lowering
branch
from
4706e06 to
8cea8ba
Compare
BatchZero
changed the title
homotopy(actual, simplified) operator for initialization (spec 3.7.4)
|
Restructured into PR 1 of 2 per your review. Force-pushed as a single commit ( Addressed here (PR A):
Deferred to the initialization follow-up (PR B): the Separate, in NonlinearSolve: "It doesn't start at lambda=0… that's a bug." — Confirmed, and it's worse than robustness: the adaptive bisection only triggers on an inner-solve failure, so a bad Deferred further: "should be split out to the SCC." — The structure supports a per-SCC sweep ( MTKBase tests green locally. Expected CI reds, pre-attributed: codecov upload (fork token), the repo-wide Runic |
Implement the `homotopy(actual, simplified)` operator as a general nonlinear-solving construct: a `System` whose equations carry `homotopy` annotations is lowered into a `SciMLBase.HomotopyProblem` whose residual is the convex blend `(1 - λ)*simplified + λ*actual`, compiled as `f(u, p, λ)` and solved by natural-parameter continuation (`NonlinearSolveBase.HomotopySweep`). - systems/homotopy_operator.jl: `@register_symbolic homotopy` with the numeric fallback `homotopy(a, s) = a`; nodewise derivatives kept as opaque `homotopy(1, 0)`/`homotopy(0, 1)` nodes so jacobians and index reduction stay consistent under the blend; the `has_*homotopy` detection family; the fixed continuation sentinel `HOMOTOPY_LAMBDA` (not `gensym`, to keep codegen precompile-stable) with a collision guard; the structure-preserving `lower_homotopy` shadow pass; and `generate_homotopy_residual`. - problems/homotopyproblem.jl: `SciMLBase.HomotopyProblem(sys, op)` builds the standard residual to obtain `u0`/`p`/observed, then swaps in the swept `f(u, p, λ)`. `λ` is a trailing argument, never added to `p`. - problems/nonlinearproblem.jl: `SciMLBase.AbstractNonlinearProblem(sys, op)` auto-selects `HomotopyProblem` when `sys` contains `homotopy` nodes, otherwise `NonlinearProblem`. - systems/codegen.jl: `generate_rhs` gains an additive `extra_args` to thread `λ` as a trailing scalar argument; empty by default, leaving the standard codegen path byte-identical. - docs (basics/Homotopy.md, problems API reference) and tests for the operator, lowering, broadcasting, problem construction, auto-selection, the codegen guards, and pressure-drop parity. Initialization wiring (injecting the continuation solver as the default `initializealg`) is intentionally left to a follow-up PR.
BatchZero
force-pushed
the
feature/homotopy-lowering
branch
from
8cea8ba to
543e24b
Compare
|
@ChrisRackauckas a nudge on this one. The split you asked for is done: this is now PR 1 — the |
3 participants
Summary
Implements the
homotopy(actual, simplified)operator as a general nonlinear-solving construct.homotopyis a way to solve a hard nonlinear system by continuation: start from an easysimplifiedform and deform it into theactualone. Initialization is the most common consumer, but it is not the only one — so this PR exposes the construct on its own, and a follow-up wires it into the initialization default path.Per the review, this is PR 1 of 2:
SciMLBase.HomotopyProblem(sys, op)constructor, and the automaticSciMLBase.AbstractNonlinearProblem(sys, op)selection.initializealg = OverrideInit(nlsolve = HomotopySweep())injection,InitializationProblemrouting through the constructor,expression = Val{true}auto-injection, theNonlinearSolveBasehard dep.MTK's part is an independent lowering pass targeting an MTK-independent API. The solver stack merged and registered upstream first:
SciMLBase.HomotopyProblem(AddHomotopyProblem: a homotopy/continuation problem type SciMLBase.jl#1375, SciMLBase v3.19.0) — residualf(u, p, λ)/f(du, u, p, λ); λ is a trailing scalar argument, never an entry ofp;NonlinearSolveBase.HomotopySweep(AddHomotopySweep: natural-parameter continuation solver forHomotopyProblemNonlinearSolve.jl#962, NonlinearSolveBase v2.31.0, re-exported by NonlinearSolve) — the natural-parameter continuation solver dispatching onHomotopyProblem. AHomotopyProblemwith no algorithm defaults to it.This PR contains lowering + targeting only — no solver code.
Review items addressed
HomotopyProblem/AbstractNonlinearProblemconstruction is the headline, with initialization mentioned as one consumer (wired in PR B).AbstractNonlinearProblemconstructor that is automatic." AddedSciMLBase.AbstractNonlinearProblem(sys, op)(re-exported), which returns aHomotopyProblemwhensyscarrieshomotopynodes and aNonlinearProblemotherwise (get_nonlinear_problem_type). The explicitSciMLBase.HomotopyProblem(sys, op)is also provided.AbstractNonlinearProblemconstructor." The detection + lowering + wrap now live in theHomotopyProblem(sys, op)constructor (problems/homotopyproblem.jl), not insideInitializationProblem. PR B routes init through it.gensymbreaks precompilation. λ is now a fixed reserved unicode sentinel (__homotopy_λₘₜₖ, matching MTK's other internalₘₜₖsentinels), so codegen is a pure function of the system and the precompile / RuntimeGeneratedFunctions Expr-hash caches hit. A cheap guard turns the otherwise-silent user-symbol collision into an error.generate_rhsgains an additiveextra_args::Tuple = ()keyword;generate_homotopy_residualnow threads λ through it and collapses to two guards plus a one-line forward. The "keep in sync" comment is gone — there is no second copy. Zero behavior change for the 9 existinggenerate_rhscallers (extra_argsdefaults to(),p_endonly forwarded when non-empty).homotopy.(actual, simplified)creates one node per element; the lowering rewrites each independently under the single shared λ.How it works
homotopy(actual, simplified)is registered and opaque throughSystemconstruction,mtkcompile, and runtime codegen. The numeric fallback isactual—simplifiedis evaluated and discarded, an honest cost borne only by annotated systems. Systems that don't use the operator are byte-unaffected.homotopy(1, 0), ∂₂ =homotopy(0, 1). At runtime these fold to the derivative ofactual; under the swept rewrite they blend to λ and 1 − λ, so differentiated equations keep the same blend structure.lower_homotopy). A shadow rewrite replaces everyhomotopy(a, s)(in equations and observed) with(1 − λ)·s + λ·a, one shared sentinel λ across all occurrences (and nested ones). Structure-preserving: same equation count/order/incidence modulo the free λ. The shadow never passes throughmtkcompile.generate_homotopy_residual). Compiled asf(u, p, λ)viagenerate_rhs'sextra_args, with λ in the t-slot; an explicitp_endkeeps λ out of theMTKParametersdestructuring. Userp/MTKParameterspass through untouched. The jacobian/sparsity of the standard build are dropped (they encode the λ = 1 system).HomotopyProblem(sys, op)(orAbstractNonlinearProblem(sys, op)) builds the swept problem;λspanis inherited from the upstreamHomotopyProblemdefault(0.0, 1.0)— this PR introduces no numeric defaults of its own. Solve it withsolve(prob)(defaults toHomotopySweep) or an explicitHomotopySweep(...).Tests (MTKBase, all green locally)
homotopy_lowering.jl— opaque-through-complete, runtime fallback, nodewise derivatives, shared-λ lowering across equations + observed, nestedhomotopy, the sentinel guard, and broadcast.homotopy_problem.jl— directHomotopyProblem(sys, op)construction + residual endpoints (simplified / actual / convex blend, oop + iip), the automaticAbstractNonlinearProblemselection (and its negative control), out-of-basin rescue through the sweep, structuredppass-through, twohomotopycalls sharing one sweep.homotopy_omc_parity.jl— a pressure-drop network pattern (homotopy reaching the residual through an observed equation), solved to the analytic root (m_flow = √5to 1e-6) by bothsolve(prob)and an explicitHomotopySweep.Dependencies & compat
SciMLBasefloored at3.19(HomotopyProblem). No new hard dependency: the solver routes toHomotopySweepthroughNonlinearSolveBase's ownsolve(::HomotopyProblem)fallback at solve time, soNonlinearSolveBaseis only a test dependency here (re-exported viaNonlinearSolve). The hardNonlinearSolveBasedep arrives with the initialization injection in PR B.SimpleNonlinearSolvefloored at2.11.1(lowest co-resolvable withSciMLBase 3.19).Deferred to follow-up PRs (listed so "why just N?" is pre-answered)
initializealginjection,InitializationProblemrouting,Val{true}auto-injection).SCCNonlinearProblem.probshas no element-type constraint;HomotopyProblemsubproblems route toHomotopySweepvia the existing fallback), but it needs a small NonlinearSolve enabler first (probvec(::HomotopyProblem), asolve(::HomotopyProblem, ::AbstractNonlinearAlgorithm)wrapping rule, Stalled-resid aggregation) and per-SCC residual codegen — its own PR.expression = Val{true}for the direct constructor.Expected CI reds (pre-attributed)
FormatCheck(red on master itself; this PR adds zero new Runic violations — every touched line is Runic 1.7.0-clean, and unrelated formatting on master was left untouched to avoid drive-by churn);