Pardalos & Kovoor's O(n) solver for singly-constrained bounded QPs

pardalos_kovoor.py

# Author: Vlad Niculae <vlad@vene.ro>
# License: Simplified BSD

import numpy as np

try:
    from numba import jit
except ImportError:
    print("numba not available")
    def jit(nopython):
        def id(f):
            return f
        return id


@jit(nopython=True)
def solve_pardalos(c, d, a, b):
    """Solve the following constrained QP

    minimize  sum_i c_i x_i^2
          st  sum_i c_i x_i = d
              a[i] <= x_i <= b[i]

    Parameters
    ----------
    c, a, b: ndarrays of length n
    d: scalar

    Returns
    -------
    x: ndarray of length n

    Notes
    -----
    Naive / straightforward implementation based on:

        An algorithm for a singly constrained class of quadratic
        programs subject to upper and lower bounds.
        P. M. Pardalos and N. Kovoor. 1990.
        https://link.springer.com/article/10.1007/BF01585748

    Assuming np.median() is efficient, this should be O(n)

    """

    UnsetV = list(range(len(c)))

    NINF = -np.inf
    PINF = np.inf

    Intervalpts = [NINF, PINF]

    for i in range(len(c)):
        Intervalpts.append(a[i])
        Intervalpts.append(b[i])

    Min = NINF
    Max = PINF
    Tightsum = 0
    Slackweight = 0

    while len(UnsetV) > 0:

        Intervalpts_np = np.array(Intervalpts)
        Mid = np.median(Intervalpts_np)
        Testsum = Tightsum
        for i in UnsetV:
            if b[i] < Mid:
                Testsum += c[i] * b[i]
            if a[i] > Mid:
                Testsum += c[i] * a[i]

        tmp = 0
        for i in UnsetV:
            if a[i] <= Mid <= b[i]:
                tmp += c[i]
        Testsum += (Slackweight + tmp) * Mid

        if Testsum <= d:
            Min = Mid
        if Testsum >= d:
            Max = Mid

        # below differs from original paper: strict inequalities.
        # I encountered infinite cycling otherwise.
        Intervalpts = [ppp for ppp in Intervalpts if  Min < ppp < Max]

        for i in UnsetV:
            if b[i] <= Min:
                Tightsum += c[i] * b[i]
            if a[i] >= Max:
                Tightsum += c[i] * a[i]

            if a[i] <= Min and b[i] >= Max:
                Slackweight += c[i]
        UnsetV = [i for i in UnsetV if (Min < a[i] < Max) or (Min < b[i] < Max)]

    if Slackweight != 0 and (d - Tightsum) != 0:
        tau = (d - Tightsum) / Slackweight
    else:
        tau = Mid
    return np.maximum(a, np.minimum(b, tau))


def project_constraints(theta, z, lower=0, upper=1):
    """Use the Pardalos & Kovoor algorithm to compute the euclidean projection:

    argmin || y - theta ||
    st     sum_i y_i = z
           lower <= y_i <= upper

    Uses a simple change of variable.

    Parameters
    ----------
    theta: ndarray, length n
    lower, upper: scalars or ndarrays of length n
    z: scalar

    Returns
    -------
    y: ndarray, length n
    """
    a = (lower - theta) / 2
    b = (upper - theta) / 2
    c = np.ones_like(theta)
    d = (z - np.sum(theta)) / 2

    out = solve_pardalos(c, d, a, b)
    return 2 * out + theta


if __name__ == '__main__':

    from time import clock
    try:
        import matplotlib.pyplot as plt
        plot = True
    except ImportError:
        plot = False

    rng = np.random.RandomState(0)

    # run once to initialize numba
    theta = rng.randn(10)
    project_constraints(theta, 4)

    n_runs = 100
    n_problems = 20

    ns = np.arange(1000, 10000, step=500)
    all_times = [[] for _ in ns]

    for n, times in zip(ns, all_times):

        z = n // 10

        for _ in range(n_problems):  # solve multiple problems
            theta = rng.randn(n)
            tic = clock()
            for _ in range(n_runs):
                project_constraints(theta, z)
            toc = clock()
            times.append((toc - tic) / n_runs)

    if plot:

        medians = []
        intervals = []

        for times in all_times:
            p5, median, p95 = np.percentile(times, [5, 50, 95])
            medians.append(median)
            intervals.append((p5, p95))

        medians = np.array(medians)
        intervals = medians - np.array(intervals).T
        intervals[1] *= -1

        plt.errorbar(ns, medians, yerr=intervals)

        plt.xlabel("n (problem size)")
        plt.ylabel("time (s)")
        plt.show()

    else:

        for n, times in zip(ns, all_times):
            p5, median, p95 = np.percentile(times, [5, 50, 95])
            print(n, median, p5, p95)

test_pardalos_kovoor.py

import numpy as np
import cvxpy as cvx

from pardalos_kovoor import project_constraints

theta_problematic = np.array([-0.24872674, 0.3, 0.3, 0.3, 0.3, -0.24872674, 0.3, 0.3, 0.3, 0.3])


def project_cvxpy(theta, z, lower, upper):
    n = len(theta)
    y = cvx.Variable(n)
    objective = cvx.Minimize(cvx.sum_squares(y - theta))
    constraints = [lower <= y, y <= upper, cvx.sum(y) == z]
    prob = cvx.Problem(objective, constraints)
    res = prob.solve()
    return y.value

    
# check projection to self
def test_feasible():
    x = np.zeros(10)
    x[0] = 1
    x[4] = 1

    print(x)
    print(project_constraints(x, z=2, lower=0, upper=1))


def test_scaled():
    x = np.zeros(10)
    x[0] = 2 
    x[4] = 2 

    print(x)
    print(project_constraints(x, z=2, lower=0, upper=1))


def test_terminate():
    project_constraints(theta_problematic, z=3, lower=0, upper=1)


def _compare(x, z, lower, upper):
    y_cvx = project_cvxpy(x, z=3, lower=0, upper=1)
    y_our = project_constraints(x, z=3, lower=0, upper=1)
    err = np.linalg.norm(y_our - y_cvx)
    if err > 1e-5:
        print("Error:", err)
        # print("x=", x)
        # print("cvxpy:", y_cvx.round(2), np.sum(y_cvx), np.linalg.norm(x - y_cvx))
        # print("our impl", y_our.round(2), np.sum(y_our), np.linalg.norm(x - y_our))
        print("cvxpy:", np.sum(y_cvx), np.linalg.norm(x - y_cvx))
        print("our impl", np.sum(y_our), np.linalg.norm(x - y_our))


def test_numeric():
    rng = np.random.RandomState(0)
    _compare(theta_problematic, z=3, lower=0, upper=1)

    for _ in range(10000):
       x = rng.randn(6)
       _compare(x, z=3, lower=0, upper=1)


test_feasible()
test_scaled()
test_terminate()
test_numeric()

(excuse the variable names: I followed the ones from the reference paper verbatim)