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
//! Karush-Kuhn-Tucker residual checks.
use thiserror::Error;
use crate::{matrix::norm_inf, problem::QpProblem};
/// Dual multipliers associated with each constraint block.
///
/// Inequality multipliers use the convention `A_ineq * x <= b_ineq` and are
/// therefore non-negative. A box multiplier is negative at a lower bound and
/// positive at an upper bound.
#[derive(Clone, Debug, Default, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct DualVariables {
/// Equality multipliers.
pub equalities: Vec<f64>,
/// Upper-inequality multipliers.
pub inequalities: Vec<f64>,
/// Combined normal-cone multipliers for variable boxes.
pub bounds: Vec<f64>,
/// Subgradient multipliers of the [`L1Term`](crate::L1Term); empty when
/// the problem has none.
///
/// Each entry lies in `[-costs[i], costs[i]]` and equals the signed cost
/// wherever the variable moved off the anchor (`+costs[i]` above it,
/// `-costs[i]` below).
pub l1: Vec<f64>,
}
/// Infinity-norm KKT diagnostics.
#[derive(Clone, Copy, Debug, Default, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct KktResiduals {
/// Maximum equality, inequality, or bound violation.
pub primal: f64,
/// Maximum stationarity or dual-cone violation.
pub dual: f64,
/// Maximum absolute complementary-slackness product.
pub complementarity: f64,
}
/// Errors from the standalone KKT checker.
#[derive(Clone, Debug, Error, PartialEq, Eq)]
pub enum KktError {
/// A primal or dual vector has the wrong length.
#[error("{field} has length {actual}; expected {expected}")]
Dimension {
/// Name of the vector.
field: &'static str,
/// Expected length.
expected: usize,
/// Supplied length.
actual: usize,
},
}
/// Evaluates primal feasibility, stationarity/dual feasibility, and
/// complementarity.
///
/// # Errors
///
/// Returns [`KktError::Dimension`] if a primal or multiplier vector has the
/// wrong length.
pub fn check_kkt(
problem: &QpProblem,
x: &[f64],
dual: &DualVariables,
) -> Result<KktResiduals, KktError> {
let n = problem.quadratic.dimension();
let l1_expected = if problem.l1.is_some() { n } else { 0 };
for (field, actual, expected) in [
("x", x.len(), n),
(
"dual.equalities",
dual.equalities.len(),
problem.equalities.len(),
),
(
"dual.inequalities",
dual.inequalities.len(),
problem.inequalities.len(),
),
("dual.bounds", dual.bounds.len(), n),
("dual.l1", dual.l1.len(), l1_expected),
] {
if actual != expected {
return Err(KktError::Dimension {
field,
expected,
actual,
});
}
}
let equality_values = problem.equalities.matrix.mul_vec(x);
let inequality_values = problem.inequalities.matrix.mul_vec(x);
let mut primal = 0.0_f64;
for (value, rhs) in equality_values.iter().zip(&problem.equalities.rhs) {
primal = primal.max((value - rhs).abs());
}
for (value, rhs) in inequality_values.iter().zip(&problem.inequalities.rhs) {
primal = primal.max((value - rhs).max(0.0));
}
for (index, value) in x.iter().enumerate() {
primal = primal
.max((problem.lower_bounds[index] - value).max(0.0))
.max((value - problem.upper_bounds[index]).max(0.0));
}
let mut stationarity = problem.quadratic.apply(x);
for (value, linear) in stationarity.iter_mut().zip(&problem.linear) {
*value += linear;
}
problem
.equalities
.matrix
.transpose_mul_add(&dual.equalities, &mut stationarity);
problem
.inequalities
.matrix
.transpose_mul_add(&dual.inequalities, &mut stationarity);
for (value, bound_dual) in stationarity.iter_mut().zip(&dual.bounds) {
*value += bound_dual;
}
for (value, l1_dual) in stationarity.iter_mut().zip(&dual.l1) {
*value += l1_dual;
}
let mut cone_violation = 0.0_f64;
let mut complementarity = 0.0_f64;
for ((value, rhs), multiplier) in inequality_values
.iter()
.zip(&problem.inequalities.rhs)
.zip(&dual.inequalities)
{
cone_violation = cone_violation.max((-multiplier).max(0.0));
complementarity = complementarity.max((multiplier * (value - rhs)).abs());
}
// Box multipliers are scored without an activity window: the positive
// part of a multiplier must pair with a finite upper bound and the
// negative part with a finite lower bound (else it violates the dual
// cone), and each part is charged the epsilon-complementarity product
// `|multiplier| * distance-to-bound`. This keeps the checker fair to
// near-active solutions from interior-point solvers, whose multipliers
// decay smoothly with the distance to the bound instead of switching
// off at machine precision.
for (index, value) in x.iter().copied().enumerate() {
let lower = problem.lower_bounds[index];
let upper = problem.upper_bounds[index];
let multiplier = dual.bounds[index];
let toward_upper = multiplier.max(0.0);
let toward_lower = (-multiplier).max(0.0);
if upper.is_finite() {
complementarity = complementarity.max(toward_upper * (upper - value).abs());
} else {
cone_violation = cone_violation.max(toward_upper);
}
if lower.is_finite() {
complementarity = complementarity.max(toward_lower * (value - lower).abs());
} else {
cone_violation = cone_violation.max(toward_lower);
}
}
// L1 multipliers must lie in the subdifferential of
// `sum_i costs[i] * |x[i] - anchor[i]|`: inside `[-c_i, c_i]` everywhere
// (dual-cone violation otherwise), pinned at the signed cost wherever the
// variable moved off the anchor. The pinning is scored continuously,
// mirroring the box scoring above: the shortfall from the required signed
// cost is charged times the distance moved in that direction.
if let Some(l1) = &problem.l1 {
for (index, multiplier) in dual.l1.iter().copied().enumerate() {
let cost = l1.costs[index];
let offset = x[index] - l1.anchor[index];
cone_violation = cone_violation.max(multiplier.abs() - cost);
complementarity = complementarity
.max(((cost - multiplier) * offset.max(0.0)).abs())
.max(((cost + multiplier) * (-offset).max(0.0)).abs());
}
}
Ok(KktResiduals {
primal,
dual: norm_inf(&stationarity).max(cone_violation),
complementarity,
})
}