Speed issue for Quadratic Program

[vaibhav.u]

I am dealing with a projection problem in a space defined by multiple linear constraints, with a problem dimension of about 1000, and I have 1000 such problems in a ‘for’ loop. I initially used the OSQP solver with JuMP, but then it was too slow. Then I removed JuMP and used SparseArrays to construct and solve the problem, but it’s now taking 60-120 seconds, which is faster than before but still slow. Is there a way to speed up the process? Or is there a better solver I could use for this projection problem?

The constraints are A_eq * x = b_eq, A_upper * x >= p.clipping_threshold .* ones(n_upper) and A_lower * x <= -p.clipping_threshold .* ones(n_lower) and I am trying to minimize norm(x - x_ellipsoid), where x is the decision variable. Originally, I had an ellipsoidal constraint as well. Still, it seems that whenever I project it onto the ellipsoid by scaling, then to the space of linear constraints as given above, it still satisfies the ellipsoidal constraint, so I am not including that in the optimization problem.

Please see the following code for the problem construction:

function build_osqp_solver(
    A_eq::Matrix{Float64},
    A_upper::Matrix{Float64},
    A_lower::Matrix{Float64},
    n::Int,
    p::ModelParams
)
    n_eq    = size(A_eq,    1)
    n_upper = size(A_upper, 1)
    n_lower = size(A_lower, 1)
    n_con   = n_eq + n_upper + n_lower

    P = sparse(2.0 * I, n, n)

    if n_con == 0
        A_con = spzeros(0, n)
    else
        blocks = Matrix{Float64}[]
        n_eq    > 0 && push!(blocks, A_eq)
        n_upper > 0 && push!(blocks, A_upper)
        n_lower > 0 && push!(blocks, A_lower)
        A_con = sparse(vcat(blocks...))
    end

    # Cache l and u — we'll mutate the equality rows in each sample
    l = zeros(n_con)
    u = zeros(n_con)

    if n_upper > 0
        l[n_eq+1 : n_eq+n_upper] .=  p.clipping_threshold
        u[n_eq+1 : n_eq+n_upper] .=  Inf
    end
    if n_lower > 0
        l[n_eq+n_upper+1 : end] .= -Inf
        u[n_eq+n_upper+1 : end] .= -p.clipping_threshold
    end

    solver = OSQP.Model()
    OSQP.setup!(solver;
        P             = P,
        q             = zeros(n),
        A             = A_con,
        l             = l,
        u             = u,
        warm_starting = true,
        verbose       = false,
        eps_abs       = 1e-6,
        eps_rel       = 1e-6,
        max_iter      = 10000,
        polish        = true
    )

    # Return cached l and u so we can mutate them per sample
    return solver, n_eq, l, u
end


function project_coeff_fast(
    x0::Vector{Float64},
    b_eq::Vector{Float64},
    solver::OSQP.Model,
    n_eq::Int,
    l::Vector{Float64},   # mutated in place
    u::Vector{Float64},   # mutated in place
    p::ModelParams
)
    x_ellipsoid = project_onto_ellipsoid(x0, p)
    q = -2.0 .* x_ellipsoid

    if n_eq > 0
        # Update equality rows in the cached vectors, then pass full l and u
        l[1:n_eq] .= b_eq
        u[1:n_eq] .= b_eq
        OSQP.update!(solver; q=q, l=l, u=u)
    else
        OSQP.update!(solver; q=q)
    end

    results = OSQP.solve!(solver)

    if results.info.status != :Solved && results.info.status != :Almost_Solved
        @warn "OSQP status: $(results.info.status). Returning ellipsoid projection."
        return x_ellipsoid
    end

    return results.x
end

[odow]

Hi @vaibhav.u, welcome to the forum

It’s hard to say what is going on without a reproducible example.

What makes you think 60 seconds I slow? Do you have a benchmark in a different language/solver/ecosystem?

Have you benchmarked the code to see where the time is being spent? Is it inside OSQP.solve!? Or is it somewhere else?

[vaibhav.u]

@odow

What makes you think 60 seconds I slow? Do you have a benchmark in a different language/solver/ecosystem?

Hi @odow

Thank you very much for your response. Yes, indeed, it is OSQP.Solve!. I ran the following benchmark code:

# ============================================================
# Reproducible example: QP projection with equality +
# inequality constraints using OSQP.jl
# ============================================================

using OSQP, SparseArrays, LinearAlgebra, Random, BenchmarkTools, Profile

struct ModelParams
    n_data::Int
    window_time::Float64
    clipping_threshold::Float64
    G_normal::Matrix{Float64}
end

function classify_data_points(data::Vector{Float64}, p::ModelParams)
    upper_idx   = findall(x ->  x >= p.clipping_threshold,  data)
    lower_idx   = findall(x ->  x <= -p.clipping_threshold, data)
    trusted_idx = findall(x -> abs(x) < p.clipping_threshold && x != 0.0, data)
    return trusted_idx, upper_idx, lower_idx
end

function project_onto_ellipsoid(x::Vector{Float64}, p::ModelParams)
    energy = x' * p.G_normal * x
    energy <= p.n_data && return x
    return x * sqrt(p.window_time / energy)
end

function build_osqp_solver(
    A_eq::Matrix{Float64},
    A_upper::Matrix{Float64},
    A_lower::Matrix{Float64},
    n::Int,
    p::ModelParams
)
    n_eq    = size(A_eq,    1)
    n_upper = size(A_upper, 1)
    n_lower = size(A_lower, 1)
    n_con   = n_eq + n_upper + n_lower

    P = sparse(2.0 * I, n, n)

    if n_con == 0
        A_con = spzeros(0, n)
    else
        blocks = Matrix{Float64}[]
        n_eq    > 0 && push!(blocks, A_eq)
        n_upper > 0 && push!(blocks, A_upper)
        n_lower > 0 && push!(blocks, A_lower)
        A_con = sparse(vcat(blocks...))
    end

    l = zeros(n_con)
    u = zeros(n_con)

    if n_upper > 0
        l[n_eq+1 : n_eq+n_upper] .=  p.clipping_threshold
        u[n_eq+1 : n_eq+n_upper] .=  Inf
    end
    if n_lower > 0
        l[n_eq+n_upper+1 : end] .= -Inf
        u[n_eq+n_upper+1 : end] .= -p.clipping_threshold
    end

    solver = OSQP.Model()
    OSQP.setup!(solver;
        P             = P,
        q             = zeros(n),
        A             = A_con,
        l             = l,
        u             = u,
        warm_starting = true,
        verbose       = false,
        eps_abs       = 1e-6,
        eps_rel       = 1e-6,
        max_iter      = 10000,
        polish        = true
    )

    return solver, n_eq, l, u
end

function project_single(
    x0::Vector{Float64},
    b_eq::Vector{Float64},
    solver::OSQP.Model,
    n_eq::Int,
    l::Vector{Float64},
    u::Vector{Float64},
    p::ModelParams
)
    x_ellipsoid = project_onto_ellipsoid(x0, p)
    q = -2.0 .* x_ellipsoid

    if n_eq > 0
        l[1:n_eq] .= b_eq
        u[1:n_eq] .= b_eq
        OSQP.update!(solver; q=q, l=l, u=u)
    else
        OSQP.update!(solver; q=q)
    end

    res = OSQP.solve!(solver)

    if res.info.status != :Solved && res.info.status != :Almost_Solved
        @warn "OSQP status: $(res.info.status). Returning ellipsoid projection."
        return x_ellipsoid
    end

    return res.x
end

function data_subspace_projection(
    coefficients::Matrix{Float64},
    data::Vector{Float64},
    basis::Matrix{Float64},
    p::ModelParams
)
    trusted_idx, upper_idx, lower_idx = classify_data_points(data, p)

    n             = 2 * p.n_data
    maxim_samples = size(coefficients, 1)
    result        = similar(coefficients)

    if isempty(trusted_idx) && isempty(upper_idx) && isempty(lower_idx)
        for j in 1:maxim_samples
            result[j, :] = project_onto_ellipsoid(coefficients[j, :], p)
        end
        return result
    end

    A_eq    = isempty(trusted_idx) ? Matrix{Float64}(undef, 0, n) : basis[trusted_idx, :]
    A_upper = isempty(upper_idx)   ? Matrix{Float64}(undef, 0, n) : basis[upper_idx,   :]
    A_lower = isempty(lower_idx)   ? Matrix{Float64}(undef, 0, n) : basis[lower_idx,   :]
    b_eq    = isempty(trusted_idx) ? Float64[]                    : data[trusted_idx]

    solver, n_eq, l, u = build_osqp_solver(A_eq, A_upper, A_lower, n, p)

    ellipsoid_violations = 0
    for j in 1:maxim_samples
        result[j, :] = project_single(coefficients[j, :], b_eq, solver, n_eq, l, u, p)
        if result[j, :]' * p.G_normal * result[j, :] > p.n_data
            ellipsoid_violations += 1
        end
    end

    if ellipsoid_violations > 0
        @warn "Ellipsoid constraint violated for $ellipsoid_violations / $maxim_samples samples."
    end

    return result
end

# ===============================================================
# MAIN
# ===============================================================

Random.seed!(42)

n_data             = 320
n                  = 2 * n_data
maxim_samples      = 2 * n_data + 1
clipping_threshold = 0.5

A_rand   = randn(n, n)
G_normal = (A_rand' * A_rand) / n

p = ModelParams(n_data, 64.0 / 1000.0, clipping_threshold, G_normal)

basis        = randn(n_data, n)
raw_signal   = 0.8 .* randn(n_data)
data         = clamp.(raw_signal, -clipping_threshold, clipping_threshold)
coefficients = 0.1 .* randn(maxim_samples, n)

println("Problem size : $maxim_samples samples × $n coefficients")
println("Trusted pts  : ", count(x -> abs(x) < clipping_threshold && x != 0.0, data))
println("Upper clipped: ", count(x -> x >= clipping_threshold,  data))
println("Lower clipped: ", count(x -> x <= -clipping_threshold, data))

# ---------------------------------------------------------------
# Explicitly build constraint matrices in main scope for benchmarking
# ---------------------------------------------------------------
trusted_idx, upper_idx, lower_idx = classify_data_points(data, p)

A_eq    = isempty(trusted_idx) ? Matrix{Float64}(undef, 0, n) : basis[trusted_idx, :]
A_upper = isempty(upper_idx)   ? Matrix{Float64}(undef, 0, n) : basis[upper_idx,   :]
A_lower = isempty(lower_idx)   ? Matrix{Float64}(undef, 0, n) : basis[lower_idx,   :]
b_eq    = isempty(trusted_idx) ? Float64[]                    : data[trusted_idx]

solver, n_eq, l, u = build_osqp_solver(A_eq, A_upper, A_lower, n, p)
x0 = coefficients[1, :]

# ---------------------------------------------------------------
# Warmup — avoid measuring compilation
# ---------------------------------------------------------------
data_subspace_projection(coefficients, data, basis, p)
project_single(x0, b_eq, solver, n_eq, l, u, p)

# ---------------------------------------------------------------
# 1. Coarse benchmarks
# ---------------------------------------------------------------
println("\n--- Coarse benchmarks ---")
print("build_osqp_solver      : "); @btime build_osqp_solver($A_eq, $A_upper, $A_lower, $n, $p)
print("project_single         : "); @btime project_single($x0, $b_eq, $solver, $n_eq, $l, $u, $p)
print("full loop              : "); @btime data_subspace_projection($coefficients, $data, $basis, $p)

# ---------------------------------------------------------------
# 2. project_single internals
# ---------------------------------------------------------------
x_ellipsoid = project_onto_ellipsoid(x0, p)
q_vec       = -2.0 .* x_ellipsoid
l[1:n_eq]  .= b_eq
u[1:n_eq]  .= b_eq

println("\n--- project_single breakdown ---")
print("  project_onto_ellipsoid : "); @btime project_onto_ellipsoid($x0, $p)
print("  OSQP.update!           : "); @btime OSQP.update!($solver; q=$q_vec, l=$l, u=$u)
print("  OSQP.solve!            : "); @btime OSQP.solve!($solver)

# ---------------------------------------------------------------
# 3. Profiler
# ---------------------------------------------------------------
Profile.clear()
@profile for _ in 1:10
    data_subspace_projection(coefficients, data, basis, p)
end
Profile.print(mincount=5, sortedby=:count)

And for a single projection, the breakdown was:

project_onto_ellipsoid :   62.800 μs (3 allocations: 5.07 KiB)
OSQP.update!           :   2.644 μs (0 allocations: 0 bytes)
OSQP.solve!            :   201.039 ms (14 allocations: 15.50 KiB)

I am solving 641 of these problems, so I was wondering if there is a good QP solver out there that can speed up the computation, or if there are any tricks for the QP formulation that would speed up the process?

[WalterMadelim]

One solve only takes 201ms, which is mild. Is your task parallelizable? i.e. samples are independent. If so, use multithreading in julia may help.

[odow]

Yes, indeed, it is OSQP.Solve!

Great! So then my follow up question:

What makes you think 60 seconds I slow? Do you have a benchmark in a different language/solver/ecosystem?

It also looks like all of the sub-solves are independent. Have you tried forming one big problem to minimize all of the samples at the same time?