diff --git a/src/algorithms/pdhg.jl b/src/algorithms/pdhg.jl index 009210a..e800030 100644 --- a/src/algorithms/pdhg.jl +++ b/src/algorithms/pdhg.jl @@ -47,7 +47,7 @@ function initialize( η = fixed_stepsize(milp, algo.step_size) ω = one(η) step_sizes = StepSizes(; η, ω) - scratch = Scratch(; x = similar(sol.x), y = similar(sol.y), r = similar(sol.x)) + scratch = Scratch(sol) stats = ConvergenceStats(T; starting_time) state = PDHGState(; sol, sol_last, step_sizes, scratch, stats) return state diff --git a/src/algorithms/pdlp.jl b/src/algorithms/pdlp.jl index a56d7c6..6ae6fe5 100644 --- a/src/algorithms/pdlp.jl +++ b/src/algorithms/pdlp.jl @@ -55,7 +55,7 @@ function initialize( η = fixed_stepsize(milp, algo.step_size) ω = primal_weight_init(milp, algo.step_size) step_sizes = StepSizes(; η, ω) - scratch = Scratch(; x = similar(sol.x), y = similar(sol.y), r = similar(sol.x)) + scratch = Scratch(sol) iteration = IterationCounter(0, 0, 0) restart_stats = RestartStats(T) stats = ConvergenceStats(T; starting_time) diff --git a/src/components/errors.jl b/src/components/errors.jl index fc49be8..8d5632e 100644 --- a/src/components/errors.jl +++ b/src/components/errors.jl @@ -71,31 +71,37 @@ function kkt_errors!( (; c, lv, uv, A, At, lc, uc, D1, D2) = milp A_x = mul!(scratch.y, A, x) - At_y = mul!(scratch.x, At, y) - r = @. scratch.r = proj_multiplier(c - At_y, lv, uv) + c_At_y = c - mul!(scratch.x, At, y) + z = @. scratch.z = proj_multiplier(c_At_y, lv, uv) primal_diff = @. scratch.y = inv(D1.diag) * (A_x - proj_box(A_x, lc, uc)) primal = norm(primal_diff) + rescaled_combined_bounds = @. scratch.y = inv(D1.diag) * combine(lc, uc) primal_scale = one(T) + norm(rescaled_combined_bounds) - dual_diff = @. scratch.x = inv(D2.diag) * (c - At_y - r) + dual_diff = @. scratch.x = inv(D2.diag) * (c_At_y - z) dual = norm(dual_diff) + rescaled_obj = @. scratch.x = inv(D2.diag) * c dual_scale = one(T) + norm(rescaled_obj) + # dual objective: lᵀ|y|⁺ - uᵀ|y|⁻ + lᵥᵀ|z|⁺ - uᵥᵀ|z|⁻ + # We reformulate to ∑ⱼ (l⋅|y|⁺ - u⋅|y|⁻)ⱼ + ∑ᵢ (lᵥ⋅|z|⁺ - uᵥ⋅|z|⁻)ᵢ + # where pc = (l⋅|y|⁺ - u⋅|y|⁻) and pv = (lᵥ⋅|z|⁺ - uᵥ⋅|z|⁻) pc = @. scratch.y = ( - safeprod_left(uc, positive_part(-y)) - safeprod_left(lc, negative_part(-y)) + safeprod_left(lc, positive_part(y)) - safeprod_left(uc, negative_part(y)) ) - pv = @. scratch.r = ( - safeprod_left(uv, positive_part(-r)) - safeprod_left(lv, negative_part(-r)) + pv = @. scratch.z = ( + safeprod_left(lv, positive_part(z)) - safeprod_left(uv, negative_part(z)) ) pc_sum = sum(pc) pv_sum = sum(pv) cx = dot(c, x) + dobj = pc_sum + pv_sum - gap = abs(cx + pc_sum + pv_sum) - gap_scale = one(T) + abs(pc_sum + pv_sum) + abs(cx) + gap = abs(cx - dobj) + gap_scale = one(T) + abs(dobj) + abs(cx) err = KKTErrors(; primal, diff --git a/src/components/scratch.jl b/src/components/scratch.jl index 834e0a3..dafaa05 100644 --- a/src/components/scratch.jl +++ b/src/components/scratch.jl @@ -4,7 +4,7 @@ "dual scratch (length `ncons`)" y::V "dual scratch (length `nvar`)" - r::V + z::V end Scratch(sol::PrimalDualSolution) = Scratch(similar(sol.x), similar(sol.y), similar(sol.x)) diff --git a/test/components/errors.jl b/test/components/errors.jl index 5d416a4..17d82ca 100644 --- a/test/components/errors.jl +++ b/test/components/errors.jl @@ -7,7 +7,7 @@ function p(y, l, u) y⁻ = CoolPDLP.negative_part.(y) u_noinf = CoolPDLP.safe.(u) l_noinf = CoolPDLP.safe.(l) - return dot(y⁺, u_noinf) - dot(y⁻, l_noinf) + return dot(y⁺, l_noinf) - dot(y⁻, u_noinf) end milp, sol = CoolPDLP.random_milp_and_sol(100, 200, 0.4) @@ -27,10 +27,10 @@ err_p = CoolPDLP.kkt_errors!(scratch, sol_p, milp_p) @testset "Correct KKT errors" begin @test err.primal ≈ norm(A * x - CoolPDLP.proj_box.(A * x, lc, uc)) @test err.dual ≈ norm(c - At * y - r) - @test err.gap ≈ abs(dot(c, x) + p(-y, lc, uc) + p(-r, lv, uv)) + @test err.gap ≈ abs(dot(c, x) - (p(y, lc, uc) + p(r, lv, uv))) @test err.primal_scale ≈ 1 + norm(CoolPDLP.combine.(lc, uc)) @test err.dual_scale ≈ 1 + norm(c) - @test err.gap_scale ≈ 1 + abs(dot(c, x)) + abs(p(-y, lc, uc) + p(-r, lv, uv)) + @test err.gap_scale ≈ 1 + abs(dot(c, x)) + abs(p(y, lc, uc) + p(r, lv, uv)) end @testset "Invariance by preconditioning" begin