Solving

Epsilon continuation

Let \(f_0(x,y)\) denote the modeled lower objective expression. BLVPY first defines the minimization-oriented expression

\[\begin{split} \widehat f_0(x,y)= \begin{cases} f_0(x,y),&\text{for }\texttt{cp.Minimize},\\ -f_0(x,y),&\text{for }\texttt{cp.Maximize}. \end{cases} \end{split}\]

The canonicalized lower problem is then a convex conic program in the form

\[\begin{split} \begin{array}{ll} \text{minimize} & c(x)^T u + d(x)\\ \text{subject to} & A(x)u + s = b(x)\\ & s\in \mathcal{K}, \end{array} \end{split}\]

where \(u,s\) are the lower primal and slack variables; \(c(x),d(x),A(x),b(x)\) are upper-dependent canonical data; and \(\mathcal{K}\) is a Cartesian product of cones. \(c(x)^Tu+d(x)\) represents \(\widehat f_0\), so for a modeled lower maximization its value is the negative of the original lower objective value. BLVPY solves the bilevel problem by continuation on a small relaxation of the lower KKT complementarity condition.

BLVPY introduces a dual vector \(\lambda\in \mathcal{K}^*\) and imposes

\[ A(x)^T\lambda+c(x)=0, \qquad s^T\lambda\leq\epsilon. \]

The default solve begins at \(\epsilon=10^{-1}\) and contracts it by \(0.1\) until \(10^{-6}\). Each accepted point warm-starts the next DNLP solve. If a scheduled step fails, BLVPY can retry from an intermediate epsilon, subject to max_retries.

An example of a solve call is

result = problem.solve(
    epsilon_initial=1e-1,
    epsilon_target=1e-7,
    contraction=0.1,
    feasibility_tolerance=1e-7,
    max_retries=8,
)

Deterministic initialization and restoration

The ordinary solve() path uses one deterministic upper point. If an upper variable already has a .value, BLVPY preserves it. Otherwise, for each scalar component \(x_i\) with native lower and upper bounds \(l_i\) and \(u_i\), it constructs

\[\begin{split} \widetilde{x}_i= \begin{cases} (l_i+u_i)/2, & l_i,u_i\text{ are finite},\\ l_i+1, & l_i\text{ is finite and }u_i=+\infty,\\ u_i-1, & l_i=-\infty\text{ and }u_i\text{ is finite},\\ 0, & l_i=-\infty\text{ and }u_i=+\infty. \end{cases} \end{split}\]

BLVPY first projects \(\widetilde{x}\) through the variable’s CVXPY attributes, such as nonnegativity or symmetry, to obtain \(\widehat{x}\). When the upper and generated linked-variable domain constraints define a DCP set \(\mathcal{U}\), BLVPY then attempts the least-distance projection

\[ x^{(0)}\in\mathop{\mathrm{argmin}}_{x\in\mathcal{U}} \sum_v\left\|x_v-\widehat{x}_v\right\|_F^2, \]

where the sum is over the upper variables. This projection is best effort: if it cannot be compiled or solved, BLVPY retains the attribute-projected point.

At \(x^{(0)}\), the selected conic backend solves the fixed-upper canonical lower problem

\[\begin{split} \begin{array}{ll} \text{minimize} & c(x^{(0)})^T u \\ \text{subject to} & A(x^{(0)})u+s=b(x^{(0)})\\ & s\in\mathcal{K} \end{array} \end{split}\]

with variables \(u, s\). (Here the objective offset \(d(x^{(0)})\) is irrelevant to the lower solution and is therefore omitted.) The conic solution supplies the initial canonical primal \(u^{(0)}\), slack \(s^{(0)}\), and equality dual \(\lambda^{(0)}\); BLVPY recovers the corresponding source lower variables from its affine recovery map.

If every automatic path fails, BLVPY raises an InitializationError naming variables for which an explicit .value may help.

When restoration=True (the default) and this initial point fails BLVPY’s independent residual check, BLVPY introduces a nonnegative restoration radius \(\rho\) and schematically solves

\[\begin{split} \begin{array}{ll} \mathop{\mathrm{minimize}} & \rho\\ \mathop{\mathrm{subject\ to}} & \text{each lifted constraint is relaxed by }\rho,\\ & s^T\lambda\leq\epsilon+\rho,\\ & \rho\geq0. \end{array} \end{split}\]

Restoration seeks a compatible feasible starting point; it does not optimize the upper objective. BLVPY proceeds only if the recomputed feasibility and relaxed-gap residuals are within feasibility_tolerance.

Best-of search for local solutions

With best_of=None, BLVPY follows one deterministic path. Explicit best_of=N, including N=1, generates exactly \(N\) upper initializations and runs a complete, independent continuation for every viable one. Eligible components are sampled; components controlled by an existing .value remain fixed unless sample_bounds overrides that value. Only runs that reach the target epsilon with acceptable residuals compete, and BLVPY selects the best final upper objective in its modeled sense: the lowest value for cp.Minimize and the highest value for cp.Maximize. Ties are broken by run index.

If no run reaches the target epsilon, BLVPY returns the best partial run. It first prefers the smallest attained epsilon, then the best finite upper objective in the modeled sense, and finally the lowest run index. A missing or nonfinite objective ranks after any finite objective at the same epsilon.

x.sample_bounds = (-2.0, 2.0)  # sampling metadata (not a constraint)
result = problem.solve(best_of=5, seed=42)

print(result.selected_run_index)
print(result.all_objectives)
for run in result.runs:
    print(run.index, run.initial_values[x], run.epsilon_history)

Sampling precedence is:

  1. finite sample_bounds, which override .value for that variable;

  2. an existing .value, reused in every run;

  3. finite two-sided native CVXPY bounds; otherwise

  4. blvpy.InitializationError, naming the variables that need sampling information.

sample_bounds must be a finite (lower, upper) pair broadcastable to the variable shape.

Numerical backends

IPOPT is the default DNLP backend. A different backend accepted by CVXPY’s nlp=True solve path can be passed to blvpy.BilevelProblem.solve() after proper installation, but alternative backends are not fully tested.

Clarabel is the default conic backend for fixed-upper lower initialization and upper-point projection within solve(). It is also the default for the fresh fixed-upper solves performed by polish() and gap_diagnostics().

The DNLP and conic backends used by solve() can be overridden independently:

result = problem.solve(
    solver=cp.IPOPT,
    conic_solver=cp.CLARABEL,
    solver_options={"max_iter": 500},
    conic_solver_options={"max_iter": 200},
)

BLVPY copies the option mappings and forwards them to the corresponding CVXPY solve. It uses the selected DNLP backend consistently for restoration and continuation. It uses the selected conic backend for initialization and upper projection. Availability is checked when CVXPY actually invokes the backend. polish() and gap_diagnostics() each take their own solver and solver_options; changing conic_solver on the originating solve() call does not configure either later operation.

Evaluating exponential- and 3D power-cone residuals can invoke fixed internal projection backends. These are not selected or configured by conic_solver or conic_solver_options; see Cone distance diagnostics for their numerical contract.

verbose=True (the default) prints concise BLVPY progress to standard error. solver_verbose=False (the default) suppresses CVXPY and native solver output on a best-effort basis. The flags are independent. For quiet IPOPT calls, BLVPY supplies print_level=0 and sb="yes" only when the user did not provide those options.

After a continuation solve, polish() can reuse the fixed-upper canonical lower-solve path to produce a fresh lower response. See Polishing for its feasibility and objective summary, non-mutating behavior, and solver options.