ACE 592 - Lecture 3.3

Complementarity Problems

Author
Affiliation

Diego S. Cardoso

University of Illinois Urbana-Champaign

1 Course Roadmap

  1. Introduction to Scientific Computing
  2. Fundamentals of numerical methods
  3. Systems of equations
    1. Linear systems
    2. Nonlinear systems
    3. Complementarity problems
  4. Optimization
  5. Function approximation
  6. Structural estimation

2 Agenda

  • Today, we will cover complementarity problems, which are a generalization of rootfinding problems that allow for constraints on the solution
  • Here, we cover how to set up and solve mixed complementarity problems in Julia

3 Main references

  • Miranda & Fackler (2002), Chs. 3.7 and 3.8
  • Lecture notes for Ivan Rudik’s Dynamic Optimization (Cornell)

4 Complementarity problems

In these problems, we have

  • A function f: \mathbb{R}^n \rightarrow \mathbb{R}^n
  • n-vectors a and b, with a<b

And we are looking for an n-vector x \in [a,b] such that for all i = 1,\dots,n

\begin{align*} x_i > a_i \Rightarrow f_i(x) \geq 0 \\ x_i < b_i \Rightarrow f_i(x) \leq 0 \end{align*}

. . .

  • This formulation is usually referred to as a Mixed Complementarity Problem because it has an upper and lower bounds
    • Standard complementarity problems have one-sided bounds

5 Complementarity problems

An economic interpretation

  • x is an n-dimensional vector of some economic action
  • a_i and b_i are the lower and upper bounds on action i
  • f_i(x) is the marginal arbitrage profit of action i

. . .

There are disequilibrium profit opportunities if

  1. x_i < b_i and f_i(x) > 0 (we can increase profits by raising x_i)

. . .

  1. x_i > a_i and f_i(x) < 0 (we can increase profits by decreasing x_i)

6 Complementarity problems

We obtain a no-arbitrage equilibrium if and only if x solves the complementary problem CP(f,a,b)

. . .

We can write out the problem as finding a vector x \in [a,b] that solves

\begin{align} x_i > a_i \Rightarrow f_i(x) \geq 0 \,\,\, \forall i = 1,...,n \notag\\ x_i < b_i \Rightarrow f_i(x) \leq 0 \,\,\, \forall i = 1,...,n \notag \end{align}

. . .

At interior solution, the function must be precisely be zero

Corner solution at the upper bound b_i for x_i \rightarrow f must be increasing in direction i

The opposite is true if we are at the lower bound

7 Complementarity problems

Most economic problems are complementarity problems where we are

  • Looking for a root of a function (e.g. marginal profit)
  • Subject to some constraint (e.g. price floors)

. . .

The Karush-Kuhn-Tucker theorem shows that x solves the constrained optimization problem ( \max_x F(x) subject to x \in [a,b]) only if it solves the complementarity problem CP(f, a, b), where f_i = \partial F/\partial x_i

. . .

Let’s use a linear f to visualize what do we mean by a solution in complementarity problems

8 Complementarity problems

Case 1

What is the solution?

x^{*} = a, with f(x^{*}) < 0

Another way of seeing it: imagine we’re trying to maximize F

What would F look like between a and b?

  • The decreasing part of a concave parabola

9 Complementarity problems

Case 2

What is the solution?

x^{*} = b, with f(x^{*}) > 0

Once again, imagine we’re trying to maximize F

What would F look like between a and b?

  • The increasing part of a concave parabola

10 Complementarity problems

Case 3

What is the solution?

Some x^{*} between a and b, with f(x^{*}) = 0

What would F look like between a and b?

  • A concave parabola with an interior maximum

11 Complementarity problems

Case 4

What is the solution?

Actually, we have 3 solutions:

  1. x^{*} = a, with f(x^{*}) < 0
  2. x^{*} = b, with f(x^{*}) > 0
  3. Some x^{*} \in (a,b), with f(x^{*}) = 0
  • And this is an unstable solution

What would F look like between a and b?

  • A convex parabola! So we end up with multiple local maxima that satisfy the 1st-order condition

12 Solving CPs

A complementarity problem CP(f, a, b) can be re-framed as a rootfinding problem

\hat{f}(x) = min(max(f(x),a-x),b-x) = 0

. . .

Let’s revisit those 4 cases to understand why this works

13 Solving CPs

Case 1

\hat{f}(x) = min(max(f(x),a-x),b-x) = 0

  • Dashed lines are b-x and a-x
  • Thin solid line is f(x)
  • Thick solid line is \hat{f}(x)

\rightarrow The solution is x^{*} = a

  • Note that f(x) < 0 for all x \in [a,b], so this is still case 1

14 Solving CPs

Case 2:

\hat{f}(x) = min(max(f(x),a-x),b-x) = 0

In this case, f(x) > 0 for all x \in [a,b]

\rightarrow The solution is x^{*} = b

15 Solving CPs

Case 3:

\hat{f}(x) = min(max(f(x),a-x),b-x) = 0

In this case, f(a) > 0, f(b) < 0

\rightarrow The solution is some x^{*} \in (a, b)

16 Solving CPs

Case 4:

\hat{f}(x) = min(max(f(x),a-x),b-x) = 0

In this case, f(a) < 0, f(b) > 0

\rightarrow We have the same 3 solutions

  1. x^{*} = a, with f(x^{*}) < 0
  2. x^{*} = b, with f(x^{*}) > 0
  3. Some x^{*} \in (a,b), with f(x^{*}) = 0

17 Solving CPs

Once \hat{f} is defined, we can use Newton’s or quasi-Newton methods to solve a CP

If using Newton, we need to define the Jacobian \hat{J}(x) with row i being

  • \hat{J}_i(x) = J_i(x), if a_i - x_i < f_i(x) < b_i - x_i
  • \hat{J}_i(x) = -I_i(x), otherwise

where I is the identity matrix

18 Solving CPs

Rootfinding with \hat{f} works well in many cases but can be problematic in others

One problem is that \hat{f} has nondiferentiable kinks. This can lead to slower convergence, cycles, and incorrect answers1

A workaround is to use an alternative function with smoother transitions, such as Fischer’s function

\tilde{f}(x) = \phi^{-}(\phi^{+}(f(x), a-x), b-x)

where \phi_i^\pm(u, v) = u_i + v_i \pm \sqrt{u_i^2 + v_i^2}

19 Solving CPs

20 Solving CPs in Julia

We can use NLsolve.mcpsolve to solve CPs for us.2 Let’s see an example of f:\mathbb{R}^2\rightarrow\mathbb{R}^2

\begin{align*} f_1(x_1, x_2) & = 3x_1^2 + 2x_1x_2 + 2x_2^2 + 4x_1 - 2 \\ f_2(x_1, x_2) & = 2x_1^2 + 5x_1x_2 - x_2^2 + x_2 + 1 \end{align*}

Code 
function f!(F, x) # We need to program the memory-efficient version here
    F[1] = 3*x[1]^2 + 2*x[1]*x[2] + 2*x[2]^2 + 4*x[1] - 2
    F[2] = 2*x[1]^2 + 5*x[1]*x[2] -   x[2]^2 + x[2] + 1
end;

21 Solving CPs in Julia

Let’s check the solution to the standard rootfinding problem using Newton’s method

Code 
using NLsolve;
NLsolve.nlsolve(f!, [1.0; 1.0], method=:newton)
Results of Nonlinear Solver Algorithm
 * Algorithm: Newton with line-search
 * Starting Point: [1.0, 1.0]
 * Zero: [-0.13736984903421828, 1.1872338075744588]
 * Inf-norm of residuals: 0.000000
 * Iterations: 7
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 8
 * Jacobian Calls (df/dx): 8

22 Solving CPs in Julia

Getting back to CPs, let’s impose non-negativity bounds on x_1 and x_2 but no upper bounds (i.e., b = \infty)

Code 
a = [0.0; 0.0]; b = [Inf; Inf];
r = NLsolve.mcpsolve(f!, a, b, [1.0; 1.0], method=:newton, reformulation=:minmax)
Results of Nonlinear Solver Algorithm
 * Algorithm: Newton with line-search
 * Starting Point: [1.0, 1.0]
 * Zero: [0.3874258867227982, 0.0]
 * Inf-norm of residuals: 0.000000
 * Iterations: 10
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 17
 * Jacobian Calls (df/dx): 10

23 Solving CPs in Julia

We can get the value of the root using r.zero. Checking f(x^*)

Code 
r.zero
2-element Vector{Float64}:
 0.3874258867227982
 0.0
Code 
return_f = zeros(2);
f!(return_f, r.zero);
return_f
2-element Vector{Float64}:
 3.241851231905457e-14
 1.300197635405893

So the non-negativity constraint is binding for x_2 but not x_1

24 Solving CPs in Julia

mcpsolve takes the same additional arguments nlsolve had: xtol, ftol, iterations, autodiff, show_trace, etc

In addition, it accepts argument reformulation, which can be smooth (default) or minmax

Code 
NLsolve.mcpsolve(f!, a, b, [1.0; 1.0], method=:newton, reformulation=:minmax, autodiff=:forward)
Results of Nonlinear Solver Algorithm
 * Algorithm: Newton with line-search
 * Starting Point: [1.0, 1.0]
 * Zero: [0.38742588672279815, 0.0]
 * Inf-norm of residuals: 0.000000
 * Iterations: 10
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 17
 * Jacobian Calls (df/dx): 10

25 Solving CPs in Julia

IMPORTANT NOTE. Miranda and Fackler’s textbook (and Matlab package) flip the sign convention for MCP problems. That’s because they are formulating it with an economic context in mind, where these problems arise from constrained optimization.

. . .

The conventional setup is (note the flipped inequalities)

\begin{align*} x_i > a_i \Rightarrow f_i(x) \leq 0 \\ x_i < b_i \Rightarrow f_i(x) \geq 0 \end{align*}

. . .

Solution. If you are using standard MCP solvers (like mcpsolve), flip the sign of your f function

  • For more details, see Miranda and Fackler’s Bibliographic Notes after Chapter 3

26 Example: Cournot duopoly with capacity constraints

Let’s revisit our duopoly example from Lecture 3.2 and add capacity constraints

. . .

Same setup as before:

  • Inverse demand: P(Q) = \theta Q^{-1/\epsilon}, where we add \theta as demand shifter
  • Costs: C_i(q_i) = \frac{1}{2}c_i q_i^2

. . .

New: each firm has a capacity constraint q_i \in [0, K_i]

  • Firm i maximizes \pi_i = P(Q)q_i - C_i(q_i) subject to 0 \leq q_i \leq K_i

27 Example: Cournot duopoly as a CP

The KKT conditions for firm i’s problem are exactly the complementarity conditions

\begin{align*} q_i > 0 &\Rightarrow \frac{\partial \pi_i}{\partial q_i} \leq 0 \\ q_i < K_i &\Rightarrow \frac{\partial \pi_i}{\partial q_i} \geq 0 \end{align*}

. . .

where marginal profit is

f_i(q_1, q_2) = \theta Q^{-1/\epsilon} - \frac{\theta}{\epsilon} Q^{-1/\epsilon-1}q_i - c_i q_i

. . .

This is CP(f, a, b) with a = (0, 0) and b = (K_1, K_2)

28 Sign convention for mcpsolve

Recall that NLsolve.mcpsolve uses the convention

\begin{align*} x_i > a_i &\Rightarrow f_i(x) \geq 0 \\ x_i < b_i &\Rightarrow f_i(x) \leq 0 \end{align*}

But our marginal profit f_i > 0 means the firm wants to increase q_i. At the upper bound q_i = K_i, we expect f_i \geq 0 (wants to produce more but can’t)

. . .

The KKT convention has the opposite sign: q_i < b_i \Rightarrow f_i \geq 0, but mcpsolve expects f_i \leq 0

We need to pass -f to mcpsolve

29 Example: best responses with capacity constraints

Without constraints, the Nash equilibrium is (q_1^*, q_2^*) \approx (0.84, 0.69)

But firm 1 (the low-cost firm) can’t produce more than K_1 = 0.7!

. . .

The constrained equilibrium: q_1^* = 0.7 (at capacity) and q_2^* \approx 0.70 (interior)

  • Firm 1’s marginal profit is positive at q_1 = K_1: it would produce more if it could
  • Firm 2’s marginal profit is zero: it is at its unconstrained optimum

. . .

This is precisely the complementarity structure: firm 1 is at the upper bound with f_1 > 0, while firm 2 is interior with f_2 = 0

30 Example: defining the CP with NLsolve

We define the functions kind of like before, but now using a vector of parameters and flipping the sign convention

Code 
epsilon = 1.6; c = [0.6; 0.8]; theta = 1.0; K = [0.7; 1.0];

function neg_foc!(F, q)
    Q = sum(q)
    for i in 1:2
        # NOTE: we flip the sign of f_i to match mcpsolve's convention
        F[i] = -(theta * Q^(-1/epsilon) - (theta/epsilon) * Q^(-1/epsilon - 1) * q[i] - c[i] * q[i])
    end
end
neg_foc! (generic function with 1 method)

31 Example: solving the CP with NLsolve

Code 
using NLsolve
a = [0.0; 0.0]; b = K;
r = NLsolve.mcpsolve(neg_foc!, a, b, [0.2, 0.2], method=:newton, autodiff=:forward)
Results of Nonlinear Solver Algorithm
 * Algorithm: Newton with line-search
 * Starting Point: [0.2, 0.2]
 * Zero: [0.7000000005741004, 0.6976530167583798]
 * Inf-norm of residuals: 0.000000
 * Iterations: 5
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 6
 * Jacobian Calls (df/dx): 5

32 Interpreting the solution

Remember to re-flip the sign of f to match the KKT convention

Code 
F_check = zeros(2); # Create a vector to store the value of f at the solution
neg_foc!(F_check, r.zero); # Compute f at the solution
println("f(q*) = $(-F_check)")
f(q*) = [0.13727104736212875, -1.489608436600065e-11]

. . .

  • q_1^* = K_1 = 0.7 with f_1 > 0: firm 1 is capacity-constrained. Marginal profit is positive: it wants to produce more but can’t
  • q_2^* \approx 0.70 < K_2 with f_2 \approx 0: firm 2 is at its interior optimum

. . .

Compare with the unconstrained solution from Lecture 3.2: [0.8396, 0.6888]

Without the constraint, firm 1 would produce \approx 0.84 > K_1

33 Best responses with capacity constraints

With capacity constraints, best response functions become truncated

34 Demand shifts and capacity constraints

Now let’s ask: what happens to the equilibrium as demand increases?

We vary \theta from 0.3 to 3.0

. . .

As \theta rises, both firms want to produce more. Three regimes emerge:

  1. Both interior (\theta < 0.75): neither constraint binds, both FOCs hold with equality
  2. Firm 1 at capacity (0.75 \leq \theta < 1.77): firm 1 maxes out first (lower cost \rightarrow higher desired output), firm 2 expands
  3. Both at capacity (\theta \geq 1.77): demand is so high that both firms are constrained

35 Equilibrium locus

36 Tracing the equilibrium locus

Code 
using Plots

thetas = range(0.3, 3.0, length=100)
q1_path = Float64[]; q2_path = Float64[];
q_init = [0.3; 0.25]

for θ in thetas
    function neg_foc_theta!(F, q)
        Q = sum(q)
        for i in 1:2
            F[i] = -(θ * Q^(-1/epsilon) - (θ/epsilon) * Q^(-1/epsilon - 1) * q[i] - c[i] * q[i])
        end
    end
    r = NLsolve.mcpsolve(neg_foc_theta!, a, b, q_init, method=:newton, autodiff=:forward)
    push!(q1_path, r.zero[1])
    push!(q2_path, r.zero[2])
    q_init = r.zero  # warm start for next iteration
end

37 Tracing the equilibrium locus

Code 
plot(collect(thetas), q1_path, label="q₁*", lw=2,
     xlabel="θ (demand shifter)", ylabel="Equilibrium output")
plot!(collect(thetas), q2_path, label="q₂*", lw=2)
hline!([0.7], label="K₁", ls=:dash, color=:blue, alpha=0.5)
hline!([1.0], label="K₂", ls=:dash, color=:red, alpha=0.5)

. . .

The kinks in the equilibrium paths correspond to transitions between regimes where different capacity constraints become binding

38 Takeaways from this example

  1. CP conditions arise naturally from constrained optimization: the KKT conditions of the firms’ problems are exactly the CP conditions

. . .

  1. Sign conventions matter: mcpsolve uses a different sign convention from the standard KKT formulation. Always check the documentation and verify the solution makes economic sense

. . .

  1. The CP framework handles regime changes automatically: as \theta varies, the solver transitions between interior and corner solutions without us having to enumerate cases manually

Footnotes

  1. See Miranda & Fackler (2002) Ch. 3.8 for a pathological example that needs smoothing.↩︎

  2. Because complementarity conditions are constraints, a more comprehensive way for larger systems is through JuMP.jl modeling interface. We will cover that in the Constrained Optimization unit.↩︎