Struct Solver 

pub struct Solver { /* private fields */ }

Session-style solver: holds an IpoptApplication, its TNLP, and the converged factor between calls.

Implementations

impl Solver

pub fn new(app: IpoptApplication, tnlp: Rc<RefCell<dyn TNLP>>) -> Self

Build a new session. The app should already have its options configured and initialize() called.

Examples found in repository

examples/sensitivity_session.rs (line 166)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

examples/parametric_mpc.rs (line 178)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    // One full IPM solve at the nominal parameter.
    let mut solver = Solver::new(make_app(), tnlp);
    let t0 = Instant::now();
    let status = solver.solve();
    let solve_dt = t0.elapsed();
    println!(
        "IPM solve at nominal (eta1=5.0, eta2=1.0): status={status:?}, {:.3} ms",
        solve_dt.as_secs_f64() * 1e3
    );
    assert!(solver.converged().is_some());

    // 10 parametric steps as eta2 sweeps from 1.0 to 1.5. Each one is
    // a back-solve against the held factor — no IPM iteration.
    let pins = vec![2 as Index, 3];
    let mut total_step_dt = std::time::Duration::ZERO;
    let mut steps = Vec::new();
    for k in 1..=10 {
        let d_eta2 = 0.05 * k as f64;
        let t = Instant::now();
        let dx = solver
            .parametric_step(&pins, &[0.0, d_eta2])
            .expect("parametric_step ok");
        total_step_dt += t.elapsed();
        steps.push((d_eta2, dx));
    }

    println!("\nparametric steps against the held factor:");
    for (d_eta2, dx) in &steps {
        println!(
            "  Δeta2 = {d_eta2:+.2}  -> Δx_primal = [{:+.5}, {:+.5}, {:+.5}]",
            dx[0], dx[1], dx[2]
        );
    }
    let avg_step_us = total_step_dt.as_secs_f64() * 1e6 / steps.len() as f64;
    println!(
        "\n10 parametric steps: total {:.3} ms, mean {avg_step_us:.1} µs/step",
        total_step_dt.as_secs_f64() * 1e3
    );
    println!(
        "\nRatio: each parametric step is roughly {:.0}x cheaper than the IPM solve.",
        solve_dt.as_secs_f64() / (total_step_dt.as_secs_f64() / steps.len() as f64),
    );
}

examples/sensitivity_factor_reuse_bench.rs (line 167)

fn main() {
    let deltas: Vec<Vec<Number>> = (1..=N).map(|k| vec![0.0, 0.01 * k as Number]).collect();
    let pins = vec![2 as Index, 3];

    // (1) Held-factor path: 1 solve + N parametric steps.
    let t = Instant::now();
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));
    let mut solver = Solver::new(make_app(), tnlp);
    solver.solve();
    let solve_dt = t.elapsed();
    let mut held_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let dx = solver
            .parametric_step(&pins, d)
            .expect("parametric_step ok");
        held_steps.push(dx);
    }
    let held_step_dt = t.elapsed();
    let held_total = solve_dt + held_step_dt;

    // (2) Cold path: N fresh SensSolve runs (one full IPM each).
    let mut cold_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
            eta1: 5.0,
            eta2: 1.0,
        }));
        let mut app = make_app();
        let r = SensSolve::new(pins.clone())
            .with_deltas(d.clone())
            .run(&mut app, tnlp);
        cold_steps.push(r.dx.expect("dx"));
    }
    let cold_total = t.elapsed();

    // Cross-check: held vs cold dx should agree.
    let mut max_err = 0.0_f64;
    for (h, c) in held_steps.iter().zip(cold_steps.iter()) {
        for (a, b) in h.iter().zip(c.iter()) {
            max_err = max_err.max((a - b).abs());
        }
    }

    println!("Workload: 1 IPM solve + {N} parametric steps (Δeta2 sweep).");
    println!();
    println!(
        "Held-factor path : solve {:.2} ms + {N} steps {:.2} ms = {:.2} ms total",
        solve_dt.as_secs_f64() * 1e3,
        held_step_dt.as_secs_f64() * 1e3,
        held_total.as_secs_f64() * 1e3,
    );
    println!(
        "Cold-restart path:                          {N} fresh solves = {:.2} ms total",
        cold_total.as_secs_f64() * 1e3,
    );
    println!();
    println!(
        "Speedup of held-factor over cold-restart: {:.1}x",
        cold_total.as_secs_f64() / held_total.as_secs_f64()
    );
    println!("Numerical agreement (held vs cold dx): max |err| = {max_err:.2e}");
}

examples/rh_orientation_check.rs (line 143)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(PinnedQuadratic { p0: 1.0, p1: 2.0 }));
    let mut app = IpoptApplication::new();
    app.options_mut()
        .set_integer_value("print_level", 0, true, false)
        .unwrap();
    app.options_mut()
        .set_string_value("sb", "yes", true, false)
        .unwrap();
    app.initialize().unwrap();

    let mut solver = Solver::new(app, tnlp);
    let status = solver.solve();
    assert!(
        matches!(
            status,
            ApplicationReturnStatus::SolveSucceeded
                | ApplicationReturnStatus::SolvedToAcceptableLevel
        ),
        "solve failed: {status:?}"
    );

    let (hr, vals, vecs) = solver
        .compute_reduced_hessian_eigen(&[0, 1], 1.0)
        .expect("reduced Hessian");

    println!("model:  min x0² + x1² + x0·x1  s.t.  x0 = 1, x1 = 2");
    println!("H    = [[2, 1], [1, 2]]   (eigenvalues 1 and 3)");
    println!();
    println!("compute_reduced_hessian(pins=[0, 1]):");
    for i in 0..2 {
        println!("  [{:>9.6}, {:>9.6}]", hr[i], hr[i + 2]);
    }
    println!(
        "  eigenvalues (ascending) = [{:>9.6}, {:>9.6}]",
        vals[0], vals[1]
    );
    for j in 0..2 {
        println!(
            "  eigenvector[{j}]          = [{:>9.6}, {:>9.6}]",
            vecs[2 * j],
            vecs[2 * j + 1]
        );
    }
    println!();

    // Three candidates, discriminated by magnitude as well as sign:
    // inv(H) = [[2/3, -1/3], [-1/3, 2/3]].
    let h_inv: [Number; 4] = [2.0 / 3.0, -1.0 / 3.0, -1.0 / 3.0, 2.0 / 3.0];
    let candidates: [(&str, [Number; 4]); 3] = [
        ("+H_R", H),
        ("-H_R", [-H[0], -H[1], -H[2], -H[3]]),
        ("H_R⁻¹", h_inv),
    ];
    for (name, want) in candidates {
        let err = (0..4)
            .map(|k| (hr[k] - want[k]).abs())
            .fold(0.0 as Number, Number::max);
        println!(
            "  vs {name:<6} max|Δ| = {err:.3e}  {}",
            if err < 1e-7 { "← MATCH" } else { "" }
        );
    }
    println!();
    println!(
        "So the ascending spectrum runs STIFFEST first: {:.6} is the curvature-3",
        vals[0]
    );
    println!(
        "mode and {:.6} the curvature-1 (soft) one — the reverse of the",
        vals[1]
    );
    println!("order a caller reading `+H_R` would assume.");
}

pub fn app(&self) -> &IpoptApplication

Borrow the underlying IpoptApplication (e.g. to read its options table after a solve). Mutation between solve calls is supported via Self::app_mut.

pub fn app_mut(&mut self) -> &mut IpoptApplication

Mutable borrow of the underlying IpoptApplication. Useful for reconfiguring options before a follow-up solve(). Note that changing options that affect the KKT linear system between calls will invalidate the cached factor; the next solve() rebuilds it.

pub fn solve(&mut self) -> ApplicationReturnStatus

Run the IPM to convergence. On a successful solve the ConvergedState (including the KKT backsolver) is stashed inside the Solver and accessible via Self::converged.

Each call to solve() overwrites the previous converged state; the previously held factor is dropped.

Examples found in repository

examples/sensitivity_session.rs (line 167)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

examples/parametric_mpc.rs (line 180)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    // One full IPM solve at the nominal parameter.
    let mut solver = Solver::new(make_app(), tnlp);
    let t0 = Instant::now();
    let status = solver.solve();
    let solve_dt = t0.elapsed();
    println!(
        "IPM solve at nominal (eta1=5.0, eta2=1.0): status={status:?}, {:.3} ms",
        solve_dt.as_secs_f64() * 1e3
    );
    assert!(solver.converged().is_some());

    // 10 parametric steps as eta2 sweeps from 1.0 to 1.5. Each one is
    // a back-solve against the held factor — no IPM iteration.
    let pins = vec![2 as Index, 3];
    let mut total_step_dt = std::time::Duration::ZERO;
    let mut steps = Vec::new();
    for k in 1..=10 {
        let d_eta2 = 0.05 * k as f64;
        let t = Instant::now();
        let dx = solver
            .parametric_step(&pins, &[0.0, d_eta2])
            .expect("parametric_step ok");
        total_step_dt += t.elapsed();
        steps.push((d_eta2, dx));
    }

    println!("\nparametric steps against the held factor:");
    for (d_eta2, dx) in &steps {
        println!(
            "  Δeta2 = {d_eta2:+.2}  -> Δx_primal = [{:+.5}, {:+.5}, {:+.5}]",
            dx[0], dx[1], dx[2]
        );
    }
    let avg_step_us = total_step_dt.as_secs_f64() * 1e6 / steps.len() as f64;
    println!(
        "\n10 parametric steps: total {:.3} ms, mean {avg_step_us:.1} µs/step",
        total_step_dt.as_secs_f64() * 1e3
    );
    println!(
        "\nRatio: each parametric step is roughly {:.0}x cheaper than the IPM solve.",
        solve_dt.as_secs_f64() / (total_step_dt.as_secs_f64() / steps.len() as f64),
    );
}

examples/sensitivity_factor_reuse_bench.rs (line 168)

fn main() {
    let deltas: Vec<Vec<Number>> = (1..=N).map(|k| vec![0.0, 0.01 * k as Number]).collect();
    let pins = vec![2 as Index, 3];

    // (1) Held-factor path: 1 solve + N parametric steps.
    let t = Instant::now();
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));
    let mut solver = Solver::new(make_app(), tnlp);
    solver.solve();
    let solve_dt = t.elapsed();
    let mut held_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let dx = solver
            .parametric_step(&pins, d)
            .expect("parametric_step ok");
        held_steps.push(dx);
    }
    let held_step_dt = t.elapsed();
    let held_total = solve_dt + held_step_dt;

    // (2) Cold path: N fresh SensSolve runs (one full IPM each).
    let mut cold_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
            eta1: 5.0,
            eta2: 1.0,
        }));
        let mut app = make_app();
        let r = SensSolve::new(pins.clone())
            .with_deltas(d.clone())
            .run(&mut app, tnlp);
        cold_steps.push(r.dx.expect("dx"));
    }
    let cold_total = t.elapsed();

    // Cross-check: held vs cold dx should agree.
    let mut max_err = 0.0_f64;
    for (h, c) in held_steps.iter().zip(cold_steps.iter()) {
        for (a, b) in h.iter().zip(c.iter()) {
            max_err = max_err.max((a - b).abs());
        }
    }

    println!("Workload: 1 IPM solve + {N} parametric steps (Δeta2 sweep).");
    println!();
    println!(
        "Held-factor path : solve {:.2} ms + {N} steps {:.2} ms = {:.2} ms total",
        solve_dt.as_secs_f64() * 1e3,
        held_step_dt.as_secs_f64() * 1e3,
        held_total.as_secs_f64() * 1e3,
    );
    println!(
        "Cold-restart path:                          {N} fresh solves = {:.2} ms total",
        cold_total.as_secs_f64() * 1e3,
    );
    println!();
    println!(
        "Speedup of held-factor over cold-restart: {:.1}x",
        cold_total.as_secs_f64() / held_total.as_secs_f64()
    );
    println!("Numerical agreement (held vs cold dx): max |err| = {max_err:.2e}");
}

examples/rh_orientation_check.rs (line 144)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(PinnedQuadratic { p0: 1.0, p1: 2.0 }));
    let mut app = IpoptApplication::new();
    app.options_mut()
        .set_integer_value("print_level", 0, true, false)
        .unwrap();
    app.options_mut()
        .set_string_value("sb", "yes", true, false)
        .unwrap();
    app.initialize().unwrap();

    let mut solver = Solver::new(app, tnlp);
    let status = solver.solve();
    assert!(
        matches!(
            status,
            ApplicationReturnStatus::SolveSucceeded
                | ApplicationReturnStatus::SolvedToAcceptableLevel
        ),
        "solve failed: {status:?}"
    );

    let (hr, vals, vecs) = solver
        .compute_reduced_hessian_eigen(&[0, 1], 1.0)
        .expect("reduced Hessian");

    println!("model:  min x0² + x1² + x0·x1  s.t.  x0 = 1, x1 = 2");
    println!("H    = [[2, 1], [1, 2]]   (eigenvalues 1 and 3)");
    println!();
    println!("compute_reduced_hessian(pins=[0, 1]):");
    for i in 0..2 {
        println!("  [{:>9.6}, {:>9.6}]", hr[i], hr[i + 2]);
    }
    println!(
        "  eigenvalues (ascending) = [{:>9.6}, {:>9.6}]",
        vals[0], vals[1]
    );
    for j in 0..2 {
        println!(
            "  eigenvector[{j}]          = [{:>9.6}, {:>9.6}]",
            vecs[2 * j],
            vecs[2 * j + 1]
        );
    }
    println!();

    // Three candidates, discriminated by magnitude as well as sign:
    // inv(H) = [[2/3, -1/3], [-1/3, 2/3]].
    let h_inv: [Number; 4] = [2.0 / 3.0, -1.0 / 3.0, -1.0 / 3.0, 2.0 / 3.0];
    let candidates: [(&str, [Number; 4]); 3] = [
        ("+H_R", H),
        ("-H_R", [-H[0], -H[1], -H[2], -H[3]]),
        ("H_R⁻¹", h_inv),
    ];
    for (name, want) in candidates {
        let err = (0..4)
            .map(|k| (hr[k] - want[k]).abs())
            .fold(0.0 as Number, Number::max);
        println!(
            "  vs {name:<6} max|Δ| = {err:.3e}  {}",
            if err < 1e-7 { "← MATCH" } else { "" }
        );
    }
    println!();
    println!(
        "So the ascending spectrum runs STIFFEST first: {:.6} is the curvature-3",
        vals[0]
    );
    println!(
        "mode and {:.6} the curvature-1 (soft) one — the reverse of the",
        vals[1]
    );
    println!("order a caller reading `+H_R` would assume.");
}

pub fn converged(&self) -> Option<Ref<'_, ConvergedState>>

Borrow the converged state, if a successful solve has been run. Returns None if no solve has run or if the most recent solve failed before reaching convergence.

Examples found in repository

examples/sensitivity_session.rs (line 169)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

examples/parametric_mpc.rs (line 186)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    // One full IPM solve at the nominal parameter.
    let mut solver = Solver::new(make_app(), tnlp);
    let t0 = Instant::now();
    let status = solver.solve();
    let solve_dt = t0.elapsed();
    println!(
        "IPM solve at nominal (eta1=5.0, eta2=1.0): status={status:?}, {:.3} ms",
        solve_dt.as_secs_f64() * 1e3
    );
    assert!(solver.converged().is_some());

    // 10 parametric steps as eta2 sweeps from 1.0 to 1.5. Each one is
    // a back-solve against the held factor — no IPM iteration.
    let pins = vec![2 as Index, 3];
    let mut total_step_dt = std::time::Duration::ZERO;
    let mut steps = Vec::new();
    for k in 1..=10 {
        let d_eta2 = 0.05 * k as f64;
        let t = Instant::now();
        let dx = solver
            .parametric_step(&pins, &[0.0, d_eta2])
            .expect("parametric_step ok");
        total_step_dt += t.elapsed();
        steps.push((d_eta2, dx));
    }

    println!("\nparametric steps against the held factor:");
    for (d_eta2, dx) in &steps {
        println!(
            "  Δeta2 = {d_eta2:+.2}  -> Δx_primal = [{:+.5}, {:+.5}, {:+.5}]",
            dx[0], dx[1], dx[2]
        );
    }
    let avg_step_us = total_step_dt.as_secs_f64() * 1e6 / steps.len() as f64;
    println!(
        "\n10 parametric steps: total {:.3} ms, mean {avg_step_us:.1} µs/step",
        total_step_dt.as_secs_f64() * 1e3
    );
    println!(
        "\nRatio: each parametric step is roughly {:.0}x cheaper than the IPM solve.",
        solve_dt.as_secs_f64() / (total_step_dt.as_secs_f64() / steps.len() as f64),
    );
}

pub fn kkt_dim(&self) -> Option<usize>

Total dimension of the compound KKT vector (sum of block_dims). Returns None if no converged factor is held.

Examples found in repository

examples/sensitivity_session.rs (line 188)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

pub fn block_dims(&self) -> Option<[usize; 8]>

Block dimensions of the compound KKT vector in (x, s, y_c, y_d, z_l, z_u, v_l, v_u) order. Returns None if no converged factor is held.

pub fn classify_activity(&self) -> Result<ActivityReport, SolverError>

Classify every bounded variable and every finite-bounded inequality row of the converged solve by activity: see crate::activity and dev-notes/covariance-information-roadmap.md item 0 (gh #362).

Requires the held solve to have run with bound_relax_factor=0 (the Ipopt default is 1e-8): with relaxed bounds the solver’s slacks are measured against perturbed bounds, and the complementarity products the classifier reads no longer track μ.

The guard reads ConvergedState::bound_relax_factor — the value that solve ran under — not the application’s current options. Setting the option after the fact neither unlocks a state whose bounds were relaxed nor invalidates one whose bounds were not; re-solve to change the answer.

Neither classes’ q is a reduced curvature

A variable’s ratio is Σ_i/|H_ii|, and at a kink the multiplier is generated by the curvature reduced along that coordinate, not by the diagonal. The two agree only where the coordinate is decoupled, so a genuine kink coupled to a neighbour reads AMBIGUOUS here at any tolerance (gh#763). Do not read that class as “probably not a kink”: use Self::reduced_activity, which normalizes by the reduced curvature at one back-solve per coordinate.

A row’s ratio divides by the curvature along the row’s own gradient instead, which is a genuine directional curvature but still not a reduced one, so the same warning and the same remedy apply there: its ratio is reduced/directional and Self::reduced_row_activity answers the kink question (gh#804).

pub fn reduced_activity( &self, user_vars: &[usize], ) -> Result<ReducedActivityReport, SolverError>

Self::classify_activity’s per-variable verdict for user_vars, re-measured against the reduced curvature along each coordinate instead of the Hessian diagonal — one back-solve against the held factor per variable (gh#763).

classify_activity normalizes a variable’s Σ by H_ii, but the multiplier at a kink is generated by the curvature left after the other free variables re-optimize. The two agree only where the coordinate is decoupled, so a genuine kink coupled to a neighbour reads AMBIGUOUS there at any tolerance — the ratio is μ-independent. Ask here and the same kink reads WEAKLY_ACTIVE.

Indices are user space (full-x), as the report’s are. The intended call is over a report’s ambiguous entries:

let report = solver.classify_activity()?;
let ask: Vec<usize> = (0..report.var_status.len())
    .filter(|&i| report.var_status[i] == AMBIGUOUS)
    .collect();
let refined = solver.reduced_activity(&ask)?;

The cost is one back-solve per index, so it is a refinement to call over the entries in question, not over every bounded variable of a large model. See [crate::activity::reduced_activity] for the algebra and the edge cases.

Requires the held solve to have run with bound_relax_factor=0 for the same reason Self::classify_activity does: both read the slacks relaxed bounds shift.

pub fn reduced_row_activity( &self, user_rows: &[usize], ) -> Result<ReducedRowActivityReport, SolverError>

Self::classify_activity’s per-ROW verdict for user_rows, re-measured against the reduced curvature along each row’s gradient instead of the directional curvature ∇dᵀH∇d/‖∇d‖² — one back-solve against the held factor per row (gh#804).

The row counterpart of Self::reduced_activity, and the same defect one block over. A row’s directional denominator is a genuine curvature along the row’s own gradient — strictly better than the variable path’s bare H_ii, which is why gh#763 fixed the variables first — but it is still not reduced: it does not account for the other free coordinates re-optimizing, and the multiplier is generated by what is left after they do. So a row’s ratio there is reduced/directional, equal to 1 only where the row’s direction is decoupled from the remaining free space, and a genuine row kink that is coupled reads AMBIGUOUS at any tolerance — the ratio is μ-independent. Ask here and the same kink reads WEAKLY_ACTIVE.

Indices are user space (full-g), as the report’s are — equality rows included, which report EQUALITY rather than being an error. The intended call is over a report’s ambiguous rows:

let report = solver.classify_activity()?;
let ask: Vec<usize> = (0..report.row_status.len())
    .filter(|&j| report.row_status[j] == AMBIGUOUS)
    .collect();
let refined = solver.reduced_row_activity(&ask)?;

The cost is one back-solve per index, so it is a refinement to call over the rows in question, not over every bounded row of a large model. See [crate::activity::reduced_row_activity] for the algebra and the edge cases.

Requires the held solve to have run with bound_relax_factor=0 for the same reason Self::classify_activity does: both read the slacks relaxed bounds shift.

pub fn row_normal(&self, user_row: usize) -> Result<Vec<Number>, SolverError>

The gradient of user constraint row user_row at the converged iterate, in user variable order (length n_full_x) and in natural (unscaled) units: the internal Jacobian row carries the solver’s per-row c_scale/d_scale, which is divided out here, so this is the gradient of the row as the user wrote it. Equality and inequality rows alike; entries for fixed (make_parameter-removed) variables are 0 because the solve dropped their columns. Errors on an out-of-range row.

Serves the covariance roadmap’s item 1: a binding row’s normal restricted to the fitted block is the projection direction.

pub fn hessian_vec(&self, v: &[Number]) -> Result<Vec<Number>, SolverError>

The exact Lagrangian Hessian times a user-space vector, in user variable order and natural units (see [crate::activity::hessian_vec]). Errors on a length mismatch.

pub fn kkt_solve( &self, rhs: &[Number], lhs: &mut [Number], ) -> Result<(), SolverError>

Solve K · lhs = rhs against the converged KKT factor. Both slices must have length kkt_dim(); the layout is the flat x || s || y_c || y_d || z_l || z_u || v_l || v_u packing.

K here is the natural-units (unscaled) KKT matrix: when the IPM solved with active NLP scaling, the backsolver scales the RHS/solution (all eight blocks, including the z/v bound-multiplier rows) so callers pass and receive data in the user’s own units (pounce#128) — see crate::PdSensBacksolver::solve. For the raw scaled-space back-solve use Self::kkt_solve_scaled.

Examples found in repository

examples/sensitivity_session.rs (line 191)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

pub fn kkt_solve_scaled( &self, rhs: &[Number], lhs: &mut [Number], ) -> Result<(), SolverError>

Self::kkt_solve without the natural-units conjugation: the back-solve runs against the factor exactly as the IPM holds it (the solver’s internal scaled space). Identical to kkt_solve when no NLP scaling is active. “Scaled space” includes a user-scaling change of variables (gh#486), so on such a solve the x and z blocks here are in the substituted coordinates x̃ = d ⊙ x, not the model’s.

pub fn kkt_solve_many( &self, rhs_flat: &[Number], lhs_flat: &mut [Number], n_rhs: usize, ) -> Result<(), SolverError>

Batched-RHS back-solve. rhs_flat and lhs_flat are row-major (n_rhs, kkt_dim) buffers; each row is solved against the same converged factor. Equivalent in result to looping Self::kkt_solve but reuses one IteratesVector for the RHS and one for the result across all n_rhs calls — see crate::algorithm_backsolver::PdSensBacksolver::solve_many.

pub fn kkt_solve_many_scaled( &self, rhs_flat: &[Number], lhs_flat: &mut [Number], n_rhs: usize, ) -> Result<(), SolverError>

Self::kkt_solve_many without the natural-units conjugation (the batched sibling of Self::kkt_solve_scaled).

pub fn parametric_step( &self, pin_constraint_indices: &[Index], deltas: &[Number], ) -> Result<Vec<Number>, SolverError>

First-order parametric step Δx ≈ ∂x*/∂p · Δp for a set of pinned equality constraints. pin_constraint_indices are 0-based indices into the user’s g(x); deltas is the perturbation Δp (same length).

Returns the n_x-long primal step. For the full KKT-space step, use Self::kkt_solve directly.

Examples found in repository

examples/sensitivity_session.rs (line 176)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

examples/parametric_mpc.rs (line 197)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    // One full IPM solve at the nominal parameter.
    let mut solver = Solver::new(make_app(), tnlp);
    let t0 = Instant::now();
    let status = solver.solve();
    let solve_dt = t0.elapsed();
    println!(
        "IPM solve at nominal (eta1=5.0, eta2=1.0): status={status:?}, {:.3} ms",
        solve_dt.as_secs_f64() * 1e3
    );
    assert!(solver.converged().is_some());

    // 10 parametric steps as eta2 sweeps from 1.0 to 1.5. Each one is
    // a back-solve against the held factor — no IPM iteration.
    let pins = vec![2 as Index, 3];
    let mut total_step_dt = std::time::Duration::ZERO;
    let mut steps = Vec::new();
    for k in 1..=10 {
        let d_eta2 = 0.05 * k as f64;
        let t = Instant::now();
        let dx = solver
            .parametric_step(&pins, &[0.0, d_eta2])
            .expect("parametric_step ok");
        total_step_dt += t.elapsed();
        steps.push((d_eta2, dx));
    }

    println!("\nparametric steps against the held factor:");
    for (d_eta2, dx) in &steps {
        println!(
            "  Δeta2 = {d_eta2:+.2}  -> Δx_primal = [{:+.5}, {:+.5}, {:+.5}]",
            dx[0], dx[1], dx[2]
        );
    }
    let avg_step_us = total_step_dt.as_secs_f64() * 1e6 / steps.len() as f64;
    println!(
        "\n10 parametric steps: total {:.3} ms, mean {avg_step_us:.1} µs/step",
        total_step_dt.as_secs_f64() * 1e3
    );
    println!(
        "\nRatio: each parametric step is roughly {:.0}x cheaper than the IPM solve.",
        solve_dt.as_secs_f64() / (total_step_dt.as_secs_f64() / steps.len() as f64),
    );
}

examples/sensitivity_factor_reuse_bench.rs (line 174)

fn main() {
    let deltas: Vec<Vec<Number>> = (1..=N).map(|k| vec![0.0, 0.01 * k as Number]).collect();
    let pins = vec![2 as Index, 3];

    // (1) Held-factor path: 1 solve + N parametric steps.
    let t = Instant::now();
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));
    let mut solver = Solver::new(make_app(), tnlp);
    solver.solve();
    let solve_dt = t.elapsed();
    let mut held_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let dx = solver
            .parametric_step(&pins, d)
            .expect("parametric_step ok");
        held_steps.push(dx);
    }
    let held_step_dt = t.elapsed();
    let held_total = solve_dt + held_step_dt;

    // (2) Cold path: N fresh SensSolve runs (one full IPM each).
    let mut cold_steps = Vec::with_capacity(N);
    let t = Instant::now();
    for d in &deltas {
        let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
            eta1: 5.0,
            eta2: 1.0,
        }));
        let mut app = make_app();
        let r = SensSolve::new(pins.clone())
            .with_deltas(d.clone())
            .run(&mut app, tnlp);
        cold_steps.push(r.dx.expect("dx"));
    }
    let cold_total = t.elapsed();

    // Cross-check: held vs cold dx should agree.
    let mut max_err = 0.0_f64;
    for (h, c) in held_steps.iter().zip(cold_steps.iter()) {
        for (a, b) in h.iter().zip(c.iter()) {
            max_err = max_err.max((a - b).abs());
        }
    }

    println!("Workload: 1 IPM solve + {N} parametric steps (Δeta2 sweep).");
    println!();
    println!(
        "Held-factor path : solve {:.2} ms + {N} steps {:.2} ms = {:.2} ms total",
        solve_dt.as_secs_f64() * 1e3,
        held_step_dt.as_secs_f64() * 1e3,
        held_total.as_secs_f64() * 1e3,
    );
    println!(
        "Cold-restart path:                          {N} fresh solves = {:.2} ms total",
        cold_total.as_secs_f64() * 1e3,
    );
    println!();
    println!(
        "Speedup of held-factor over cold-restart: {:.1}x",
        cold_total.as_secs_f64() / held_total.as_secs_f64()
    );
    println!("Numerical agreement (held vs cold dx): max |err| = {max_err:.2e}");
}

pub fn parametric_step_bounded( &self, pin_constraint_indices: &[Index], deltas: &[Number], max_iter: usize, bound_eps: Option<Number>, ) -> Result<(Vec<Number>, Vec<Index>, RefineStop), SolverError>

Parametric step with the bounds respected by pinning, not by clamping. Returns the n_x-long primal step, the rows it constrained to reach it, and why the refinement stopped.

Self::parametric_step answers where the linear predictor points, which can be outside the box. Clamping a coordinate back to its bound leaves every other coordinate at its predictor value, so the answer is feasible but no longer consistent with the KKT relations. This instead adds a row pinning each offending coordinate at its bound and re-solves, so the others move to stay consistent under the pins, which is the refinement upstream runs under sens_boundcheck.

A pass takes every crossing it can see, pins them together, and re-solves, so the loop ends when nothing is left outside rather than when the passes run out. Each pass rebuilds the Schur complement over the pins so far, so a pass carrying k of them costs one dense k × k solve and k + 1 back-solves; the factorization itself is never rebuilt for a pin.

What counts as outside a bound is the eps argument when the caller passes one, and the solve’s own margin when it passes None: the solve was willing to leave a converged point bound_relax_factor outside its bound, so anything within that is on the bound. An unrelaxed solve gets a roundoff floor.

Passes stop when nothing is outside its bound by that much, when a pin cannot be achieved because the pins have exhausted the problem’s degrees of freedom, or at max_iter, which is a safety limit rather than a budget: it took one pin per pass until gh#732, where a model with more crossings than passes had its answer picked by the limit. None of those is an error, and the returned crate::boundcheck::RefineStop says which happened.

pub fn parametric_step_path( &self, pin_constraint_indices: &[Index], deltas: &[Number], max_iter: usize, ) -> Result<(Vec<Number>, Vec<PathSegment>), SolverError>

Parametric step applied a little at a time instead of taken whole, stopping wherever the active set changes and continuing from there under the new one. Returns the primal step and the breakpoints crossed.

Self::parametric_step_bounded decides every condition at the base point, which is upstream’s fix-relax. This is past it: the result is piecewise linear in the parameter, exact for a QP because a QP’s solution is piecewise affine in the parameter, and still a predictor for an NLP because nothing is re-linearized between breakpoints.

max_iter caps the breakpoints crossed. It is in practice a budget on factorizations, since a pin is a back-solve against the held factor while a release re-factors.

pub fn parametric_step_path_decided( &self, pin_constraint_indices: &[Index], deltas: &[Number], max_iter: usize, held_var_rows: &[Index], ) -> Result<(Vec<Number>, Vec<PathSegment>), SolverError>

Self::parametric_step_path with the weak-row decision supplied by the caller instead of searched for. held_var_rows names the var-x rows of the weakly active bounds the direction holds; every other weakly active bound is forced into the walk’s base-activity table as a leaving row. A row left there is still reachable, so a caller that hands in an empty held list — every weak row declared a leaver, which is what an undecided study of the all-released step does — gets the bound back at the fraction the walk finds the direction pressing into it, rather than an answer outside the box (gh#852). Study surface for an externally solved eq. 14 QP.

pub fn correct_step( &self, pin_constraint_indices: &[Index], deltas: &[Number], step: &[Number], max_iter: usize, ) -> Result<(Vec<Number>, CorrectorReport), SolverError>

Newton iterations on the barrier system, refining a step that some mode already produced.

step is a full compound step, the shape Self::parametric_step_full returns, so any mode’s result can be handed in. Every correction pays one derivative evaluation and one factorization at the predicted point, and each iteration after that costs one back-solve. Returns the refined step and a [CorrectorReport] saying what the iterations bought.

The corrector aims at the barrier solution at the μ the solve finished on, not at a re-solve, so the accuracy it can reach is bounded by that offset. Its operator is assembled at the PREDICTED point, every block: the Hessian, the constraint Jacobians, and the barrier diagonal all evaluated at the stepped iterate with the step’s own multipliers, and the predictor’s active set applied to the diagonal in that frame. A base solve the sigma ceiling (gh#737) touched, or one that crossed over into the declared frame (gh#654), is no exception: both rules are re-derived at the predicted point rather than read from the base-point diagonals stored for the held factor’s own back-solves. A chord iteration contracts at the rate the distance between its operator and the true Jacobian sets, and the predicted point is where the truth is. Under a limited-memory solve the quasi-Newton matrix is kept as is, since no exact Hessian exists to evaluate elsewhere. Where the perturbation needs a bound to leave the active set that the step’s endpoint does not show, no released row is applied: the step’s clamped multiplier leaves a weak diagonal entry there, the iterations can move the coordinate partway off the bound, and the answer is not the re-solve. The release-deciding modes are the ones that cross exactly. CorrectorReport::improved reports whether the residual fell; when it did not, the step handed back is the caller’s own.

The returned point always satisfies the variable bounds, since the barrier residual is undefined outside them and the fraction-to-boundary rule keeps every iterate inside. A step that arrives pointing out of the box is therefore put back in before the first iteration, which means max_iter = 0 is not a no-op: it costs the derivative evaluation and the residual evaluation, no back-solve, and reports the residual the caller’s step leaves.

Errors with SolverError::SensComputationFailed when the barrier residual at that starting point is not finite, which is what a predicted point outside the domain of one of the model’s functions gives (gh#845). A declared bound is protection here, since the clamp above puts the coordinate back inside it; a variable held in a function’s domain by a constraint has no bound to be put back inside, and an ordinary log, sqrt or reciprocal is then reachable by a large enough perturbation. There is no correction to make from such a point, so it is an error rather than a report – and never a step full of NaN carrying residual = 0.0 and converged = true.

pub fn parametric_step_bounded_decided( &self, pin_constraint_indices: &[Index], deltas: &[Number], max_iter: usize, held_var_rows: &[Index], bound_eps: Option<Number>, ) -> Result<(Vec<Number>, Vec<Index>, RefineStop), SolverError>

Self::parametric_step_bounded with the weak-row decision supplied by the caller instead of searched for. The direction is computed for the given working set (all weak rows released, the held variables pinned through Schur rows), then refined onto the bounds exactly as the searched variant does. Study surface for an externally solved eq. 14 QP.

pub fn path_direction_decided( &self, pin_constraint_indices: &[Index], deltas: &[Number], released_bound_rows: &[Index], held_primal_rows: &[Index], operator: PathOperator, ) -> Result<(Vec<Number>, Vec<Number>), SolverError>

crate::boundcheck::path_direction for a working set the caller names, and the force each held row carries under it. Study surface, and the seam the two pins are measured against each other through.

released_bound_rows are compound bound-multiplier rows, as Self::weakly_active_bounds reports them; held_primal_rows are primal rows – x block for a variable, s block for a constraint’s own limit (gh#928).

operator picks which operator the walk’s (single, exact) Schur pin is applied to. Plain is the released system as the factorization already holds it – one factorization for the whole segment, and no inverse at all when releasing two curvature-free variables that share a row leaves the stationarity rows dependent (gh#930). Regularized raises the pinned diagonals until it is invertible, at the cost of rebuilding the diagonal per solve, so the factorization cache misses. Preferred is what the walk itself runs: the plain one, falling back when it fails or when its pins do not take, which is not the same test – see crate::boundcheck::path_direction.

Both operators give the same answer: the Schur row enforces Eᵀ w = 0, which annihilates the added diagonal, so the system solved is the released one either way and both return values agree in value, frame and units. issue_930_two_curvature_free_releases.rs measures that.

pub fn parametric_step_release_all( &self, pin_constraint_indices: &[Index], deltas: &[Number], ) -> Result<(Vec<Number>, usize), SolverError>

The all-released step: the plain parametric step solved with every weakly active bound’s row released, and nothing decided.

This is Self::parametric_step_directional’s first back-solve returned as the answer instead of refined. The caller trades the engagement’s budget for whatever violations the released direction carries at weak bounds the perturbation actually holds, which come back as crossings for the mode’s clamp, pins, or path segments, or for a correction, to handle. A clean base point takes the plain step. Returns the direction over the model’s variables and the number of rows released.

pub fn parametric_step_directional( &self, pin_constraint_indices: &[Index], deltas: &[Number], max_iter: usize, ) -> Result<(Vec<Number>, Vec<usize>, usize), SolverError>

The eq. 14 directional derivative, decided by pounce-qp over the weak rows the direction engages.

One released factorization serves the whole decision: the released Σ is built once and every solve passes the same object, so the factorization cache reuses the factor across the all-released direction and the basis columns. The decision itself is the dual of eq. 14 restricted to the weak rows the direction engages: with a_k the signed unit vector of weak row k (positive for a lower bound), X_k = K_rel^{-1} a_k, S = aᵀX and m = aᵀd0, the pin forces λ solve

    min  ½ λᵀ S λ + mᵀ λ    s.t.  λ ≥ 0

whose KKT conditions are eq. 14’s complementarity: a released row moves to its feasible side (the QP gradient Sλ + m ≥ 0) and a held row’s pin force is nonnegative. Rows outside the engaged set are verified against the decided direction and the set expands until no new row violates. Nothing reads the perturbation’s size, so the decision is linear in the step.

An engaged row is decided only when its bound is at a kink, read off kappa = sigma * S_kk, the barrier weight times the row’s own diagonal of the reduced matrix. sigma equals the curvature reduced along the coordinate at an exact kink, and S_kk is that reduced curvature’s inverse, so kappa is 1 there at any curvature, coupling, or scaling, and it falls as the squared ratio of kink width to slack away from one. A row below KAPPA_MIN is dropped from the engaged set and its plain movement stands: its bound is too far from a kink for a pin force to decide, and the error of leaving it undecided is bounded by its own slack, order sqrt(mu) at the threshold. A coordinate an equality pins is the limiting case, S_kk exactly zero, dropped by the same test.

max_iter is the total back-solve budget: the all-released solve, every basis column, and the combined solve that recovers the direction all count against it. A budget of zero errs before any work. Any budget above that pays the all-released factorization first, because which rows engage is only known once that solve has run, and the shortfall is reported when the basis columns cannot fit. Either way the caller falls back to the one-sided step. Returns the direction, the var-x rows held, and the back-solves spent.

pub fn weakly_active_bounds(&self) -> Result<Vec<WeakBound>, SolverError>

The bounds the activity classifier could not call at the base point: on the bound with a multiplier of the same order as the slack. Each entry is a bound row present in the held factorization, with the side taken from the smaller slack, which is the only side an ambiguous label can come from.

Both WEAKLY_ACTIVE and AMBIGUOUS count as weak here, deliberately: the ambiguous class contains genuine kinks whose coordinate is coupled to a neighbour (gh#763), so treating it as “not a kink” would drop real weak rows. That is why the mislabeling is not a wrong answer in the step path — see Self::reduced_activity for the class itself.

The classifier reports per user variable, in full-x, while the bound context and the factor’s rows are var-x, and the two index spaces diverge from the first fixed variable on. Each var-x row’s status is read through the same map the classifier scattered through, so a fixed variable shifts nothing. Using the full-x index as a factor row instead returns a NEIGHBORING variable’s answer, plausible and wrong, which is the gh#450 hazard the primal_row discipline exists to prevent.

pub fn parametric_step_full( &self, pin_constraint_indices: &[Index], deltas: &[Number], ) -> Result<Vec<Number>, SolverError>

Full KKT-space parametric step for a set of pinned equality constraints: the same computation as Self::parametric_step, returned WITHOUT truncating to the primal block. The layout is the compound KKT vector (x, s, y_c, y_d, z_l, z_u, v_l, v_u); use Self::block_dims for the block sizes and Self::g_multiplier_rows to locate a constraint’s multiplier row. This exposes the multiplier sensitivities ∂λ*/∂p alongside the primal step.

pub fn g_multiplier_rows( &self, g_indices: &[Index], ) -> Result<Vec<Option<Index>>, SolverError>

Flat rows of the compound KKT vector holding the equality multipliers y_c for the given 0-based full-g constraint indices. None for inequalities — their multipliers live in the y_d block, which Self::d_multiplier_rows addresses. Row r of a Self::parametric_step_full result is then ∂λ_g/∂p · Δp.

pub fn d_multiplier_rows( &self, g_indices: &[Index], ) -> Result<Vec<Option<Index>>, SolverError>

Flat rows of the compound KKT vector holding the inequality multipliers y_d for the given 0-based full-g constraint indices. None for equalities (those are Self::g_multiplier_rows’s). The y_d counterpart of that accessor, added by gh#910: parametric_step_full already returned the y_d block, and this is the map that says which row of it belongs to which user constraint.

Reading the row is not the same as the row being a derivative. y_d holds a number for every inequality, in all three activity regimes, and only one of them has a two-sided ∂λ/∂p at all:

• strictly active (s ≈ 0, λ > 0): the row behaves as an equality over a neighbourhood of the solved point, and this row is the same back-solve an equality gets. Well defined.
• inactive (s > 0, λ ≈ 0): the derivative is a structural zero over a neighbourhood; the KKT row carries the barrier’s residue rather than a derivative of anything.
• weakly active (a kink: s ≈ 0 and λ ≈ 0): the two one-sided derivatives differ and no two-sided value exists. The entry holds whichever side the factorization landed on — the silently-wrong-while-reporting-success class. Self::parametric_step_directional is what answers a kink, and it needs a direction.

So a caller reading these rows as ∂λ/∂p must gate on the regime first, and the classifier that answers it is Self::reduced_row_activity, not Self::classify_activity: a genuine kink whose row couples to the remaining free space reports INACTIVE on the directional normalizer at strong enough coupling (gh#804), and INACTIVE is the one class whose derivative a caller may legitimately read as a structural zero. Gating on the cheap classifier would therefore answer “it does not move” about a kink — a wrong answer wearing a refusal’s clothes. That inference, reading an activity class as a proxy for kink-ness, is what shipped gh#756.

pub fn d_slack_rows( &self, g_indices: &[Index], ) -> Result<Vec<Option<Index>>, SolverError>

Flat rows of the compound KKT vector holding an inequality’s slack s, for the given 0-based full-g row indices; None for a row that is not an inequality.

The primal counterpart of Self::d_multiplier_rows, and the discriminator a consumer of Self::parametric_step_path needs. A limit written as g(x) <= cap bounds this slack rather than any variable, so a breakpoint on it carries a primal KKT row in the s block (gh#928). Reading such a row as a var-x index returns a neighbouring variable’s answer, the gh#450 hazard, so a caller that maps rows back to model objects resolves them here.

The s block sits immediately after x, and is indexed by d-block position exactly as y_d is, so this row and Self::d_multiplier_rows’s row name the same inequality from the two sides of its complementarity pair.

pub fn x_primal_rows( &self, x_indices: &[Index], ) -> Result<Vec<Option<Index>>, SolverError>

Flat rows of the compound KKT vector holding the primal values x for the given 0-based full-x variable indices. None where the solve removed the column (x_l == x_u under fixed_variable_treatment = make_parameter), which has no row in the factor at all.

The x counterpart of Self::g_multiplier_rows, and needed for the same reason: a caller holding user-space indices — from the .col file, from Self::classify_activity, from Self::row_normal — cannot index the factor with them directly. Row r of a Self::parametric_step_full result is then ∂x/∂p · Δp for that variable, and e_r is the unit vector selecting its column in a Self::kkt_solve.

pub fn n_full_x(&self) -> Result<usize, SolverError>

The user TNLP’s variable count: the length of a full-x report and the domain of Self::x_primal_rows.

pub fn n_full_g(&self) -> Result<usize, SolverError>

The user TNLP’s constraint count: the length of a full-g report and the domain of Self::reduced_row_activity.

pub fn compute_reduced_hessian( &self, pin_constraint_indices: &[Index], obj_scal: Number, ) -> Result<Vec<Number>, SolverError>

Reduced Hessian over the pinned equality-constraint rows: obj_scal · B K⁻¹ Bᵀ, where B selects the pin_constraint_indices rows of the y_c block and K is the natural-units (unscaled) KKT matrix — active NLP scaling is undone by the backsolver, so −inv of the returned matrix is directly the parameter covariance regardless of nlp_scaling_method (pounce#128). obj_scal survives as a plain extra multiplier (default 1.0); it is no longer needed to recover natural units. Returns the n²-long column-major dense matrix (n = pin_constraint_indices.len()).

Sign convention: this returns −H_R, not H_R (gh#937)

The matrix is the negated reduced Hessian. On a model whose objective Hessian is [[2, 1], [1, 2]] with both variables pinned, this returns [[−2, −1], [−1, −2]]. So a well-posed minimum reports an all-negative spectrum; that is the convention, not an indefiniteness or convergence bug.

The minus is the augmented system’s, and it is why the covariance recipe above negates: pin indices map to the y_c multiplier block, and for K = [[H, Aᵀ], [A, 0]] the (y_c, y_c) block of K⁻¹ is −(A H⁻¹ Aᵀ)⁻¹ — so over pin rows B K⁻¹ Bᵀ is the multiplier sensitivity ∂λ/∂p = −∂²f*/∂p², i.e. ±H_R itself and not a submatrix of an inverse. (The x block of K⁻¹ is an inverse. The two blocks sit on opposite sides of one inversion, which is what makes the CLI’s red_hessian suffix path — upstream sIPOPT’s, selecting x rows — a different quantity rather than the same one with a different sign.)

Negate to read curvature: −hr is H_R, and −inv(hr) is the covariance. Pinned by tests/issue_937_reduced_hessian_sign.rs; demonstrated by examples/rh_orientation_check.rs.

Equivalent to crate::SensSolve::with_reduced_hessian but usable post-hoc on a held Solver. For the solver-space (pre-#128) value use Self::compute_reduced_hessian_scaled; the factors themselves are exposed via Self::nlp_scaling / Self::pin_g_scaling.

Examples found in repository

examples/sensitivity_session.rs (line 183)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(ParametricTNLP {
        eta1: 5.0,
        eta2: 1.0,
    }));

    let mut solver = Solver::new(make_app(), tnlp);
    let status = solver.solve();
    println!("solve status: {status:?}");
    assert!(solver.converged().is_some(), "solver did not converge");

    let pins = vec![2 as Index, 3];

    // Two cheap parametric steps against the same factor.
    for deltas in &[vec![-0.5, 0.0], vec![0.0, 0.2]] {
        let dx = solver
            .parametric_step(&pins, deltas)
            .expect("parametric_step ok");
        println!("parametric_step(deltas={deltas:?}) -> dx = {dx:?}");
    }

    // Reduced Hessian over the same pinned-row set.
    let hr = solver
        .compute_reduced_hessian(&pins, 1.0)
        .expect("reduced Hessian ok");
    println!("reduced Hessian (2x2, column-major) = {hr:?}");

    // Raw back-solve against a zero RHS — must come back zero.
    let dim = solver.kkt_dim().expect("kkt_dim available");
    let rhs = vec![0.0; dim];
    let mut lhs = vec![1.0; dim];
    solver.kkt_solve(&rhs, &mut lhs).expect("kkt_solve ok");
    let max_abs = lhs.iter().fold(0.0_f64, |a, b| a.max(b.abs()));
    println!("kkt_solve(0) max |lhs| = {max_abs:e}");
    assert!(max_abs < 1e-10);
}

pub fn compute_reduced_hessian_eigen( &self, pin_constraint_indices: &[Index], obj_scal: Number, ) -> Result<(Vec<Number>, Vec<Number>, Vec<Number>), SolverError>

Self::compute_reduced_hessian plus its eigendecomposition — (H_R, eigenvalues, eigenvectors).

The curvature on the null space of the active constraints is the question; its spectrum is what answers “is this parameter identifiable, and along which direction”. SensSolve has offered that since gh#561 (crate::SensSolve::with_reduced_hessian_eigen), and the session API did not — so a caller holding a Solver had to re-solve the whole NLP through the one-shot builder to get a decomposition of a matrix it already had. That is the gap this closes; the numbers are the one-shot path’s, from the same pounce_linalg::symmetric_eigen.

Eigenvectors are column-major, length n², column j belonging to eigenvalue j, and sign-pinned by symmetric_eigen so a column read as a direction reproduces across builds.

This is the spectrum of −H_R, so ascending runs stiffest first

Self::compute_reduced_hessian returns the negated reduced Hessian (gh#937, and see its docs for why). The eigenvalues here are that matrix’s, in ascending order — which on −H_R runs from most negative to least, i.e. stiffest mode first and softest last, the reverse of what the identifiability reading wants. On H = [[2, 1], [1, 2]] fully pinned they come back [−3, −1]: the leading column is the curvature-3 stiff direction, the trailing one the curvature-1 soft direction.

So a caller taking the leading columns as the least-identifiable directions gets the best-identified ones, and nothing looks wrong — the vectors are unit-norm, sign-pinned and entirely plausible. Either negate the eigenvalues and reverse the order, or read the trailing columns as the soft modes. Pinned by tests/issue_937_reduced_hessian_sign.rs.

Examples found in repository

examples/rh_orientation_check.rs (line 155)

fn main() {
    let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(PinnedQuadratic { p0: 1.0, p1: 2.0 }));
    let mut app = IpoptApplication::new();
    app.options_mut()
        .set_integer_value("print_level", 0, true, false)
        .unwrap();
    app.options_mut()
        .set_string_value("sb", "yes", true, false)
        .unwrap();
    app.initialize().unwrap();

    let mut solver = Solver::new(app, tnlp);
    let status = solver.solve();
    assert!(
        matches!(
            status,
            ApplicationReturnStatus::SolveSucceeded
                | ApplicationReturnStatus::SolvedToAcceptableLevel
        ),
        "solve failed: {status:?}"
    );

    let (hr, vals, vecs) = solver
        .compute_reduced_hessian_eigen(&[0, 1], 1.0)
        .expect("reduced Hessian");

    println!("model:  min x0² + x1² + x0·x1  s.t.  x0 = 1, x1 = 2");
    println!("H    = [[2, 1], [1, 2]]   (eigenvalues 1 and 3)");
    println!();
    println!("compute_reduced_hessian(pins=[0, 1]):");
    for i in 0..2 {
        println!("  [{:>9.6}, {:>9.6}]", hr[i], hr[i + 2]);
    }
    println!(
        "  eigenvalues (ascending) = [{:>9.6}, {:>9.6}]",
        vals[0], vals[1]
    );
    for j in 0..2 {
        println!(
            "  eigenvector[{j}]          = [{:>9.6}, {:>9.6}]",
            vecs[2 * j],
            vecs[2 * j + 1]
        );
    }
    println!();

    // Three candidates, discriminated by magnitude as well as sign:
    // inv(H) = [[2/3, -1/3], [-1/3, 2/3]].
    let h_inv: [Number; 4] = [2.0 / 3.0, -1.0 / 3.0, -1.0 / 3.0, 2.0 / 3.0];
    let candidates: [(&str, [Number; 4]); 3] = [
        ("+H_R", H),
        ("-H_R", [-H[0], -H[1], -H[2], -H[3]]),
        ("H_R⁻¹", h_inv),
    ];
    for (name, want) in candidates {
        let err = (0..4)
            .map(|k| (hr[k] - want[k]).abs())
            .fold(0.0 as Number, Number::max);
        println!(
            "  vs {name:<6} max|Δ| = {err:.3e}  {}",
            if err < 1e-7 { "← MATCH" } else { "" }
        );
    }
    println!();
    println!(
        "So the ascending spectrum runs STIFFEST first: {:.6} is the curvature-3",
        vals[0]
    );
    println!(
        "mode and {:.6} the curvature-1 (soft) one — the reverse of the",
        vals[1]
    );
    println!("order a caller reading `+H_R` would assume.");
}

pub fn compute_reduced_hessian_scaled( &self, pin_constraint_indices: &[Index], obj_scal: Number, ) -> Result<Vec<Number>, SolverError>

The reduced Hessian as the solver’s internal scaled space sees it — the value Self::compute_reduced_hessian returned before pounce#128: H̃_ij = (df / (dc_i·dc_j)) · H_ij. Identical to compute_reduced_hessian when no NLP scaling is active.

Sign: this is Self::compute_reduced_hessian’s −H_R multiplied through by df / (dc_i·dc_j) (gh#937), so unlike the natural-units value its orientation is not fixed. Measured on a fully pinned [[2, 1], [1, 2]]: [[−2, −1], [−1, −2]] by default, but [[2, 1], [1, 2]] under obj_scaling_factor = −1, where df carries the minus that makes a maximization a minimization. Read the sign off the reported factors (Self::nlp_scaling, Self::pin_g_scaling) rather than assuming it, or use the natural-units value, whose −H_R holds whatever the scaling.

pub fn nlp_scaling( &self, ) -> Result<(Number, Option<Vec<Number>>, Option<Vec<Number>>), SolverError>

Effective NLP scaling the IPM applied on the most recent converged solve: (obj_scaling_factor, c_scale, d_scale). (1.0, None, None) ⇔ no scaling was active. The vectors are per-row factors over the algorithm’s equality (c) and inequality (d) blocks.

pub fn variable_scaling(&self) -> Result<Option<Vec<Number>>, SolverError>

The per-variable user-scaling factors d the held solve ran under (gh#486), in the user TNLP’s full-x space, or None when the solve applied no change of variables.

Every accessor on this type already reports natural units, so this is diagnostic rather than a correction a caller has to apply — it answers “was this solve conditioned, and by how much”, the x-axis counterpart of Self::nlp_scaling.

pub fn kkt_perturbations(&self) -> Result<[Number; 4], SolverError>

Inertia-correction perturbations (δ_x, δ_s, δ_c, δ_d) baked into the held KKT factor. All zero ⇔ the final factorization was unregularized and the natural-units back-solves invert the exact KKT matrix — see crate::PdSensBacksolver::kkt_perturbations.

pub fn pin_g_scaling( &self, pin_constraint_indices: &[Index], ) -> Result<Vec<Number>, SolverError>

Per-pin equality-row scaling factors dc_i (1.0 entries when no constraint scaling is active), ordered like pin_constraint_indices.