diff --git a/.gitignore b/.gitignore
index 2d24bf50..71cf35e6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@ build/
 *.mexa64
 pull_requests
 *.so
+*.csv
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 42d71efe..40a9ee42 100644
--- a/CMakeLists.txt
+++ b/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"
diff --git a/codegen/codegen.c b/codegen/codegen.c
index f326acef..2cfa4d1a 100644
--- a/codegen/codegen.c
+++ b/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
diff --git a/docs/docs/settings.md b/docs/docs/settings.md
index 2c48590b..d5bc7cd6 100644
--- a/docs/docs/settings.md
+++ b/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 
diff --git a/include/constants.h b/include/constants.h
index f5bb40aa..2dfd4e43 100644
--- a/include/constants.h
+++ b/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)
diff --git a/include/types.h b/include/types.h
index 8ffaf6df..dd731141 100644
--- a/include/types.h
+++ b/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;
diff --git a/interfaces/daqp-julia/src/api.jl b/interfaces/daqp-julia/src/api.jl
index 07034c14..a65610a1 100644
--- a/interfaces/daqp-julia/src/api.jl
+++ b/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 
diff --git a/interfaces/daqp-julia/src/types.jl b/interfaces/daqp-julia/src/types.jl
index f70532c2..ccac72fa 100644
--- a/interfaces/daqp-julia/src/types.jl
+++ b/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}
diff --git a/interfaces/daqp-julia/test/benchmark.jl b/interfaces/daqp-julia/test/benchmark.jl
index 9f172a36..4a7bb082 100644
--- a/interfaces/daqp-julia/test/benchmark.jl
+++ b/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
diff --git a/interfaces/daqp-julia/test/core_tests.jl b/interfaces/daqp-julia/test/core_tests.jl
index 0452d7ac..3a71388b 100644
--- a/interfaces/daqp-julia/test/core_tests.jl
+++ b/interfaces/daqp-julia/test/core_tests.jl
@@ -431,3 +431,66 @@ end
     @test norm(xref-x) < tol
 end
 
+@testset "Semi-proximal method" begin
+    # --- PD Hessian: n_prox must be 0, result matches eps_prox=0 solve ---
+    n2 = 5; m2 = 20; ms2 = 5; nAct2 = 4
+    xref,H,f,A,bupper,blower,sense = generate_test_QP(n2,m2,ms2,nAct2,1e2)
+
+    # Reference (no proximal)
+    xref2,fval_ref,ef_ref,_ = quadprog(H,f,A,bupper,blower,sense)
+    @test ef_ref == DAQPBase.OPTIMAL
+
+    # With eps_prox > 0 on PD H: solver should still reach the optimum.
+    s_prox = settings(DAQPBase.Model(), Dict(:eps_prox => 1e-4))
+    x_prox,_,ef_prox,_ = quadprog(H,f,A,bupper,blower,sense; settings=s_prox)
+    @test ef_prox == DAQPBase.OPTIMAL
+    @test norm(xref2 - x_prox) < tol
+
+    # Check that the Workspace struct layout is correct by reading n_prox via
+    # unsafe_load -- it must be 0 for a PD Hessian (no regularisation needed).
+    d = DAQPBase.Model()
+    DAQPBase.settings(d, Dict(:eps_prox => 1e-4))
+    DAQPBase.setup(d, H, f, A, bupper, blower, sense)
+    p = d.work  # raw workspace pointer (Ptr{Cvoid})
+    ws = unsafe_load(Ptr{DAQPBase.Workspace}(p))
+    @test ws.n_prox == 0
+
+    # --- Rank-1 Hessian: x2 direction is singular -> n_prox == 1 ---
+    H_sing = [1.0 0.0; 0.0 0.0]
+    f_sing = [1.0; 1.0]
+    A_sing = zeros(0, 2)
+    bu_sing = [2.0; 2.0]
+    bl_sing = [-2.0; -2.0]
+    sense_sing = zeros(Cint, 2)
+
+    d2 = DAQPBase.Model()
+    DAQPBase.settings(d2, Dict(:eps_prox => 1e-3))
+    DAQPBase.setup(d2, H_sing, f_sing, A_sing, bu_sing, bl_sing, sense_sing)
+    ws2 = unsafe_load(Ptr{DAQPBase.Workspace}(d2.work))
+    @test ws2.n_prox == 1   # only x2 direction needed regularisation
+
+    x2,_,ef2,_ = DAQPBase.solve(d2)
+    @test ef2 == DAQPBase.OPTIMAL
+    @test abs(x2[1] - (-1.0)) < tol   # x1* = -1
+    @test abs(x2[2] - (-2.0)) < tol   # x2* = -2 (lower bound)
+
+    # --- Zero Hessian: all directions singular -> n_prox == n ---
+    n3 = 3
+    H_zero = zeros(n3, n3)
+    f_zero = ones(n3)
+    A_zero = zeros(0, n3)
+    bu_zero =  5.0*ones(n3)
+    bl_zero = -5.0*ones(n3)
+    sense_zero = zeros(Cint, n3)
+
+    d3 = DAQPBase.Model()
+    DAQPBase.settings(d3, Dict(:eps_prox => 1e-2))
+    DAQPBase.setup(d3, H_zero, f_zero, A_zero, bu_zero, bl_zero, sense_zero)
+    ws3 = unsafe_load(Ptr{DAQPBase.Workspace}(d3.work))
+    @test ws3.n_prox == n3
+
+    x3,_,ef3,_ = DAQPBase.solve(d3)
+    @test ef3 > 0
+    @test norm(x3 .- (-5.0)) < tol  # all at lower bound
+end
+
diff --git a/interfaces/daqp-python/daqp.pyx b/interfaces/daqp-python/daqp.pyx
index 3906d5d2..8c849beb 100644
--- a/interfaces/daqp-python/daqp.pyx
+++ b/interfaces/daqp-python/daqp.pyx
@@ -367,12 +367,20 @@ cdef class Model:
             free_daqp_workspace(self._work)   # also frees settings
             free_daqp_ldp(self._work)
         else:
-            # Free the settings allocated in __cinit__ so setup_daqp can
-            # allocate its own (preventing the 'own_settings=0' leak path).
             if self._work.settings != NULL:
                 free(self._work.settings)
                 self._work.settings = NULL
 
+        # ---- pre-allocate settings with user values so setup_daqp sees them ----
+        # This is critical: daqp_update_Rinv (called inside setup_daqp) reads
+        # eps_prox from settings to decide which Cholesky diagonals need
+        # regularisation.  If settings were NULL or held default (eps_prox=0)
+        # when the Cholesky runs, a singular Hessian would be rejected as
+        # non-convex even when the caller supplied a non-zero eps_prox.
+        allocate_daqp_settings(self._work)   # fresh allocation, defaults
+        if restore_settings:
+            self._work.settings[0] = old_settings_val  # apply user values NOW
+
         # ---- build the problem struct ----
         cdef double* H_ptr  = NULL if H is None else &self._H[0, 0]
         cdef double* f_ptr  = NULL if f is None else &self._f[0]
@@ -405,14 +413,7 @@ cdef class Model:
                 with nogil:
                     daqp_primal_init_active(&self._qp, &primal_start_c[0])
 
-        # ---- enable proximal iterations for pure LPs ----
-        # (mirrors the Julia interface behaviour)
-        if H is None and f is not None:
-            # settings may be NULL at this point; setup_daqp will allocate them
-            pass  # eps_prox handled after successful setup below
-
         # ---- call setup_daqp ----
-        # work->settings is NULL here; setup_daqp will allocate (own_settings=1)
         with nogil:
             setup_flag = setup_daqp(&self._qp, self._work, &setup_time_c)
 
@@ -428,11 +429,6 @@ cdef class Model:
         if restore_settings:
             self._work.settings[0] = old_settings_val
 
-        # ---- LP: ensure proximal-point iterations are active ----
-        if H is None and f is not None:
-            if self._work.settings.eps_prox == 0:
-                self._work.settings.eps_prox = 1
-
         self._has_model = True
 
         # ---- set primal iterate ----
diff --git a/interfaces/daqp-python/test/example_test.py b/interfaces/daqp-python/test/example_test.py
index d294d223..e786ff5a 100644
--- a/interfaces/daqp-python/test/example_test.py
+++ b/interfaces/daqp-python/test/example_test.py
@@ -277,5 +277,88 @@ def test_model_setup_does_not_mutate_sense(self):
                                       err_msg="setup primal_start must not mutate sense")
 
 
+class TestSemiProximal(unittest.TestCase):
+    """Tests for the semi-proximal method (eps_prox > 0)."""
+
+    def test_pd_hessian_n_prox_zero(self):
+        """PD Hessian with eps_prox > 0: no direction needs regularisation (n_prox=0)."""
+        H = np.array([[1.0, 0.0], [0.0, 1.0]], dtype=c_double)
+        f = np.array([1.0, 1.0], dtype=c_double)
+        A = np.zeros((0, 2), dtype=c_double)
+        bupper = np.array([1.0, 1.0], dtype=c_double)
+        blower = np.array([-1.0, -1.0], dtype=c_double)
+        sense = np.array([0, 0], dtype=c_int)
+
+        d = daqp.Model()
+        d.settings = {'eps_prox': 1e-4}
+        d.setup(H, f, A, bupper, blower, sense)
+        x, fval, ef, info = d.solve()
+
+        self.assertEqual(ef, 1)
+        np.testing.assert_allclose(x, [-1.0, -1.0], atol=1e-4)
+
+    def test_singular_hessian_semi_proximal(self):
+        """Rank-1 Hessian: only the singular direction gets regularised."""
+        # H = diag(1, 0): x2 direction is singular
+        H = np.array([[1.0, 0.0], [0.0, 0.0]], dtype=c_double)
+        f = np.array([1.0, 1.0], dtype=c_double)
+        # Use simple-bound columns as A rows so A is non-empty (ms=0, mA=2)
+        A = np.eye(2, dtype=c_double)
+        bupper = np.array([2.0, 2.0], dtype=c_double)
+        blower = np.array([-2.0, -2.0], dtype=c_double)
+        sense = np.array([0, 0], dtype=c_int)
+
+        d = daqp.Model()
+        d.settings = {'eps_prox': 1e-3}
+        d.setup(H, f, A, bupper, blower, sense)
+        x, fval, ef, info = d.solve()
+
+        self.assertEqual(ef, 1)
+        # x1 minimises 0.5*x1^2 + x1 unconstrained -> x1* = -1
+        # x2 is singular; regularisation + f[1]=1 pushes to lower bound
+        self.assertAlmostEqual(x[0], -1.0, places=3)
+        self.assertAlmostEqual(x[1], -2.0, places=3)
+
+    def test_zero_hessian_all_singular(self):
+        """Zero Hessian: all directions are singular, full proximal applied."""
+        H = np.array([[0.0, 0.0], [0.0, 0.0]], dtype=c_double)
+        f = np.array([1.0, 2.0], dtype=c_double)
+        # Use simple-bound columns as A rows so A is non-empty
+        A = np.eye(2, dtype=c_double)
+        bupper = np.array([3.0, 3.0], dtype=c_double)
+        blower = np.array([-3.0, -3.0], dtype=c_double)
+        sense = np.array([0, 0], dtype=c_int)
+
+        d = daqp.Model()
+        d.settings = {'eps_prox': 1e-2}
+        d.setup(H, f, A, bupper, blower, sense)
+        x, fval, ef, info = d.solve()
+
+        self.assertGreater(ef, 0)
+        # Regularised problem: min eps/2*||x||^2 + f'*x  -> x* = -f/eps (clipped to bounds)
+        np.testing.assert_allclose(x, [-3.0, -3.0], atol=1e-3)
+
+    def test_pd_hessian_matches_no_prox(self):
+        """With a PD Hessian, eps_prox > 0 still gives the correct solution."""
+        H = np.array([[4.0, 1.0], [1.0, 3.0]], dtype=c_double)
+        f = np.array([1.0, 2.0], dtype=c_double)
+        A = np.zeros((0, 2), dtype=c_double)
+        bupper = np.array([5.0, 5.0], dtype=c_double)
+        blower = np.array([-5.0, -5.0], dtype=c_double)
+        sense = np.array([0, 0], dtype=c_int)
+
+        # Reference: solve without proximal
+        x_ref, _, ef_ref, _ = daqp.solve(H, f, A, bupper, blower, sense)
+        self.assertEqual(ef_ref, 1)
+
+        # Solve with proximal on PD H -- should converge to same solution
+        d = daqp.Model()
+        d.settings = {'eps_prox': 1e-4}
+        d.setup(H, f, A, bupper, blower, sense)
+        x_prox, _, ef_prox, _ = d.solve()
+        self.assertEqual(ef_prox, 1)
+        np.testing.assert_allclose(x_prox, x_ref, atol=1e-3)
+
+
 if __name__ == '__main__':
     unittest.main()
diff --git a/src/api.c b/src/api.c
index 986dc49d..8d00ddfe 100644
--- a/src/api.c
+++ b/src/api.c
@@ -12,7 +12,7 @@ void daqp_solve(DAQPResult *res, DAQPWorkspace *work){
     if(work->settings->time_limit > 0)  work->timer = &timer;
 #endif
     // Select algorithm
-    if(work->settings->eps_prox==0){
+    if(work->n_prox==0){
         if(work->avi == NULL){
             if(work->bnb != NULL)
                 res->exitflag = daqp_bnb(work);
@@ -150,8 +150,8 @@ int setup_daqp_ldp(DAQPWorkspace *work, DAQPProblem *qp){
         alloc_R = 1;
         update_mask+=DAQP_UPDATE_Rinv;
     }
-    // Only allocate v if f is not NULL, or if proximal
-    if(qp->f!=NULL || work->settings->eps_prox != 0){
+    // Only allocate v if f is not NULL (QP linear term or LP objective).
+    if(qp->f!=NULL){
         alloc_v = 1;
         update_mask+=DAQP_UPDATE_v;
     }
@@ -168,6 +168,10 @@ int setup_daqp_ldp(DAQPWorkspace *work, DAQPProblem *qp){
         free_daqp_ldp(work);
         return error_flag;
     }
+    // For LPs (no Hessian), mark all directions as needing proximal regularisation.
+    // This lets daqp_solve dispatch to daqp_prox based on n_prox rather than eps_prox.
+    if(qp->H == NULL && qp->f!=NULL)
+        work->n_prox = work->n;
     return 1;
 }
 
@@ -278,6 +282,9 @@ void allocate_daqp_workspace(DAQPWorkspace *work, int n, int ns){
 
     work->xold= malloc(work->n*sizeof(c_float));
 
+    work->prox_mask = calloc(work->n, sizeof(int)); // all zeros initially
+    work->n_prox = 0;
+
 #ifdef SOFT_WEIGHTS
     work->d_ls= NULL;
     work->d_us= NULL;
@@ -362,6 +369,8 @@ void free_daqp_workspace(DAQPWorkspace *work){
 
         free(work->xold);
 
+        free(work->prox_mask);
+
         work->lam = NULL;
     }
 
@@ -413,14 +422,15 @@ void daqp_extract_result(DAQPResult* res, DAQPWorkspace* work){
     }
 
     // Shift back function value
-    if(work->v != NULL && work->avi == NULL && (work->settings->eps_prox == 0
-                || work->Rinv != NULL || work->RinvD != NULL)){ // Normal QP
+    if(work->v != NULL && work->avi == NULL && (work->Rinv != NULL || work->RinvD != NULL)){ // QP
         res->fval = work->fval;
         for(i=0;i<work->n;i++) res->fval-=work->v[i]*work->v[i];
         res->fval *=0.5;
-        if(work->settings->eps_prox != 0)
-            for(i=0;i<work->n;i++) // compensate for proximal iterations
-                res->fval+= work->settings->eps_prox*work->x[i]*work->x[i];
+        if(work->n_prox > 0)
+            for(i=0;i<work->n;i++) // remove proximal bias: at fixed point x≈x_old so
+                // true_fval = perturbed_fval + 0.5*eps*||x_mask||^2
+                if(work->prox_mask == NULL || work->prox_mask[i])
+                    res->fval+= 0.5*work->settings->eps_prox*work->x[i]*work->x[i];
     }
     else if(work->qp != NULL && work->qp->f != NULL ){ // LP
         res->fval = 0;
@@ -444,7 +454,7 @@ void daqp_default_settings(DAQPSettings* settings){
     settings->iter_limit = DAQP_DEFAULT_ITER_LIMIT;
     settings->fval_bound = DAQP_INF; 
 
-    settings->eps_prox = 0;
+    settings->eps_prox = DAQP_DEFAULT_EPS_PROX;
     settings->eta_prox = DAQP_DEFAULT_ETA;
 
     settings->rho_soft = DAQP_DEFAULT_RHO_SOFT; 
diff --git a/src/daqp_prox.c b/src/daqp_prox.c
index 2e272d3d..bc58b5b7 100644
--- a/src/daqp_prox.c
+++ b/src/daqp_prox.c
@@ -3,151 +3,242 @@
 
 static int gradient_step(DAQPWorkspace* work);
 
+/* --------------------------------------------------------------------------
+ * daqp_prox  --  outer proximal-point / semi-proximal loop
+ *
+ * QP problems (Rinv or RinvD is set):
+ *   Implements a Semi-Proximal Method.  Only the directions i where the
+ *   Cholesky diagonal would have been non-positive without regularisation
+ *   (prox_mask[i] == 1) receive the eps·x_old[i] perturbation in the
+ *   right-hand side.  If H is already positive definite everywhere
+ *   (n_prox == 0) the inner QP equals the original and we exit after one
+ *   solve.
+ *
+ * LP problems (Rinv == NULL && RinvD == NULL):
+ *   Classical regularisation-based smoothing with adaptive eps.
+ * --------------------------------------------------------------------------*/
 int daqp_prox(DAQPWorkspace *work){
-    int i,total_iter=0;
-    const int nx=work->n;
+    int i, total_iter = 0;
+    const int nx = work->n;
     int exitflag;
     c_float *swp_ptr;
-    c_float diff,eps=work->settings->eps_prox;
-    c_float tol_stat, eta=work->settings->eta_prox;
-    int cycle_counter = 0;
-    c_float best_fval = DAQP_INF;
-
-    while(total_iter  <  work->settings->iter_limit){
-        // ** Perturb problem **
-        // Compute v = R'\(f-eps*x) (FWS Skipped if LP since R = I)
-        if(work->Rinv== NULL && work->RinvD == NULL){
-            eps*= (work->iterations==1) ? 10 : 0.9; // Adapt epsilon TODO: add to settings
-            eps = (eps > 1e3) ? 1e3 : eps; // TODO: add saturation option to settings
-            for(i = 0; i<nx;i++)
-                work->v[i] = work->qp->f[i]*eps-work->x[i];
+    c_float max_diff, tol_stat;
+    c_float eta = work->settings->eta_prox;
+    int  cycle_counter = 0;
+    c_float best_fval  = DAQP_INF;
+
+    const int is_lp = (work->Rinv == NULL && work->RinvD == NULL);
+    c_float eps = is_lp ? 1.0 : work->settings->eps_prox;
+
+    // For a QP whose Hessian is already positive definite (n_prox == 0),
+    // no direction needs a proximal shift.  The inner QP equals the
+    // original problem, so one solve gives the exact solution.
+    const int all_pd = (!is_lp) && (work->n_prox == 0);
+
+    while(total_iter < work->settings->iter_limit){
+
+        /* ----------------------------------------------------------------
+         * Perturb the problem: form v = R'\(f - eps_mask * x_old)
+         * ----------------------------------------------------------------*/
+        if(is_lp){
+            // No Hessian factor.  Adapt eps heuristically: grow when the
+            // inner LP stalls (iterations==1), shrink otherwise to improve
+            // accuracy.
+            eps *= (work->iterations == 1) ? 10.0 : 0.9;
+            if(eps > 1e3) eps = 1e3;
+            for(i = 0; i < nx; i++)
+                work->v[i] = work->qp->f[i]*eps - work->x[i];
         }
         else{
-            for(i = 0; i<nx;i++)
-                work->v[i] = work->qp->f[i]-eps*work->x[i];
-            daqp_update_v(work->v,work,0);
+            // Semi-proximal: shift f only for directions that needed
+            // regularisation to make the Cholesky factor non-singular.
+            if(work->prox_mask != NULL){
+                for(i = 0; i < nx; i++)
+                    work->v[i] = work->qp->f[i]
+                                 - (work->prox_mask[i] ? eps : 0.0) * work->x[i];
+            }
+            else{
+                // prox_mask unavailable -- fall back to full proximal
+                for(i = 0; i < nx; i++)
+                    work->v[i] = work->qp->f[i] - eps * work->x[i];
+            }
+            daqp_update_v(work->v, work, 0);
         }
-        // Perturb RHS of constraints
-        daqp_update_d(work, work->qp->bupper,work->qp->blower);
 
-        // xold <-- x
+        // Perturb RHS of constraints to match the shifted objective
+        daqp_update_d(work, work->qp->bupper, work->qp->blower);
+
+        // xold <-- x  (pointer swap avoids copying)
         swp_ptr = work->xold; work->xold = work->x; work->x = swp_ptr;
 
-        // ** Solve least-distance problem **
+        /* ----------------------------------------------------------------
+         * Solve the (regularised) least-distance problem
+         * ----------------------------------------------------------------*/
         work->u = work->x;
         exitflag = daqp_ldp(work);
 
         total_iter += work->iterations;
-        if(exitflag<0) 
-            return exitflag; // Could not solve LDP -> return
-        else 
-            ldp2qp_solution(work); // Get qp solution 
-
-        if(eps==0) break; // No regularization -> no outer iterations
-
-        // ** Check convergence **
-        if(work->iterations==1){ // No changes to the working set 
-            tol_stat = (work->Rinv == NULL && work->RinvD == NULL) ? eta*eps : eta/eps;
-            for(i=0, diff= 0;i<nx;i++){ // ||x_old - x|| > eta  ?
-                diff= work->x[i] - work->xold[i];
-                if((diff> tol_stat) || (diff< -tol_stat)) break;
+        if(exitflag < 0)
+            break;              // Inner solver failed -- propagate error
+        else
+            ldp2qp_solution(work); // Recover QP primal from LDP dual
+
+        if(eps == 0) break;     // No regularisation -> single outer step
+
+        /* ----------------------------------------------------------------
+         * If H is fully positive definite, the inner QP is the original
+         * problem.  The first solve gives the exact solution.
+         * ----------------------------------------------------------------*/
+        if(all_pd){
+            exitflag = DAQP_EXIT_OPTIMAL;
+            break;
+        }
+
+        /* ----------------------------------------------------------------
+         * Convergence check: fixed point  ||x - x_old||_inf < tol_stat
+         * Only test when the active set did not change (iterations == 1),
+         * because that is the cheapest indicator that a stationary point
+         * has been reached for the inner QP.
+         * ----------------------------------------------------------------*/
+        if(work->iterations == 1){
+            tol_stat = is_lp ? eta*eps : eta/eps;
+            for(i = 0; i < nx; i++){
+                max_diff = work->x[i] - work->xold[i];
+                if(max_diff > tol_stat || max_diff < -tol_stat) break;
             }
-            if(i==nx){
-                exitflag = DAQP_EXIT_OPTIMAL; // Fix point reached
+            if(i == nx){
+                exitflag = DAQP_EXIT_OPTIMAL; // Fixed point reached
                 break;
             }
-            // Take gradient step if LP (and the iterate is not constrained to a vertex)
-            if((work->Rinv == NULL && work->RinvD ==NULL )&&(work->n_active != work->n)){
-                if(gradient_step(work)==DAQP_EMPTY_IND){
-                    exitflag= DAQP_EXIT_UNBOUNDED;
+            // LP: when not at a vertex take a gradient step toward the
+            // nearest constraint to escape from the interior.
+            if(is_lp && work->n_active != nx){
+                if(gradient_step(work) == DAQP_EMPTY_IND){
+                    exitflag = DAQP_EXIT_UNBOUNDED;
+                    break;
+                }
+            }
+        }
+
+        /* ----------------------------------------------------------------
+         * Stagnation / cycle detection
+         * ----------------------------------------------------------------*/
+        if(is_lp){
+            // Track LP objective f'x; no improvement over several
+            // consecutive iterations signals a fixed point.
+            // Use DAQP_DEFAULT_LP_PROG_TOL (1e-10) rather than the QP
+            // progress_tol (1e-14): the LP objective is in raw problem
+            // units and a much larger threshold is appropriate.
+            c_float lp_obj = 0.0;
+            for(i = 0; i < nx; i++) lp_obj += work->qp->f[i] * work->x[i];
+            max_diff = best_fval - lp_obj;
+            if(max_diff < DAQP_DEFAULT_LP_PROG_TOL){
+                if(cycle_counter++ > work->settings->cycle_tol){
+                    exitflag = DAQP_EXIT_OPTIMAL; // Stagnated at fixed point
                     break;
                 }
             }
+            else{
+                best_fval    = lp_obj;
+                cycle_counter = 0;
+            }
         }
-        // For LPs, compute objective function value to detect progress
-        // TODO: add paramters to settings
-        if(work->Rinv == NULL && work->RinvD == NULL){
-            for(i=0, diff=best_fval;i<nx;i++)
-                diff-=work->qp->f[i]*work->x[i];
-            if(diff<1e-10){
-                if(cycle_counter++ > 10) return DAQP_EXIT_OPTIMAL; // assume fix point
+        else{
+            // QP: track the inner LDP dual objective ||u||^2 (work->fval).
+            // When the outer iterate converges, v stabilises and so does
+            // work->fval; stagnation therefore indicates a fixed point.
+            max_diff = best_fval - work->fval;
+            if(max_diff < work->settings->progress_tol){
+                if(++cycle_counter > work->settings->cycle_tol){
+                    exitflag = DAQP_EXIT_OPTIMAL; // Stagnated at fixed point
+                    break;
+                }
             }
-            else{ // Progress -> update objective function value
-                best_fval -=diff ;
-                cycle_counter=0;
+            else{
+                best_fval    = work->fval;
+                cycle_counter = 0;
             }
         }
-
     }
-    // Finalize results
-    if(total_iter >= work->settings->iter_limit) exitflag = DAQP_EXIT_ITERLIMIT; 
-    if(work->Rinv == NULL && work->RinvD == NULL){
-        for(i = 0; i<work->n_active;i++) 
-            work->lam_star[i]/=eps;// Rescale dual variables
+
+    // Finalize
+    if(total_iter >= work->settings->iter_limit) exitflag = DAQP_EXIT_ITERLIMIT;
+    if(is_lp){
+        for(i = 0; i < work->n_active; i++)
+            work->lam_star[i] /= eps; // Rescale dual variables
     }
     work->iterations = total_iter;
     return exitflag;
 }
-// Gradient step
-// TODO: could probably reuse code from daqp
+
+/* --------------------------------------------------------------------------
+ * gradient_step  --  line-search toward the first blocking constraint
+ *
+ * Used for LP problems when the current iterate is not at a vertex (the
+ * active set does not span all n variables).  Advances x along the
+ * direction delta_x = x - x_old until the nearest constraint boundary is
+ * hit, then returns the index of that constraint.
+ * Returns DAQP_EMPTY_IND if the problem is unbounded in that direction.
+ * --------------------------------------------------------------------------*/
 static int gradient_step(DAQPWorkspace* work){
-    int j,k,disp,add_ind=DAQP_EMPTY_IND;
-    const int nx=work->n;
-    const int m=work->m;
-    const int ms=work->ms;
-    c_float Ax,delta_s, min_alpha=DAQP_INF;
-    // Find constraint j to add: j =  argmin_j s_j
-    // Simple bounds
-    for(j=0, disp=0;j<ms;j++){
-        if(work->sense[j]&(DAQP_ACTIVE+DAQP_IMMUTABLE)) continue;
-        delta_s = work->x[j]-work->xold[j];
-        if(delta_s>0 && //Feasible descent direction
-                work->qp->bupper[j]<DAQP_INF && // Not single-sided
-                work->qp->bupper[j]-work->x[j]<min_alpha*delta_s){
-            add_ind = j;
-            min_alpha = (work->qp->bupper[j]-work->x[j])/delta_s;
+    int j, k, disp, add_ind = DAQP_EMPTY_IND;
+    const int nx = work->n;
+    const int m  = work->m;
+    const int ms = work->ms;
+    c_float Ax, delta_s, min_alpha = DAQP_INF;
+
+    // Simple bounds: find first blocking constraint along delta_x
+    for(j = 0; j < ms; j++){
+        if(work->sense[j] & (DAQP_ACTIVE + DAQP_IMMUTABLE)) continue;
+        delta_s = work->x[j] - work->xold[j];
+        if(delta_s > 0 &&
+                work->qp->bupper[j] < DAQP_INF &&
+                work->qp->bupper[j] - work->x[j] < min_alpha*delta_s){
+            add_ind   = j;
+            min_alpha = (work->qp->bupper[j] - work->x[j]) / delta_s;
         }
-        else if(delta_s < 0 && //Feasible descent direction
-                work->qp->blower[j]>-DAQP_INF && // Not single-sided
-                work->qp->blower[j]-work->x[j]>min_alpha*delta_s){
-            add_ind = j;
-            min_alpha = (work->qp->blower[j]-work->x[j])/delta_s;
+        else if(delta_s < 0 &&
+                work->qp->blower[j] > -DAQP_INF &&
+                work->qp->blower[j] - work->x[j] > min_alpha*delta_s){
+            add_ind   = j;
+            min_alpha = (work->qp->blower[j] - work->x[j]) / delta_s;
         }
     }
-    //General bounds
-    for(j=ms, disp=0;j<m;j++){
-        if(work->sense[j]&(DAQP_ACTIVE+DAQP_IMMUTABLE)){
-            disp+=nx;// Skip ahead in A
+
+    // General bounds
+    for(j = ms, disp = 0; j < m; j++){
+        if(work->sense[j] & (DAQP_ACTIVE + DAQP_IMMUTABLE)){
+            disp += nx;
             continue;
         }
-        //delta_s[j] = A[j,:]*delta_x
-        for(k=0,delta_s=0,Ax=0;k<nx;k++){ // compute s = A(x-xold) and Ax
-            Ax += work->M[disp]*work->x[k];
-            delta_s-=work->M[disp++]*work->xold[k];
+        for(k = 0, delta_s = 0, Ax = 0; k < nx; k++){
+            Ax      += work->M[disp]   * work->x[k];
+            delta_s -= work->M[disp++] * work->xold[k];
         }
-        delta_s +=Ax;
+        delta_s += Ax;
 
         if(work->scaling != NULL){
-            Ax /= work->scaling[j];
+            Ax      /= work->scaling[j];
             delta_s /= work->scaling[j];
         }
-        if(delta_s>0 && // Feasible descent direction
-                work->qp->bupper[j]<DAQP_INF && // Not single-sided
-                work->qp->bupper[j]-Ax < delta_s*min_alpha){
-            add_ind = j;
-            min_alpha=(work->qp->bupper[j]-Ax)/delta_s;
+        if(delta_s > 0 &&
+                work->qp->bupper[j] < DAQP_INF &&
+                work->qp->bupper[j] - Ax < delta_s*min_alpha){
+            add_ind   = j;
+            min_alpha = (work->qp->bupper[j] - Ax) / delta_s;
         }
-        else if(delta_s<0 && // Feasible descent direction
-                work->qp->blower[j]>-DAQP_INF && // Not single-sided
-                work->qp->blower[j]-Ax > delta_s*min_alpha){
-            add_ind = j;
-            min_alpha=(work->qp->blower[j]-Ax)/delta_s;
+        else if(delta_s < 0 &&
+                work->qp->blower[j] > -DAQP_INF &&
+                work->qp->blower[j] - Ax > delta_s*min_alpha){
+            add_ind   = j;
+            min_alpha = (work->qp->blower[j] - Ax) / delta_s;
         }
     }
-    // update iterate
+
+    // Advance: x <-- x + min_alpha * (x - x_old)
     if(add_ind != DAQP_EMPTY_IND)
-        for(k=0;k<nx;k++) // x <-- x+alpha deltax
-            work->x[k]+=min_alpha*(work->x[k]-work->xold[k]);
+        for(k = 0; k < nx; k++)
+            work->x[k] += min_alpha * (work->x[k] - work->xold[k]);
 
     return add_ind;
 }
diff --git a/src/utils.c b/src/utils.c
index 64e78e7a..2fee54fc 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -101,6 +101,13 @@ int daqp_update_Rinv(DAQPWorkspace *work, c_float* H, int is_factored){
     int i, j, k, disp, disp2, disp3;
     const int n = work->n; 
     c_float eps = work->settings->eps_prox;
+    c_float zero_tol = work->settings->zero_tol;
+
+    // Reset the semi-proximal mask for this factorization
+    if(work->prox_mask != NULL){
+        for(i = 0; i < n; i++) work->prox_mask[i] = 0;
+    }
+    work->n_prox = 0;
 
     // Check if Diagonal
     int is_diagonal = 1;
@@ -122,12 +129,21 @@ int daqp_update_Rinv(DAQPWorkspace *work, c_float* H, int is_factored){
             c_float Hi = H[disp];
             // If factored, skip eps and sqrt, else apply them
             if(!is_factored){
-                Hi += eps;
-                if (Hi <= 0) return DAQP_EXIT_NONCONVEX;
+                // Only add eps if this direction would be singular without it
+                if(Hi <= zero_tol){
+                    if(work->prox_mask != NULL) work->prox_mask[i] = 1;
+                    work->n_prox++;
+                    Hi += eps;
+                }
+                if (Hi <= zero_tol) return DAQP_EXIT_NONCONVEX;
                 Hi = sqrt(Hi);
                 disp += n+1;
             } else {
-                if (eps != 0.0)  Hi = sqrt(Hi*Hi + eps); // Regularization for factors
+                if(Hi <= zero_tol){
+                    if(work->prox_mask != NULL) work->prox_mask[i] = 1;
+                    work->n_prox++;
+                    Hi = sqrt(Hi*Hi + eps); // Regularization for factors
+                }
                 disp += n-i;
             }
 
@@ -148,14 +164,23 @@ int daqp_update_Rinv(DAQPWorkspace *work, c_float* H, int is_factored){
                 work->Rinv[disp] = H[disp];
         }
     } else {
-        // Standard Cholesky (H -> R)
+        // Standard Cholesky (H -> R), adding eps only where the diagonal
+        // of the factor would otherwise be non-positive (semi-proximal).
         for (i=0, disp=0, disp3=0; i<n; disp+=n-i, i++, disp3+=i) {
-            work->Rinv[disp] = H[disp3++] + eps;
+            // Accumulate diagonal contribution without eps first
+            c_float diag_i = H[disp3++];
             for (k=0, disp2=i; k<i; k++, disp2+=n-k)
-                work->Rinv[disp] -= work->Rinv[disp2] * work->Rinv[disp2];
+                diag_i -= work->Rinv[disp2] * work->Rinv[disp2];
+
+            // Only regularize if this direction is singular
+            if(diag_i <= zero_tol){
+                if(work->prox_mask != NULL) work->prox_mask[i] = 1;
+                work->n_prox++;
+                diag_i += eps;
+            }
 
-            if (work->Rinv[disp] <= 0) return DAQP_EXIT_NONCONVEX;
-            work->Rinv[disp] = sqrt(work->Rinv[disp]);
+            if (diag_i <= zero_tol) return DAQP_EXIT_NONCONVEX;
+            work->Rinv[disp] = sqrt(diag_i);
 
             for (j=1; j<n-i; j++) {
                 work->Rinv[disp+j] = H[disp3++];