ledge_core/
problem.rs

1//! Convex QP data structures.
2
3use thiserror::Error;
4
5use crate::matrix::{dot, Matrix};
6
7/// Factor covariance matrix \(\Omega\).
8#[derive(Clone, Debug, PartialEq)]
9#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
10pub enum FactorCovariance {
11    /// Diagonal factor covariance.
12    Diagonal(Vec<f64>),
13    /// Dense symmetric positive-semidefinite factor covariance.
14    Dense(Matrix),
15}
16
17impl FactorCovariance {
18    /// Number of factors represented by this covariance.
19    #[must_use]
20    pub fn dimension(&self) -> usize {
21        match self {
22            Self::Diagonal(diagonal) => diagonal.len(),
23            Self::Dense(matrix) => matrix.rows(),
24        }
25    }
26}
27
28/// A positive-semidefinite quadratic represented as
29/// \(Q = F\Omega F^\mathsf{T} + \operatorname{diag}(d)\).
30#[derive(Clone, Debug, PartialEq)]
31#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
32pub struct FactorQuad {
33    /// Asset-by-factor exposure matrix.
34    pub factors: Matrix,
35    /// Factor covariance.
36    pub omega: FactorCovariance,
37    /// Non-negative idiosyncratic diagonal.
38    pub diagonal: Vec<f64>,
39}
40
41impl FactorQuad {
42    /// Creates a factor quadratic after checking dimensions and convexity data.
43    ///
44    /// Positive semidefiniteness of a dense `omega` is checked during solver
45    /// setup, where its Cholesky-like factor is needed.
46    ///
47    /// # Errors
48    ///
49    /// Returns a [`ProblemError`] for inconsistent dimensions, non-finite
50    /// values, asymmetry, or negative diagonal entries.
51    pub fn new(
52        factors: Matrix,
53        omega: FactorCovariance,
54        diagonal: Vec<f64>,
55    ) -> Result<Self, ProblemError> {
56        let quadratic = Self {
57            factors,
58            omega,
59            diagonal,
60        };
61        quadratic.validate()?;
62        Ok(quadratic)
63    }
64
65    /// Number of decision variables.
66    #[must_use]
67    pub fn dimension(&self) -> usize {
68        self.factors.rows()
69    }
70
71    /// Number of latent factors.
72    #[must_use]
73    pub fn factor_count(&self) -> usize {
74        self.factors.cols()
75    }
76
77    /// Computes `Q * x` without materializing `Q`.
78    #[must_use]
79    pub fn apply(&self, x: &[f64]) -> Vec<f64> {
80        debug_assert_eq!(x.len(), self.dimension());
81        let n = self.dimension();
82        let k = self.factor_count();
83        let mut factor_projection = vec![0.0; k];
84        self.factors.transpose_mul_add(x, &mut factor_projection);
85
86        let weighted = match &self.omega {
87            FactorCovariance::Diagonal(diagonal) => factor_projection
88                .iter()
89                .zip(diagonal)
90                .map(|(value, weight)| value * weight)
91                .collect(),
92            FactorCovariance::Dense(matrix) => matrix.mul_vec(&factor_projection),
93        };
94
95        let mut result: Vec<f64> = self
96            .diagonal
97            .iter()
98            .zip(x)
99            .map(|(diagonal, value)| diagonal * value)
100            .collect();
101        for (row, value) in result.iter_mut().enumerate().take(n) {
102            *value += dot(self.factors.row(row), &weighted);
103        }
104        result
105    }
106
107    fn validate(&self) -> Result<(), ProblemError> {
108        let n = self.factors.rows();
109        let k = self.factors.cols();
110        if self.diagonal.len() != n {
111            return Err(ProblemError::Dimension {
112                field: "quadratic.diagonal",
113                expected: n,
114                actual: self.diagonal.len(),
115            });
116        }
117        if self.omega.dimension() != k {
118            return Err(ProblemError::Dimension {
119                field: "quadratic.omega",
120                expected: k,
121                actual: self.omega.dimension(),
122            });
123        }
124        if self
125            .factors
126            .as_slice()
127            .iter()
128            .any(|value| !value.is_finite())
129        {
130            return Err(ProblemError::NonFinite("quadratic.factors"));
131        }
132        if self
133            .diagonal
134            .iter()
135            .any(|value| !value.is_finite() || *value < 0.0)
136        {
137            return Err(ProblemError::NonConvex(
138                "quadratic diagonal must be finite and non-negative",
139            ));
140        }
141        match &self.omega {
142            FactorCovariance::Diagonal(diagonal) => {
143                if diagonal
144                    .iter()
145                    .any(|value| !value.is_finite() || *value < 0.0)
146                {
147                    return Err(ProblemError::NonConvex(
148                        "factor covariance diagonal must be finite and non-negative",
149                    ));
150                }
151            }
152            FactorCovariance::Dense(matrix) => {
153                if matrix.rows() != matrix.cols() {
154                    return Err(ProblemError::NotSquare("quadratic.omega"));
155                }
156                if matrix.as_slice().iter().any(|value| !value.is_finite()) {
157                    return Err(ProblemError::NonFinite("quadratic.omega"));
158                }
159                for row in 0..k {
160                    for col in 0..row {
161                        let scale = 1.0_f64
162                            .max(matrix[(row, col)].abs())
163                            .max(matrix[(col, row)].abs());
164                        if (matrix[(row, col)] - matrix[(col, row)]).abs() > 1.0e-12 * scale {
165                            return Err(ProblemError::NotSymmetric("quadratic.omega"));
166                        }
167                    }
168                }
169            }
170        }
171        Ok(())
172    }
173}
174
175/// A piecewise-linear proportional cost
176/// `sum_i costs[i] * |x[i] - anchor[i]|` added to the smooth objective
177/// (for portfolios: exact L1 turnover around the previous weights).
178///
179/// Costs must be non-negative. The term is handled by a dedicated
180/// soft-threshold proximal block inside the solver, so it never grows the
181/// SMW-reduced dimension — unlike an epigraph reformulation, which would add
182/// \(2n\) general constraint rows (see `docs/algorithm.md` §4).
183#[derive(Clone, Debug, PartialEq)]
184#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
185pub struct L1Term {
186    /// Non-negative per-variable cost coefficients \(c\).
187    pub costs: Vec<f64>,
188    /// Anchor point \(a\) the absolute deviations are measured from.
189    pub anchor: Vec<f64>,
190}
191
192impl L1Term {
193    /// Evaluates `sum_i costs[i] * |x[i] - anchor[i]|`.
194    #[must_use]
195    pub fn evaluate(&self, x: &[f64]) -> f64 {
196        self.costs
197            .iter()
198            .zip(&self.anchor)
199            .zip(x)
200            .map(|((cost, anchor), value)| cost * (value - anchor).abs())
201            .sum()
202    }
203
204    fn validate(&self, dimension: usize) -> Result<(), ProblemError> {
205        if self.costs.len() != dimension {
206            return Err(ProblemError::Dimension {
207                field: "l1.costs",
208                expected: dimension,
209                actual: self.costs.len(),
210            });
211        }
212        if self.anchor.len() != dimension {
213            return Err(ProblemError::Dimension {
214                field: "l1.anchor",
215                expected: dimension,
216                actual: self.anchor.len(),
217            });
218        }
219        if self
220            .costs
221            .iter()
222            .any(|value| !value.is_finite() || *value < 0.0)
223        {
224            return Err(ProblemError::NonConvex(
225                "l1 costs must be finite and non-negative",
226            ));
227        }
228        if self.anchor.iter().any(|value| !value.is_finite()) {
229            return Err(ProblemError::NonFinite("l1.anchor"));
230        }
231        Ok(())
232    }
233}
234
235/// A block of constraints `matrix * x = rhs` or `matrix * x <= rhs`.
236#[derive(Clone, Debug, PartialEq)]
237#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
238pub struct LinearConstraints {
239    /// Constraint matrix.
240    pub matrix: Matrix,
241    /// Right-hand side.
242    pub rhs: Vec<f64>,
243}
244
245impl LinearConstraints {
246    /// Creates a validated constraint block.
247    ///
248    /// # Errors
249    ///
250    /// Returns [`ProblemError::Dimension`] if row and RHS counts differ.
251    pub fn new(matrix: Matrix, rhs: Vec<f64>) -> Result<Self, ProblemError> {
252        if matrix.rows() != rhs.len() {
253            return Err(ProblemError::Dimension {
254                field: "constraints.rhs",
255                expected: matrix.rows(),
256                actual: rhs.len(),
257            });
258        }
259        Ok(Self { matrix, rhs })
260    }
261
262    /// An empty block with the requested decision dimension.
263    #[must_use]
264    pub fn empty(dimension: usize) -> Self {
265        Self {
266            matrix: Matrix::zeros(0, dimension),
267            rhs: Vec::new(),
268        }
269    }
270
271    /// Number of constraint rows.
272    #[must_use]
273    pub fn len(&self) -> usize {
274        self.rhs.len()
275    }
276
277    /// Whether this block contains no rows.
278    #[must_use]
279    pub fn is_empty(&self) -> bool {
280        self.rhs.is_empty()
281    }
282}
283
284/// Standard-form convex QP accepted by Ledge.
285///
286/// The objective is `0.5 * x' Q x + linear' x` plus an optional
287/// piecewise-linear term `sum_i costs[i] * |x[i] - anchor[i]|`
288/// ([`L1Term`]), with equality, upper-inequality, and box constraints.
289#[derive(Clone, Debug, PartialEq)]
290#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
291pub struct QpProblem {
292    /// Structured positive-semidefinite quadratic.
293    pub quadratic: FactorQuad,
294    /// Linear objective coefficient.
295    pub linear: Vec<f64>,
296    /// Optional proportional-cost term `sum_i costs[i] * |x[i] - anchor[i]|`.
297    pub l1: Option<L1Term>,
298    /// Constraints `A_eq * x = b_eq`.
299    pub equalities: LinearConstraints,
300    /// Constraints `A_ineq * x <= b_ineq`.
301    pub inequalities: LinearConstraints,
302    /// Variable lower bounds; `-inf` means unbounded.
303    ///
304    /// Serialized as `Option<f64>` entries (`null` in JSON) for unbounded
305    /// sides, because JSON cannot represent infinities — see
306    /// `serde_support`.
307    #[cfg_attr(feature = "serde", serde(with = "crate::serde_support::lower_bounds"))]
308    pub lower_bounds: Vec<f64>,
309    /// Variable upper bounds; `+inf` means unbounded.
310    ///
311    /// Serialized as `Option<f64>` entries (`null` in JSON) for unbounded
312    /// sides, because JSON cannot represent infinities — see
313    /// `serde_support`.
314    #[cfg_attr(feature = "serde", serde(with = "crate::serde_support::upper_bounds"))]
315    pub upper_bounds: Vec<f64>,
316}
317
318impl QpProblem {
319    /// Validates dimensions, finite coefficients, and bound consistency.
320    ///
321    /// # Errors
322    ///
323    /// Returns a [`ProblemError`] when the problem is malformed.
324    pub fn validate(&self) -> Result<(), ProblemError> {
325        self.quadratic.validate()?;
326        let n = self.quadratic.dimension();
327        for (field, actual) in [
328            ("linear", self.linear.len()),
329            ("lower_bounds", self.lower_bounds.len()),
330            ("upper_bounds", self.upper_bounds.len()),
331        ] {
332            if actual != n {
333                return Err(ProblemError::Dimension {
334                    field,
335                    expected: n,
336                    actual,
337                });
338            }
339        }
340        for (field, constraints) in [
341            ("equalities", &self.equalities),
342            ("inequalities", &self.inequalities),
343        ] {
344            if constraints.matrix.cols() != n {
345                return Err(ProblemError::Dimension {
346                    field,
347                    expected: n,
348                    actual: constraints.matrix.cols(),
349                });
350            }
351            if constraints.matrix.rows() != constraints.rhs.len() {
352                return Err(ProblemError::Dimension {
353                    field,
354                    expected: constraints.matrix.rows(),
355                    actual: constraints.rhs.len(),
356                });
357            }
358            if constraints
359                .matrix
360                .as_slice()
361                .iter()
362                .chain(&constraints.rhs)
363                .any(|value| !value.is_finite())
364            {
365                return Err(ProblemError::NonFinite(field));
366            }
367        }
368        if self.linear.iter().any(|value| !value.is_finite()) {
369            return Err(ProblemError::NonFinite("linear"));
370        }
371        if let Some(l1) = &self.l1 {
372            l1.validate(n)?;
373        }
374        for index in 0..n {
375            let lower = self.lower_bounds[index];
376            let upper = self.upper_bounds[index];
377            if lower.is_nan() || upper.is_nan() {
378                return Err(ProblemError::NonFinite("bounds"));
379            }
380            if lower > upper {
381                return Err(ProblemError::InvalidBounds {
382                    index,
383                    lower,
384                    upper,
385                });
386            }
387        }
388        Ok(())
389    }
390
391    /// Evaluates the objective at `x`, including the L1 term when present.
392    #[must_use]
393    pub fn objective(&self, x: &[f64]) -> f64 {
394        let qx = self.quadratic.apply(x);
395        let smooth = 0.5 * dot(x, &qx) + dot(&self.linear, x);
396        smooth + self.l1.as_ref().map_or(0.0, |term| term.evaluate(x))
397    }
398}
399
400/// Input validation errors.
401#[derive(Debug, Error, Clone, PartialEq)]
402pub enum ProblemError {
403    /// A vector or matrix dimension is inconsistent.
404    #[error("{field} has dimension {actual}; expected {expected}")]
405    Dimension {
406        /// Name of the invalid field.
407        field: &'static str,
408        /// Expected dimension.
409        expected: usize,
410        /// Actual dimension.
411        actual: usize,
412    },
413    /// A coefficient that must be finite is NaN or infinite.
414    #[error("{0} contains a non-finite coefficient")]
415    NonFinite(&'static str),
416    /// Convexity data is invalid.
417    #[error("non-convex quadratic data: {0}")]
418    NonConvex(&'static str),
419    /// A matrix expected to be square is not.
420    #[error("{0} must be square")]
421    NotSquare(&'static str),
422    /// A matrix expected to be symmetric is not.
423    #[error("{0} must be symmetric")]
424    NotSymmetric(&'static str),
425    /// Lower and upper bounds are inconsistent.
426    #[error("invalid bounds at index {index}: lower {lower} exceeds upper {upper}")]
427    InvalidBounds {
428        /// Variable index.
429        index: usize,
430        /// Supplied lower bound.
431        lower: f64,
432        /// Supplied upper bound.
433        upper: f64,
434    },
435}
ledge_core/problem.rs

ledge_core/
problem.rs
