# Speed issue for Quadratic Program

**URL:** <https://discourse.julialang.org/t/speed-issue-for-quadratic-program/136898>\
**Category:** Optimization (Mathematical)\
**Created:** [April 27, 2026, 11:04am UTC](https://discourse.julialang.org/t/speed-issue-for-quadratic-program/136898 "2026-04-27T11:04:52Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![vaibhav.u](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhav.u/32/210789_2.png) [@vaibhav.u](https://discourse.julialang.org/u/vaibhav.u)\
**Post date:** [April 27, 2026, 11:04am UTC](https://discourse.julialang.org/t/speed-issue-for-quadratic-program/136898/1 "2026-04-27T11:04:52Z")

</div>

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:

```julia-auto
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

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 27, 2026, 8:45pm UTC](https://discourse.julialang.org/t/speed-issue-for-quadratic-program/136898/2 "2026-04-27T20:45:37Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![vaibhav.u](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhav.u/32/210789_2.png) [@vaibhav.u](https://discourse.julialang.org/u/vaibhav.u)\
**Post date:** [April 30, 2026, 8:06am UTC](https://discourse.julialang.org/t/speed-issue-for-quadratic-program/136898/3 "2026-04-30T08:06:02Z")

</div>

> [@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:

```julia-auto
# ============================================================
# 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:

```julia-auto
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?

---

<div class="post-metadata">

**Author:** ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)\
**Post date:** [April 30, 2026, 10:52am UTC](https://discourse.julialang.org/t/speed-issue-for-quadratic-program/136898/4 "2026-04-30T10:52:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 30, 2026, 9:16pm UTC](https://discourse.julialang.org/t/speed-issue-for-quadratic-program/136898/5 "2026-04-30T21:16:30Z")

</div>

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