1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
//! gh#605 — the least-square primal initializer must not hand the
//! algorithm a *worse* starting point than the user supplied.
//!
//! `least_square_init_primal` replaces `x0` with the minimum-norm
//! solution of the **linearized** constraints. Where the Jacobian is
//! small relative to the residual, that linearization asks for a huge
//! correction and the true nonlinear violation at the far end is far
//! worse than where it started.
//!
//! This model is the smallest case that shows it:
//!
//! ```text
//! min (x0 - 3)^2 + (x1 - 3)^2
//! s.t. x0^2 + x1^2 = 1 started at (0.05, 0.05)
//! ```
//!
//! At `x0 = (0.05, 0.05)` the constraint residual is `-0.995` and the
//! Jacobian is `(0.1, 0.1)`. The min-norm linearized correction is
//! `d = (4.975, 4.975)`, landing at `(5.025, 5.025)` where the true
//! violation is `2*5.025^2 - 1 ≈ 49.5` — **50x worse** than the
//! `0.995` it started with.
//!
//! The safeguard added in gh#605 scores every trial step on the true
//! nonlinear violation and backtracks (or declines outright) when it
//! does not improve, so the algorithm never starts from the far point.
use pounce_algorithm::application::IpoptApplication;
use pounce_common::types::Number;
use pounce_nlp::tnlp::{
BoundsInfo, IndexStyle, IpoptCq, IpoptData, NlpInfo, Solution, SparsityRequest, StartingPoint,
TNLP,
};
use std::cell::RefCell;
use std::rc::Rc;
#[derive(Default)]
struct PoorLinearization;
impl TNLP for PoorLinearization {
fn get_nlp_info(&mut self) -> Option<NlpInfo> {
Some(NlpInfo {
n: 2,
m: 1,
nnz_jac_g: 2,
nnz_h_lag: 2,
index_style: IndexStyle::C,
})
}
fn get_bounds_info(&mut self, b: BoundsInfo<'_>) -> bool {
b.x_l.copy_from_slice(&[-10.0, -10.0]);
b.x_u.copy_from_slice(&[10.0, 10.0]);
b.g_l[0] = 1.0;
b.g_u[0] = 1.0;
true
}
fn get_starting_point(&mut self, sp: StartingPoint<'_>) -> bool {
sp.x.copy_from_slice(&[0.05, 0.05]);
true
}
fn eval_f(&mut self, x: &[Number], _new_x: bool) -> Option<Number> {
Some((x[0] - 3.0).powi(2) + (x[1] - 3.0).powi(2))
}
fn eval_grad_f(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
g[0] = 2.0 * (x[0] - 3.0);
g[1] = 2.0 * (x[1] - 3.0);
true
}
fn eval_g(&mut self, x: &[Number], _new_x: bool, g: &mut [Number]) -> bool {
g[0] = x[0] * x[0] + x[1] * x[1];
true
}
fn eval_jac_g(
&mut self,
x: Option<&[Number]>,
_new_x: bool,
mode: SparsityRequest<'_>,
) -> bool {
match mode {
SparsityRequest::Structure { irow, jcol } => {
irow.copy_from_slice(&[0, 0]);
jcol.copy_from_slice(&[0, 1]);
}
SparsityRequest::Values { values } => {
let x = x.expect("eval_jac_g(Values) without x");
values[0] = 2.0 * x[0];
values[1] = 2.0 * x[1];
}
}
true
}
fn eval_h(
&mut self,
_x: Option<&[Number]>,
_new_x: bool,
obj_factor: Number,
lambda: Option<&[Number]>,
_new_lambda: bool,
mode: SparsityRequest<'_>,
) -> bool {
match mode {
SparsityRequest::Structure { irow, jcol } => {
irow.copy_from_slice(&[0, 1]);
jcol.copy_from_slice(&[0, 1]);
}
SparsityRequest::Values { values } => {
let lam = lambda.expect("eval_h(Values) without lambda");
values[0] = obj_factor * 2.0 + lam[0] * 2.0;
values[1] = obj_factor * 2.0 + lam[0] * 2.0;
}
}
true
}
fn finalize_solution(&mut self, _sol: Solution<'_>, _d: &IpoptData, _q: &IpoptCq) {}
}
fn solve_with_ls_init() -> IpoptApplication {
let mut app = IpoptApplication::new();
// Set the option under test directly. Before gh#604 registered it
// this had to go through `mehrotra_algorithm=yes`, which turns on
// `least_square_init_primal` as part of a cascade that also
// rewrites `bound_push`, `bound_frac` and `bound_mult_init_val` —
// three confounds in a test about one code path. (The report is
// the same either way here: this model's `x` bounds are far enough
// from every point the safeguard visits that neither push moves
// it. The point is that the test no longer depends on that.)
app.options_mut()
.set_string_value("least_square_init_primal", "yes", true, false)
.unwrap();
// The initializer runs once, before iteration 1, so the report is
// fully populated no matter where the solve goes afterwards. Cap
// the iterations and silence the log: this test is about the
// starting point, not about what Mehrotra does with it.
app.options_mut()
.set_integer_value("max_iter", 5, true, false)
.unwrap();
app.options_mut()
.set_integer_value("print_level", 0, true, false)
.unwrap();
app.initialize().unwrap();
let tnlp: Rc<RefCell<dyn TNLP>> = Rc::new(RefCell::new(PoorLinearization));
let _ = app.optimize_tnlp(tnlp);
app
}
/// The unsafeguarded step lands at `(5.025, 5.025)`, whose true
/// violation is `2*5.025^2 - 1 = 49.50125`. Anything at or above this
/// means the raw linearized point was taken.
const UNSAFEGUARDED_VIOLATION: Number = 49.50125;
#[test]
fn least_square_init_never_worsens_the_nonlinear_violation() {
let app = solve_with_ls_init();
let report = app
.least_square_init_report()
.expect("least_square_init_primal ran, so it must report");
eprintln!("gh605: {report:?}");
assert!(
(report.violation_initial - 0.995).abs() < 1e-6,
"expected theta0 = 0.995 at (0.05, 0.05), got {}",
report.violation_initial,
);
// The contract: the point handed to the algorithm is never worse
// than the one the user gave.
assert!(
report.violation_final <= report.violation_initial,
"initializer made the starting point WORSE: {} -> {}",
report.violation_initial,
report.violation_final,
);
// And specifically, nowhere near the unsafeguarded landing point.
assert!(
report.violation_final < UNSAFEGUARDED_VIOLATION * 0.5,
"initializer accepted (or nearly accepted) the raw linearized \
step: violation_final = {} (unsafeguarded lands at {})",
report.violation_final,
UNSAFEGUARDED_VIOLATION,
);
}
#[test]
fn poor_linearization_backtracks_before_accepting() {
let app = solve_with_ls_init();
let report = app
.least_square_init_report()
.expect("least_square_init_primal ran, so it must report");
// alpha = 1 lands at violation ~49.5, alpha = 1/2 at ~11.8,
// alpha = 1/4 at ~2.4 -- all worse than 0.995. The first trial
// that improves is alpha = 1/8.
assert!(
report.rejected_trials > 0,
"expected the full-length step to be rejected, but {} trials \
were rejected (report: {report:?})",
report.rejected_trials,
);
assert!(
report.alpha < 1.0,
"expected a backtracked step, got alpha = {}",
report.alpha,
);
assert_eq!(report.termination, "accepted");
assert!(
report.step_norm > 0.0,
"an accepted step must have a positive norm",
);
}