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

use super::operator::QuadraticOperator;
use super::types::LipschitzMethod;
use super::utils::{DenseVec, l2_norm};

const AUTO_POWER_ITERS: usize = 30;
const AUTO_POWER_SAFETY: f64 = 1.2;

/// Estimate a safe gradient Lipschitz constant for the quadratic objective.
pub fn estimate_lipschitz(q: &dyn QuadraticOperator, method: &LipschitzMethod) -> f64 {
    match method {
        LipschitzMethod::Gershgorin => estimate_lipschitz_gershgorin(q),
        LipschitzMethod::Auto => {
            let power = AUTO_POWER_SAFETY * estimate_lipschitz_power(q, AUTO_POWER_ITERS);
            match q.gershgorin_upper_bound() {
                Some(gershgorin) => gershgorin.min(power),
                None => power,
            }
        }
    }
}

/// Gershgorin upper bound for the largest eigenvalue of `Q`.
pub fn estimate_lipschitz_gershgorin(q: &dyn QuadraticOperator) -> f64 {
    q.gershgorin_upper_bound()
        .unwrap_or_else(|| estimate_lipschitz_power(q, AUTO_POWER_ITERS))
}

// Refine the largest-eigenvalue estimate with a short power iteration.
fn estimate_lipschitz_power(q: &dyn QuadraticOperator, n_iter: usize) -> f64 {
    let n = q.n();
    let mut v = DenseVec::from_elem(n, 1.0 / (n as f64).sqrt());
    let mut w = DenseVec::zeros(n);
    for _ in 0..n_iter {
        q.matvec_into(&v, &mut w);
        let nw = l2_norm(&w);
        if nw == 0.0 {
            return 1.0;
        }
        for i in 0..n {
            v[i] = w[i] / nw;
        }
    }
    q.matvec_into(&v, &mut w);
    v.dot(&w).max(1e-12)
}