Improve DAQP Prox

.gitignore

@@ -2,3 +2,4 @@ build/
 *.mexa64
 pull_requests
 *.so
+*.csv

CMakeLists.txt

@@ -108,7 +108,7 @@ if(JULIA)
             NAME JuliaBenchmark
             COMMAND bash "${JULIA_INTERFACE_DIR}/test/benchmark_comparison_git.sh"
                 "${JULIA_BINARY_DIR}"
-                "v0.7.2"
+                "v0.8.0"
                 "medium"
                 "5"
                 "${CMAKE_CURRENT_BINARY_DIR}/benchmark_results"

codegen/codegen.c

@@ -136,6 +136,8 @@ void write_daqp_workspace_h(FILE *f, DAQPWorkspace* work, const char* prefix){
 
     fprintf(f, "extern int %sWS[%d];\n\n", prefix, ntot+1);
 
+    fprintf(f, "extern int %sprox_mask[%d];\n\n", prefix, n);
+
     fprintf(f, "extern DAQPWorkspace %swork;\n\n", prefix);
 }
 
@@ -185,6 +187,8 @@ void write_daqp_workspace_src(FILE* f, DAQPWorkspace* work, const char* prefix){
 
     fprintf(f, "int %sWS[%d];\n\n", prefix, ntot+1);
 
+    fprintf(f, "int %sprox_mask[%d];\n\n", prefix, n);
+
     //Workspace struct
     fprintf(f, "DAQPWorkspace %swork= {\n", prefix);
     fprintf(f, "NULL,\n"); // DAQPProblem
@@ -199,6 +203,7 @@ void write_daqp_workspace_src(FILE* f, DAQPWorkspace* work, const char* prefix){
             prefix,prefix,prefix,prefix, 0); // reuse_ind
     fprintf(f, "%sWS, %d,\n", prefix, 0); //n_active
     fprintf(f, "%d,%d,\n",0,-1); //iterations + sing_id
+    fprintf(f, "%sprox_mask, %d,\n", prefix, n); // prox_mask
     fprintf(f, "%f,\n",0.0); // Soft slack
     fprintf(f, "&%ssettings, \n", prefix);
     // BnB

docs/docs/settings.md

@@ -26,14 +26,14 @@ Table of contents
 | `cycle_tol` | Allowed number of iterations without progress before terminating| 10 |
 | `iter_limit` | Maximum number of iterations before terminating| 10000 |
 | `fval_bound` | Maximum allowed objective function value. The solver terminates if the dual objective exceeds this value (since it is a lower bound of the optimal value). | 1e30|
-| `eps_prox` | Regularization parameter used for proximal-point iterations (0 means that no proximal-point iterations are performed) | 0|
+| `eps_prox` | Regularization parameter used for proximal-point iterations | 1e-6|
 | `eta_prox` | Tolerance that determines if a fix-point has been reached during proximal-point iterations | 1e-6|
 | `rho_soft` | Weight used for soft constraints (higher enables more violations) | 1e-6|
 | `rel_subopt` | Allowed relative suboptimality in branch and bound | 0 |
 | `abs_subopt` | Allowed absolute suboptimality in branch and bound | 0 |
 | `sing_tol` | Tolerance for checking if the LDL' factorization is singular| 3.7e-11 |
 | `refactor_tol` | Tolerance for refactoring the LDL' factorization before terminating | 1e-9 |
-| `time_limit` | Maximum allowed wall-clock time in seconds before terminating (0 means no limit) | 0 |
+| `time_limit` | Maximum wall-clock time in seconds before terminating (0 means no limit) | 0 |
 
 
 ## Exit flags 

include/constants.h

@@ -15,6 +15,7 @@ extern "C" {
 #define DAQP_DEFAULT_DUAL_TOL 1e-12
 #define DAQP_DEFAULT_ZERO_TOL 1e-11
 #define DAQP_DEFAULT_PROG_TOL 1e-14
+#define DAQP_DEFAULT_LP_PROG_TOL 1e-10
 #define DAQP_DEFAULT_PIVOT_TOL 1e-6
 #define DAQP_DEFAULT_CYCLE_TOL 10
 #define DAQP_DEFAULT_ETA 1e-6
@@ -24,6 +25,7 @@ extern "C" {
 #define DAQP_DEFAULT_ABS_SUBOPT 0
 #define DAQP_DEFAULT_SING_TOL (3.7e-11)
 #define DAQP_DEFAULT_REFACTOR_TOL 1e-9
+#define DAQP_DEFAULT_EPS_PROX 1e-6
 
 // MACROS
 #define DAQP_ARSUM(x) ((x)*(x+1)/2)

include/types.h

@@ -157,6 +157,11 @@ typedef struct{
     int iterations;
     int sing_ind; // Flag for denoting whether Mk Mk' is singular or not 
 
+    // Semi-proximal support: prox_mask[i] == 1 iff direction i needed eps
+    // regularization to make the Cholesky factor non-singular.
+    int* prox_mask;
+    int  n_prox; // Number of directions that needed regularization
+
 
     // Soft constraint
     c_float soft_slack;

interfaces/daqp-julia/src/api.jl

@@ -246,11 +246,9 @@ function setup(daqp::DAQPBase.Model, qp::DAQPBase.QPj;
     end
 
     if(isempty(qp.H) && !isempty(qp.f))# LP
-        # ensure their is no binary constraint
+        # ensure there is no binary constraint
         @assert(!any((qp.sense.&BINARY).==BINARY),
                 "DAQP requires the objective to be strictly convex to support binary variables")
-        # ensure proximal-point iterations are used for LPs
-        (old_settings.eps_prox == 0) && settings(daqp,Dict(:eps_prox=>1))
     end
 
     # Handle AVI with other setup 

interfaces/daqp-julia/src/types.jl

@@ -159,6 +159,9 @@ struct Workspace
     iterations::Cint
     sing_ind::Cint
 
+    prox_mask::Ptr{Cint}
+    n_prox::Cint
+
     soft_slack::Cdouble
 
     settings::Ptr{DAQPSettings}

interfaces/daqp-julia/test/benchmark.jl

@@ -5,15 +5,20 @@ Performance benchmark script for DAQP solver.
 Measures solve time for representative QP and LP problems at different scales.
 Results are saved to CSV for easy comparison and regression detection.
 
+Also benchmarks the semi-proximal approach vs the old full-proximal (eps*I)
+approach for QPs with rank-deficient Hessians.
+
 Usage:
     julia benchmark.jl                          # Run benchmarks with default settings
     julia benchmark.jl --output results.csv    # Specify output file
     julia benchmark.jl --suite small            # Run only small problems
+    julia benchmark.jl --prox                   # Run semi-proximal vs full-proximal comparison
 """
 
 using LinearAlgebra
 using Random
 using Statistics
+using Printf
 using DAQPBase
 using DAQP_jll
 using Dates
@@ -94,6 +99,18 @@ function benchmark_lp(n, m, ms; num_problems=10, num_repeats=5)
     setup_medians = Float64[]
     solve_medians = Float64[]
     iter_medians  = Float64[]
+    #
+    # Build an LP settings struct with eps_prox=1.  We use the old-style
+    # quadprog(QPj; settings=...) API so that eps_prox is passed directly to
+    # the C function (daqp_quadprog) rather than through the workspace struct.
+    # This makes the benchmark work correctly with both the current C library
+    # (which dispatches LPs via n_prox) and older libraries (which dispatch via
+    # eps_prox != 0), avoiding a struct-layout mismatch when swapping .so files.
+    _default_s = DAQPBase.DAQPSettings()
+    _lp_settings = DAQPBase.DAQPSettings(
+        [f == :eps_prox ? 1.0 : getfield(_default_s, f)
+         for f in fieldnames(DAQPBase.DAQPSettings)]...)
+
 
     for _ in 1:num_problems
         # Generate one test problem with known solution
@@ -104,7 +121,8 @@ function benchmark_lp(n, m, ms; num_problems=10, num_repeats=5)
         prob_iters = Int[]
 
         for _ in 1:num_repeats
-            x, fval, exitflag, info = linprog(f, A, bupper, blower, sense)
+            qpj = DAQPBase.QPj(zeros(0,0), f, A, bupper, blower, sense)
+            x, fval, exitflag, info = DAQPBase.quadprog(qpj; settings=_lp_settings)
             if abs(f' * (xref - x)) > CORRECTNESS_TOL
                 @warn "Solution quality issue: |f'*(xref - x)| = $(abs(f' * (xref - x)))"
             end
@@ -135,6 +153,159 @@ function benchmark_lp(n, m, ms; num_problems=10, num_repeats=5)
     )
 end
 
+"""
+    generate_rank_deficient_QP(n, m, ms, rank, nAct, kappa)
+
+Generate a QP whose Hessian has exactly `rank` positive eigenvalues; the
+remaining (n - rank) eigenvalues are zero.  The linear term f and constraints
+are retained from the PD reference problem so that the constrained minimum is
+well-defined (the constraints clip the null-space directions).
+
+Returns (H, f, A, bupper, blower, sense) — no `xref` since the true solution
+of the rank-deficient problem differs from the PD one.
+"""
+function generate_rank_deficient_QP(n, m, ms, rank, nAct, kappa)
+    _, H_pd, f, A, bupper, blower, sense = generate_test_QP(n, m, ms, nAct, kappa)
+    F = eigen(Symmetric(H_pd))
+    vals = max.(F.values, 0.0)      # remove floating-point noise
+    vals[rank+1:end] .= 0.0         # zero out the (n - rank) smallest eigenvalues
+    H = Symmetric(F.vectors * Diagonal(vals) * F.vectors')
+    return Matrix(H), f, A, bupper, blower, sense
+end
+
+"""
+    benchmark_prox_comparison(n, m, ms, rank, nAct, kappa; num_runs=10)
+
+Compare the semi-proximal method (new: eps applied only to singular directions,
+selected by `prox_mask`) against the old full-proximal method (eps·I applied to
+ALL variables regardless of whether they are singular).
+
+The old behaviour is simulated by adding a small `η·I` to H before factorisation
+so that all Cholesky diagonals are strictly positive and `prox_mask` is all-zero
+— meaning n_prox = 0 and daqp_ldp would run once without any outer proximal
+loop.  That single-shot solve corresponds to what the old code produced when the
+user had to supply a pre-regularised H.
+
+The difference between the methods therefore becomes visible through the *number
+of outer proximal iterations* the solver performs: the semi-proximal method
+perturbs only the truly singular directions, keeping the inner QP closer to the
+original problem, which typically requires fewer outer iterations to converge.
+"""
+function benchmark_prox_comparison(n, m, ms, rank, nAct, kappa; num_runs=10)
+    semi_times  = Float64[]
+    old_times   = Float64[]
+    semi_iters  = Int[]
+    old_iters   = Int[]
+    n_prox_vals = Int[]
+
+    for _ in 1:num_runs
+        H, f, A, bupper, blower, sense = generate_rank_deficient_QP(n, m, ms, rank, nAct, kappa)
+
+        # --- Semi-proximal (new): n_prox identifies and perturbs only singular dirs ---
+        x_semi, _, ef_semi, info_semi = quadprog(H, f, A, bupper, blower, sense)
+        push!(semi_times, info_semi.solve_time)
+        push!(semi_iters, info_semi.iterations)
+
+        # Check n_prox for this problem (all runs should be identical rank)
+        if isempty(n_prox_vals)
+            d_tmp = DAQPBase.Model()
+            setup(d_tmp, H, f, A, bupper, blower, sense)
+            ws_tmp = unsafe_load(Ptr{DAQPBase.Workspace}(d_tmp.work))
+            push!(n_prox_vals, ws_tmp.n_prox)
+        end
+
+        # --- Old full-proximal (simulated): pre-regularise H with a small η·I
+        # so all directions look PD to the Cholesky.  This mimics the old code
+        # which applied eps to every direction unconditionally.  The resulting
+        # problem is then solved with eps_prox = 0 (direct daqp_ldp), giving the
+        # "one-shot" result the old code would have produced. ---
+        η = 1.0  # regularisation amount added to H (matches the eps used in the old full-proximal approach)
+        H_reg = H + η * I
+        s_old = settings(DAQPBase.Model(), Dict(:eps_prox => 0.0))
+        x_old, _, ef_old, info_old = quadprog(H_reg, f, A, bupper, blower, sense; settings=s_old)
+        push!(old_times, info_old.solve_time)
+        push!(old_iters, info_old.iterations)
+    end
+
+    return (
+        semi_time_mean = mean(semi_times),
+        semi_time_std  = std(semi_times),
+        semi_iter_mean = mean(semi_iters),
+        semi_iter_std  = std(semi_iters),
+        old_time_mean  = mean(old_times),
+        old_time_std   = std(old_times),
+        old_iter_mean  = mean(old_iters),
+        old_iter_std   = std(old_iters),
+        n_prox         = isempty(n_prox_vals) ? -1 : n_prox_vals[1],
+        num_runs       = num_runs,
+    )
+end
+
+"""
+    run_prox_benchmark(; output_file="prox_comparison.csv", use_local=false)
+
+Run semi-proximal vs old full-proximal comparison and print a summary table.
+
+Columns:
+  n          — problem dimension
+  rank       — number of positive Hessian eigenvalues
+  n_sing     — number of singular directions (n - rank)
+  n_prox     — directions actually regularised by semi-proximal
+  semi_µs    — semi-proximal solve time (µs)
+  old_µs     — old full-proximal solve time (µs)
+  semi_iters — outer iterations for semi-proximal
+  old_iters  — outer iterations for old approach (always 1 since it's a
+               direct single solve of the pre-regularised problem)
+"""
+function run_prox_benchmark(; output_file="prox_comparison.csv", use_local=false)
+    _libdaqp = joinpath(pkgdir(DAQPBase),"libdaqp."*Libc.Libdl.dlext)
+    if isfile(_libdaqp) && use_local
+        @info "Using local libdaqp"
+        DAQP_jll.libdaqp = _libdaqp
+    end
+
+    # (n, m, ms, rank, nAct, kappa)
+    cases = [
+        (20,  100,  10,  15,  16,  1e2),   # small,  1 singular direction
+        (20,  100,  10,  10,  16,  1e2),   # small,  half-rank
+        (20,  100,  10,   5,  16,  1e2),   # small,  quarter-rank
+        (50,  250,  25,  40,  40,  1e2),   # medium, 10 singular directions
+        (50,  250,  25,  25,  40,  1e2),   # medium, half-rank
+        (100, 500,  50,  80,  80,  1e2),   # large,  20 singular directions
+    ]
+
+    println("\n=== Semi-proximal vs Old Full-proximal Benchmark ===")
+    println("  semi: new code — eps only on singular directions (n_prox ≤ n)")
+    println("  old:  old code — pre-regularised H with eps·I on all directions")
+    println()
+    header = @sprintf("%-6s %-5s %-6s %-7s  %10s %10s  %11s %11s",
+                      "n", "rank", "n_sing", "n_prox",
+                      "semi_µs", "old_µs",
+                      "semi_iters", "old_iters")
+    println(header)
+    println("-"^length(header))
+
+    csv_rows = String["n,rank,n_sing,n_prox,semi_time_mean_us,semi_time_std_us,old_time_mean_us,old_time_std_us,semi_iter_mean,old_iter_mean,num_runs"]
+
+    for (n, m, ms, rank, nAct, kappa) in cases
+        r = benchmark_prox_comparison(n, m, ms, rank, nAct, kappa; num_runs=10)
+        n_sing = n - rank
+        println(@sprintf("%-6d %-5d %-6d %-7d  %10.1f %10.1f  %11.1f %11.1f",
+                         n, rank, n_sing, r.n_prox,
+                         r.semi_time_mean*1e6, r.old_time_mean*1e6,
+                         r.semi_iter_mean, r.old_iter_mean))
+        push!(csv_rows, join([n, rank, n_sing, r.n_prox,
+                              r.semi_time_mean*1e6, r.semi_time_std*1e6,
+                              r.old_time_mean*1e6,  r.old_time_std*1e6,
+                              r.semi_iter_mean, r.old_iter_mean, r.num_runs], ","))
+    end
+
+    open(output_file, "w") do io
+        foreach(l -> println(io, l), csv_rows)
+    end
+    println("\n✓ Results saved to: $output_file")
+end
+
 function run_benchmarks(; suite="all", output_file="daqp_benchmark_results.csv", use_local=false)
     """
     Run all benchmarks and save results to CSV.
@@ -252,6 +423,7 @@ if abspath(PROGRAM_FILE) == @__FILE__
     local suite = "all"
     local output_file = "daqp_benchmark_results.csv"
     local use_local = false
+    local run_prox = false
     local i = 1
     while i <= length(ARGS)
         if ARGS[i] == "--output" && i < length(ARGS)
@@ -263,10 +435,17 @@ if abspath(PROGRAM_FILE) == @__FILE__
         elseif ARGS[i] == "--local"
             use_local = true 
             i += 1 
+        elseif ARGS[i] == "--prox"
+            run_prox = true
+            i += 1
         else
             i += 1
         end
     end
 
-    run_benchmarks(; suite=suite, output_file=output_file, use_local=use_local)
+    if run_prox
+        run_prox_benchmark(; output_file="prox_comparison.csv", use_local=use_local)
+    else
+        run_benchmarks(; suite=suite, output_file=output_file, use_local=use_local)
+    end
 end