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
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

import math

import numpy as np
import pytest
import scipy.sparse as sp

import herculesabqp


class DiagonalOperator:
    def __init__(self, diag):
        self.diag = np.asarray(diag, dtype=np.float64)
        self.n = int(self.diag.size)

    def matvec(self, x):
        x = np.asarray(x, dtype=np.float64)
        return self.diag * x

    def diagonal(self):
        return self.diag

    def gershgorin_upper_bound(self):
        return float(np.max(self.diag))


class IdentityOperator:
    def __init__(self, n):
        self.n = int(n)

    def matvec(self, x):
        return np.asarray(x, dtype=np.float64)


@pytest.fixture
def explicit_problem():
    q = np.array([[4.0, 1.0], [1.0, 3.0]], dtype=np.float64)
    c = np.array([-1.0, -0.5], dtype=np.float64)
    lb = np.zeros(2, dtype=np.float64)
    ub = np.ones(2, dtype=np.float64)
    return q, c, lb, ub


@pytest.fixture
def sparse_explicit_problem(explicit_problem):
    q, c, lb, ub = explicit_problem
    return sp.csr_matrix(q), c, lb, ub


@pytest.fixture
def implicit_problem():
    op = DiagonalOperator([2.0, 2.0, 2.0])
    c = np.array([-1.0, 0.5, -3.0], dtype=np.float64)
    lb = np.zeros(3, dtype=np.float64)
    ub = np.ones(3, dtype=np.float64)
    return op, c, lb, ub


@pytest.fixture
def prepared_implicit_problem():
    op = DiagonalOperator([4.0, 9.0])
    c = np.array([-2.0, 3.0], dtype=np.float64)
    lb = np.zeros(2, dtype=np.float64)
    ub = np.ones(2, dtype=np.float64)
    return op, c, lb, ub


def test_explicit_dense_solve(explicit_problem):
    q, c, lb, ub = explicit_problem

    result = herculesabqp.solve_box_qp(
        q,
        c,
        lb,
        ub,
        assume_symmetric=True,
        scaling="none",
    )

    assert result["kkt_inf"] < 1e-5
    assert result["gap"] >= 0.0
    assert math.isfinite(result["certified_lower_bound"])
    assert abs(
        result["objective"] - (result["certified_lower_bound"] + result["gap"])
    ) < 1e-8
    assert 0.0 <= result["x"][0] <= 1.0
    assert 0.0 <= result["x"][1] <= 1.0


def test_explicit_prepared_sparse_solve(sparse_explicit_problem):
    q, c, lb, ub = sparse_explicit_problem

    solver = herculesabqp.PreparedSolver(
        q,
        c,
        assume_symmetric=True,
        scaling="none",
    )
    result = solver.solve(lb, ub)

    assert result["kkt_inf"] < 1e-5
    assert len(result["x"]) == 2


def test_implicit_solve(implicit_problem):
    op, c, lb, ub = implicit_problem

    result = herculesabqp.solve_box_qp_implicit(
        op,
        c,
        lb,
        ub,
        assume_symmetric=True,
        scaling="hessian_diag",
    )

    assert result["kkt_inf"] < 1e-5
    assert result["scaling_applied"]
    assert result["scaling_name"] == "hessian_diag"


def test_implicit_falls_back_to_unscaled_without_diagonal():
    op = IdentityOperator(2)
    c = np.array([-0.25, 0.5], dtype=np.float64)
    lb = np.zeros(2, dtype=np.float64)
    ub = np.ones(2, dtype=np.float64)

    result = herculesabqp.solve_box_qp_implicit(
        op,
        c,
        lb,
        ub,
        assume_symmetric=True,
        scaling="hessian_diag",
    )

    assert not result["scaling_applied"]
    assert result["scaling_name"] == "none"


def test_prepared_implicit_solve(prepared_implicit_problem):
    op, c, lb, ub = prepared_implicit_problem

    solver = herculesabqp.PreparedImplicitSolver(
        op,
        c,
        assume_symmetric=True,
        scaling="hessian_diag",
    )
    result = solver.solve(lb, ub)

    assert math.isfinite(result["objective"])
    assert result["scaling_applied"]