herculesabqp 0.1.2

A convex box-constrained quadratic programming solver with warm starts and active-set polishing.
Documentation 
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

use ndarray::Array1;

use crate::matrix::QuadraticMatrix;

/// Matrix-vector interface for quadratic operators used by the first-order solver.
///
/// This trait is enough to support accelerated projected-gradient solves without
/// requiring explicit access to matrix entries or reduced submatrix extraction.
/// Implementations may optionally provide diagonal data or a Gershgorin bound
/// to enable Hessian-diagonal scaling and cheaper Lipschitz estimation.
///
/// The intended contract is:
/// - `matvec_into(x, out)` writes `Qx` into `out`
/// - `n()` returns the dimension of the square operator
/// - `diagonal()` is optional, but enables Hessian-diagonal scaling for the
///   implicit solve path
/// - `gershgorin_upper_bound()` is optional, and when absent the solver falls
///   back to the power-iteration estimate for the Lipschitz constant
pub trait QuadraticOperator: Send + Sync {
    /// Return the dimension of the square quadratic operator.
    fn n(&self) -> usize;

    /// Compute `out = Q x`.
    fn matvec_into(&self, x: &Array1<f64>, out: &mut Array1<f64>);

    /// Return the diagonal of `Q` when cheaply available.
    fn diagonal(&self) -> Option<Array1<f64>> {
        None
    }

    /// Return a Gershgorin upper bound for the largest eigenvalue when available.
    fn gershgorin_upper_bound(&self) -> Option<f64> {
        None
    }
}

impl QuadraticOperator for QuadraticMatrix {
    fn n(&self) -> usize {
        self.n()
    }

    fn matvec_into(&self, x: &Array1<f64>, out: &mut Array1<f64>) {
        self.matvec_into(x, out);
    }

    fn diagonal(&self) -> Option<Array1<f64>> {
        Some(self.diagonal())
    }

    fn gershgorin_upper_bound(&self) -> Option<f64> {
        Some(match self {
            QuadraticMatrix::Dense(m) => {
                let mut best: f64 = 1e-12;
                for i in 0..m.nrows() {
                    let diag = m[[i, i]];
                    let abs_row_sum = (0..m.ncols()).map(|j| m[[i, j]].abs()).sum::<f64>();
                    best = best.max(diag + abs_row_sum - diag.abs());
                }
                best.max(1e-12)
            }
            QuadraticMatrix::Sparse(m) => {
                let diag = m.diag();
                let mut best: f64 = 1e-12;
                for (i, row) in m.outer_iterator().enumerate() {
                    let abs_row_sum = row.iter().map(|(_, v)| v.abs()).sum::<f64>();
                    best = best.max(diag[i] + abs_row_sum - diag[i].abs());
                }
                best.max(1e-12)
            }
        })
    }
}