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/solver.jl b/src/solver.jl
index a10559f..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}.
+    - `regularize_linear_solve::Symbol = :identity`: 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 = LinearAlgebra.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) &&
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