diff --git a/benchmark/Project.toml b/benchmark/Project.toml
new file mode 100644
index 0000000..b781345
--- /dev/null
+++ b/benchmark/Project.toml
@@ -0,0 +1,11 @@
+[deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+MixedComplementarityProblems = "6c9e26cb-9263-41b8-a6c6-f4ca104ccdcd"
+PATHSolver = "f5f7c340-0bb3-5c69-969a-41884d311d1b"
+ParametricMCPs = "9b992ff8-05bb-4ea1-b9d2-5ef72d82f7ad"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[sources]
+MixedComplementarityProblems = {path = ".."}
diff --git a/benchmark/SolverBenchmarks.jl b/benchmark/SolverBenchmarks.jl
new file mode 100644
index 0000000..a93b40a
--- /dev/null
+++ b/benchmark/SolverBenchmarks.jl
@@ -0,0 +1,14 @@
+"Module for benchmarking different solvers against one another."
+module SolverBenchmarks
+
+using MixedComplementarityProblems: MixedComplementarityProblems
+using ParametricMCPs: ParametricMCPs
+using Random: Random
+using Statistics: Statistics
+using Distributions: Distributions
+using PATHSolver: PATHSolver
+using ProgressMeter: @showprogress
+
+include("path.jl")
+
+end # module SolverBenchmarks
diff --git a/benchmark/path.jl b/benchmark/path.jl
new file mode 100644
index 0000000..9f5e075
--- /dev/null
+++ b/benchmark/path.jl
@@ -0,0 +1,168 @@
+""" Generate a random large (convex) quadratic problem of the form
+                               min_x 0.5 xᵀ M x - ϕᵀ x
+                               s.t.  Ax - b ≥ 0.
+
+NOTE: the problem may not be feasible!
+"""
+function generate_test_problem(; num_primals, num_inequalities)
+    G(x, y; θ) =
+        let
+            (; M, A, ϕ) = unpack_parameters(θ; num_primals, num_inequalities)
+            M * x - ϕ - A' * y
+        end
+
+    H(x, y; θ) =
+        let
+            (; A, b) = unpack_parameters(θ; num_primals, num_inequalities)
+            A * x - b
+        end
+
+    K(z, θ) =
+        let
+            x = z[1:num_primals]
+            y = z[(num_primals + 1):end]
+
+            [G(x, y; θ); H(x, y; θ)]
+        end
+
+    (; G, H, K)
+end
+
+"Generate a random parameter vector Θ corresponding to a convex QP."
+function generate_random_parameter(rng; num_primals, num_inequalities, sparsity_rate)
+    bernoulli = Distributions.Bernoulli(1 - sparsity_rate)
+
+    M = let
+        P =
+            randn(rng, num_primals, num_primals) .*
+            rand(rng, bernoulli, num_primals, num_primals)
+        P' * P
+    end
+
+    A =
+        randn(rng, num_inequalities, num_primals) .*
+        rand(rng, bernoulli, num_inequalities, num_primals)
+    b = randn(rng, num_inequalities)
+    ϕ = randn(rng, num_primals)
+
+    [reshape(M, length(M)); reshape(A, length(A)); b; ϕ]
+end
+
+"Unpack a parameter vector θ into the components of a convex QP."
+function unpack_parameters(θ; num_primals, num_inequalities)
+    M = reshape(θ[1:(num_primals^2)], num_primals, num_primals)
+    A = reshape(
+        θ[(num_primals^2 + 1):(num_primals^2 + num_inequalities * num_primals)],
+        num_inequalities,
+        num_primals,
+    )
+
+    b =
+        θ[(num_primals^2 + num_inequalities * num_primals + 1):(num_primals^2 + num_inequalities * (num_primals + 1))]
+    ϕ = θ[(num_primals^2 + num_inequalities * (num_primals + 1) + 1):end]
+
+    (; M, A, b, ϕ)
+end
+
+"Benchmark interior point solver against PATH on a bunch of random QPs."
+function benchmark(;
+    num_samples = 1000,
+    num_primals = 100,
+    num_inequalities = 100,
+    sparsity_rate = 0.9,
+    ip_mcp = nothing,
+    path_mcp = nothing,
+    ip_kwargs = (;),
+)
+    rng = Random.MersenneTwister(1)
+
+    # Generate problem and random parameters.
+    @info "Generating random problems..."
+    problem = generate_test_problem(; num_primals, num_inequalities)
+
+    θs = map(1:num_samples) do _
+        generate_random_parameter(rng; num_primals, num_inequalities, sparsity_rate)
+    end
+
+    # Generate corresponding MCPs.
+    @info "Generating IP MCP..."
+    parameter_dimension = length(first(θs))
+    ip_mcp =
+        !isnothing(ip_mcp) ? ip_mcp :
+        MixedComplementarityProblems.PrimalDualMCP(
+            problem.G,
+            problem.H;
+            unconstrained_dimension = num_primals,
+            constrained_dimension = num_inequalities,
+            parameter_dimension,
+        )
+
+    @info "Generating PATH MCP..."
+    lower_bounds = [fill(-Inf, num_primals); fill(0, num_inequalities)]
+    upper_bounds = fill(Inf, num_primals + num_inequalities)
+    path_mcp =
+        !isnothing(path_mcp) ? path_mcp :
+        ParametricMCPs.ParametricMCP(
+            problem.K,
+            lower_bounds,
+            upper_bounds,
+            parameter_dimension,
+        )
+
+    # Warm up the solvers.
+    @info "Warming up IP solver..."
+    MixedComplementarityProblems.solve(
+        MixedComplementarityProblems.InteriorPoint(),
+        ip_mcp,
+        first(θs);
+        ip_kwargs...,
+    )
+
+    @info "Warming up PATH solver..."
+    ParametricMCPs.solve(path_mcp, first(θs); warn_on_convergence_failure = false)
+
+    # Solve and time.
+    ip_data = @showprogress desc = "Solving IP MCPs..." map(θs) do θ
+        elapsed_time = @elapsed sol = MixedComplementarityProblems.solve(
+            MixedComplementarityProblems.InteriorPoint(),
+            ip_mcp,
+            θ,
+        )
+
+        (; elapsed_time, success = sol.status == :solved)
+    end
+
+    path_data = @showprogress desc = "Solving PATH MCPs..." map(θs) do θ
+        # Solve and time.
+        elapsed_time = @elapsed sol =
+            ParametricMCPs.solve(path_mcp, θ; warn_on_convergence_failure = false)
+
+        (; elapsed_time, success = sol.status == PATHSolver.MCP_Solved)
+    end
+
+    (; ip_mcp, path_mcp, ip_data, path_data)
+end
+
+"Compute summary statistics from solver benchmark data."
+function summary_statistics(data)
+    accumulate_stats(solver_data) = begin
+        (; success_rate = fraction_solved(solver_data), runtime_stats(solver_data)...)
+    end
+
+    (; ip = accumulate_stats(data.ip_data), path = accumulate_stats(data.path_data))
+end
+
+"Estimate mean and standard deviation of runtimes for all problems."
+function runtime_stats(solver_data)
+    filtered_times =
+        map(datum -> datum.elapsed_time, filter(datum -> datum.success, solver_data))
+    μ = Statistics.mean(filtered_times)
+    σ = Statistics.stdm(filtered_times, μ)
+
+    (; μ, σ)
+end
+
+"Compute fraction of problems solved."
+function fraction_solved(solver_data)
+    Statistics.mean(datum -> datum.success, solver_data)
+end
diff --git a/src/solver.jl b/src/solver.jl
index 1ba0410..e43db3f 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -25,6 +25,10 @@ function solve(
     tol = 1e-4,
     max_inner_iters = 20,
     max_outer_iters = 50,
+    tightening_rate = 0.1,
+    loosening_rate = 0.5,
+    min_stepsize = 1e-2,
+    verbose = false,
 )
     # Set up common memory.
     ∇F = mcp.∇F_z!.result_buffer
@@ -52,14 +56,14 @@ function solve(
             # TODO! Can add some adaptive regularization.
             mcp.F!(F, x, y, s; θ, ϵ)
             mcp.∇F_z!(∇F, x, y, s; θ, ϵ)
-            δz .= (∇F \ F) .* -1
+            δz .= ((∇F + tol * I) \ F) .* -1
 
             # Fraction to the boundary linesearch.
-            α_s = fraction_to_the_boundary_linesearch(s, δs; tol)
-            α_y = fraction_to_the_boundary_linesearch(y, δy; tol)
+            α_s = fraction_to_the_boundary_linesearch(s, δs; tol = min_stepsize)
+            α_y = fraction_to_the_boundary_linesearch(y, δy; tol = min_stepsize)
 
             if isnan(α_s) || isnan(α_y)
-                @warn "Linesearch failed. Exiting prematurely."
+                verbose && @warn "Linesearch failed. Exiting prematurely."
                 status = :failed
                 break
             end
@@ -73,7 +77,9 @@ function solve(
             inner_iters += 1
         end
 
-        ϵ *= (status == :solved) ? 1 - exp(-inner_iters) : 1 + exp(-inner_iters)
+        ϵ *=
+            (status == :solved) ? 1 - exp(-tightening_rate * inner_iters) :
+            1 + exp(-loosening_rate * inner_iters)
         outer_iters += 1
     end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 716815f..1c2fd84 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -18,7 +18,7 @@ using FiniteDiff: FiniteDiff
     b = [1; 1]
     θ = rand(2)
 
-    G(x, y; θ) = M * x - A' * y - θ
+    G(x, y; θ) = M * x - θ - A' * y
     H(x, y; θ) = A * x - b
     K(z; θ) = begin
         x = z[1:size(M, 1)]