Nonlinear complementarity problem methods

Nonlinear complementarity problem methods¶

Randall Romero Aguilar, PhD

This demo is based on the original Matlab demo accompanying the Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.

Original (Matlab) CompEcon file: demslv08.m

Running this file requires the Python version of CompEcon. This can be installed with pip by running

!pip install compecon --upgrade

Last updated: 2022-Sept-04

About¶

Solve nonlinear complementarity problem on \(R^2\) using semismooth and minmax methods.

Function to be solved is

\[\begin{equation*} f(x,y) = \begin{bmatrix} 200x(y - x^2) + 1 - x\\100(x^2 - y)\end{bmatrix} \end{equation*}\]

Set up the problem¶

The class MCP is used to represent mixed-complementarity problems. To create one instance, we define the objective function and the boundaries \(a\) and \(b\) such that for \(a \leq x \leq b\). Apart from the required parameters, we can specify options to be used when solving the problem.

import numpy as np
import pandas as pd
from compecon import MCP, tic, toc

def f(z):
    x, y = z
    fval = [200*x*(y - x**2) + 1 - x,
            100*(x**2 - y)]
    fjac = [[200*(y - x**2) - 400*x**2 - 1, 200*x],
            [200*x, -100]]

    return fval, fjac

Generate problem test data¶

z = 2 * np.random.randn(2, 2)
a = 1 + np.min(z, 0)
b = 1 + np.max(z, 0)

F = MCP(f, a, b, maxit=1500)

Solve by applying Newton method¶

We’ll use a random initial guess \(x\)

F.x0 = np.random.randn(2)

Check the Jacobian¶

F.user_provides_jacobian = True
F.check_jacobian()

All numerical derivatives differ from
the user-provided ones by less than 8 decimal digits.

The maximum error is 7.30e-09, for row 0 and column 0.

Using minmax formulation¶

t0 = tic()
x1 = F.newton(transform='minmax')
t1 = 100 * toc(t0)
n1 = F.fnorm

Using semismooth formulation¶

t0 = tic()
x2 = F.newton(transform='ssmooth')
t2 = 100*toc(t0)
n2 = F.fnorm

Print results¶

print(
    'Hundreds of seconds required to solve nonlinear complementarity \n' +
    'problem on R^2 using minmax and semismooth formulations, with \n'  +
    'randomly generated bounds \n' +
    f'\ta = [{a[0]:4.2f}, {a[1]:4.2f}] \n' + 
    f'\tb = [{b[0]:4.2f}, {b[1]:4.2f}]'
    )

results = pd.DataFrame.from_records(
    [('Newton minmax',t1, n1, *x1), 
     ('Newton semismooth', t2, n2, *x2)],
    columns = ['Algorithm', 'Time', 'Norm', 'x1', 'x2']
    ).set_index('Algorithm')

results.style.format(precision=3, subset=['Time', 'x1', 'x2'])

Hundreds of seconds required to solve nonlinear complementarity 
problem on R^2 using minmax and semismooth formulations, with 
randomly generated bounds 
	a = [2.08, -0.08] 
	b = [4.60, 4.42]

	Time	Norm	x1	x2
Algorithm
Newton minmax	1.100	0.000000	2.083	4.337
Newton semismooth	0.500	0.000000	2.083	4.337

Linear complementarity problem methods

Hard nonlinear complementarity problem with Billup’s function

A Python Implementation of CompEcon