From 160f02345797c0fb428ee92fcb6763f263539139 Mon Sep 17 00:00:00 2001 From: Dong Ho Lee Date: Thu, 25 Jun 2026 14:06:08 -0500 Subject: [PATCH 1/2] updated regularization schemes by adding Tikhonov and Marquardt --- src/MixedComplementarityProblems.jl | 2 +- src/solver.jl | 20 +++++++++++++++----- 2 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/MixedComplementarityProblems.jl b/src/MixedComplementarityProblems.jl index 98c98fa..12bfb42 100644 --- a/src/MixedComplementarityProblems.jl +++ b/src/MixedComplementarityProblems.jl @@ -1,7 +1,7 @@ module MixedComplementarityProblems using SparseArrays: SparseArrays -using LinearAlgebra: LinearAlgebra, I, norm, eigvals +using LinearAlgebra: LinearAlgebra, I, norm, eigvals, Diagonal using BlockArrays: blocks, blocksizes using TrajectoryGamesBase: to_blockvector using SymbolicTracingUtils: SymbolicTracingUtils as SymbolicTracingUtils diff --git a/src/solver.jl b/src/solver.jl index a10559f..a85e9fd 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -32,7 +32,7 @@ Keyword arguments: - `min_stepsize::Real = 1e-2`: the minimum step size for the linesearch. - `verbose::Bool = false`: whether to print debug information. - `linear_solve_algorithm::LinearSolve.SciMLLinearSolveAlgorithm`: the linear solve algorithm to use. Any solver from `LinearSolve.jl` can be used. - - `regularize_linear_solve::Symbol = :none`: scheme for regularizing the linear system matrix ∇F. Options are {:none, :identity, :internal}. + - `regularize_linear_solve::Symbol = :none`: scheme for regularizing the linear system matrix ∇F. Options are {:none, :identity, :internal, :Tikhonov, :Marquardt}. """ function solve( ::InteriorPoint, @@ -107,16 +107,26 @@ function solve( end if regularize_linear_solve === :identity - if size(∇F, 1) == size(∇F, 2) - linsolve.A = ∇F + η * I - else + if size(∇F, 1) != size(∇F, 2) @warn "Cannot use identity regularization on a nonsquare problem." + linsolve.A = ∇F + else + linsolve.A = ∇F + η * I end + linsolve.b = -F + elseif regularize_linear_solve === :Tikhonov + linsolve.A = (∇F' * ∇F) + η * I + linsolve.b = -∇F' * F + elseif regularize_linear_solve === :Marquardt + d = vec(sum(abs2, ∇F; dims = 1)) + D = Diagonal(d) + linsolve.A = (∇F' * ∇F) + η * D + linsolve.b = -∇F' * F else linsolve.A = ∇F + linsolve.b = -F end - linsolve.b = -F solution = solve!(linsolve) if !SciMLBase.successful_retcode(solution) && From 55ed94dd87022fc23d77be223a6bf4acc4ee535d Mon Sep 17 00:00:00 2001 From: Dong Ho Lee Date: Thu, 25 Jun 2026 16:50:29 -0500 Subject: [PATCH 2/2] update test to exercise different regularization options and ran the QP benchmarks --- benchmark/path.jl | 6 +++++- src/MixedComplementarityProblems.jl | 2 +- src/solver.jl | 4 ++-- test/runtests.jl | 30 +++++++++++++++++++++++++---- 4 files changed, 34 insertions(+), 8 deletions(-) diff --git a/benchmark/path.jl b/benchmark/path.jl index 12009dd..a67d65f 100644 --- a/benchmark/path.jl +++ b/benchmark/path.jl @@ -1,4 +1,8 @@ -"Benchmark interior point solver against PATH on a bunch of random test problems." +" +Benchmark interior point solver against PATH on a bunch of random test problems. +- benchmark_type: QuadraticProgramBenchmark(), TrajectoryGameBenchmark() +" + function benchmark( benchmark_type; num_samples = 100, diff --git a/src/MixedComplementarityProblems.jl b/src/MixedComplementarityProblems.jl index 12bfb42..98c98fa 100644 --- a/src/MixedComplementarityProblems.jl +++ b/src/MixedComplementarityProblems.jl @@ -1,7 +1,7 @@ module MixedComplementarityProblems using SparseArrays: SparseArrays -using LinearAlgebra: LinearAlgebra, I, norm, eigvals, Diagonal +using LinearAlgebra: LinearAlgebra, I, norm, eigvals using BlockArrays: blocks, blocksizes using TrajectoryGamesBase: to_blockvector using SymbolicTracingUtils: SymbolicTracingUtils as SymbolicTracingUtils diff --git a/src/solver.jl b/src/solver.jl index a85e9fd..d104773 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -32,7 +32,7 @@ Keyword arguments: - `min_stepsize::Real = 1e-2`: the minimum step size for the linesearch. - `verbose::Bool = false`: whether to print debug information. - `linear_solve_algorithm::LinearSolve.SciMLLinearSolveAlgorithm`: the linear solve algorithm to use. Any solver from `LinearSolve.jl` can be used. - - `regularize_linear_solve::Symbol = :none`: scheme for regularizing the linear system matrix ∇F. Options are {:none, :identity, :internal, :Tikhonov, :Marquardt}. + - `regularize_linear_solve::Symbol = :identity`: scheme for regularizing the linear system matrix ∇F. Options are {:none, :identity, :internal, :Tikhonov, :Marquardt}. """ function solve( ::InteriorPoint, @@ -119,7 +119,7 @@ function solve( linsolve.b = -∇F' * F elseif regularize_linear_solve === :Marquardt d = vec(sum(abs2, ∇F; dims = 1)) - D = Diagonal(d) + D = LinearAlgebra.Diagonal(d) linsolve.A = (∇F' * ∇F) + η * D linsolve.b = -∇F' * F else diff --git a/test/runtests.jl b/test/runtests.jl index 1b87ec9..2ce20cf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,6 +37,8 @@ using FiniteDiff: FiniteDiff @test sol.status == :solved end + regularize_linear_solve_options = (:none, :identity, :Tikhonov, :Marquardt) + @testset "BasicCallableConstructor" begin mcp = MixedComplementarityProblems.PrimalDualMCP( G, @@ -45,9 +47,19 @@ using FiniteDiff: FiniteDiff constrained_dimension = length(b), parameter_dimension = size(M, 1), ) - sol = MixedComplementarityProblems.solve(MixedComplementarityProblems.InteriorPoint(), mcp, θ) - check_solution(sol) + for regularize_linear_solve in regularize_linear_solve_options + @testset "$(regularize_linear_solve)" begin + sol = MixedComplementarityProblems.solve( + MixedComplementarityProblems.InteriorPoint(), + mcp, + θ; + regularize_linear_solve, + ) + + check_solution(sol) + end + end end @testset "AlternativeCallableConstructor" begin @@ -57,9 +69,19 @@ using FiniteDiff: FiniteDiff fill(Inf, size(M, 1) + length(b)); parameter_dimension = size(M, 1), ) - sol = MixedComplementarityProblems.solve(MixedComplementarityProblems.InteriorPoint(), mcp, θ) - check_solution(sol) + for regularize_linear_solve in regularize_linear_solve_options + @testset "$(regularize_linear_solve)" begin + sol = MixedComplementarityProblems.solve( + MixedComplementarityProblems.InteriorPoint(), + mcp, + θ; + regularize_linear_solve, + ) + + check_solution(sol) + end + end end @testset "AutodifferentationTests" begin