Solving a linear complementarity problem#
This tutorial walks through solving a linear complementarity problem (LCP) with Lemke’s algorithm, from creating the problem to interpreting the result and watching the algorithm run.
Given a vector \(q\) and a matrix \(M\), an LCP asks to find a vector \(z \ge 0\) such that
To get started, Lemke’s algorithm also needs a covering vector \(d\), which multiplies the extra variable \(z_0\). So, a problem is described by the triple \((M, q, d)\). All arithmetic is exact, using rational numbers.
Creating an LCP#
An LCP is represented by the lcp class.
It can be created in two ways.
Let’s consider an LCP, where
1) From Python#
Pass the matrix M (as a list of its rows) and the vectors q and d:
>>> from fractions import Fraction as F
>>> from lemke.lemke import lcp
>>> problem = lcp(
... M=[
... [F(1, 2), F(1)],
... [F(-1), F(5, 2)],
... ],
... q=[F(-1), F(-1)],
... d=[F(1), F(1)],
... )
>>> print(problem.n)
2
Note
Entries must be exact rationals, such as Python integers or
fractions.Fraction instances. Floats are not supported.
2) From a file#
Alternatively, describe the problem in a text file and load it with
from_file().
First, create a file such as examples/lcp_2x2:
n= 2
M= 1/2 1
-1 2.5
q= -1 -1
d= 1 1
then load it:
>>> problem = lcp.from_file("examples/lcp_2x2")
>>> print(problem.n)
2
See details about how the input file should be formatted on the file formats page.
Solving the problem#
Pass the problem to runlemke():
>>> from lemke.lemke import runlemke
>>> result = runlemke(lcp=problem)
>>> print(result.success)
True
runlemke() always returns an LcpResult,
whether or not a solution was found. Its fields are the same in both cases,
so you must check success field first.
Note
If every entry of q is non-negative, the trivial solution
\(z = 0,\ w = q\) is returned immediately without any pivoting.
Reading the result#
The solution values are available as tuples of fractions.Fraction:
>>> print(result.z)
(Fraction(2, 3), Fraction(2, 3))
>>> print(result.w)
(Fraction(0, 1), Fraction(0, 1))
Printing the result shows the basis and the values in aligned columns, so that the complementarity condition (\(z_i w_i = 0\)) can be checked by eye:
>>> print(result)
Process finished successfully after 3 pivots.
Solution found:
basis= z1 z2
z= 0 2/3 2/3
w= 0 0
The first column of z= is the artificial variable \(z_0\), which is
always 0 for a successful run. See LcpResult for all
fields, such as basis and num_pivots.
Degenerate solutions#
The solution above is non-degenerate: for every \(i\), exactly one of \(z_i\) and \(w_i\) is positive.
A solution is degenerate if some basic variable has value zero, meaning both \(z_i\) and \(w_i\) are zero for some \(i\). It is still a valid solution. For example:
>>> problem = lcp(
... M=[
... [F(1), F(0)],
... [F(0), F(1)],
... ],
... q=[F(-1), F(0)],
... d=[F(1), F(1)],
... )
>>> result = runlemke(lcp=problem)
>>> print(result.z)
(Fraction(1, 1), Fraction(0, 1))
>>> print(result.w)
(Fraction(0, 1), Fraction(0, 1))
Here \(z_2 = w_2 = 0\), so the second complementary pair is degenerate.
Degeneracy is handled by the lexicographic minimum ratio test, which
prevents the algorithm from cycling. You can inspect how often that test
was needed with the lexstats option of the printing callback (see below).
Since in this case it is not clear which variables are basic,
you can explicitly check the basis field:
>>> print(sorted(result.basis))
['w2', 'z1']
Ray termination#
Lemke’s algorithm can also end on a secondary ray, when the variable that should enter the basis has no positive pivot candidate. In that case no valid solution is found.
Consider the following example, where \(w = -1 + 0 * z\) can never be non-negative:
>>> problem = lcp(
... M=[[F(0)]],
... q=[F(-1)],
... d=[F(1)],
... )
>>> result = runlemke(lcp=problem)
>>> print(result.success)
False
>>> print(result.ray_entering_variable)
z1
When success is False:
ray_entering_variablenames the entering variable that led to the secondary ray,z0,zandwdescribe the final values reached, which is not a solution of the LCP.
>>> print(result)
Terminated on a secondary ray after 1 pivot, when trying to enter z1.
Current basis not an LCP solution:
basis= z0
z= 1 0
w= 0
Note
Ray termination does not, by itself, prove that the LCP has no solution. It only means the algorithm could not find one from this covering vector.
Watching the algorithm#
To see what happens at each pivot, pass a callback to runlemke().
The ready-made PrintingCallback prints the LCP, the
tableaus and the chosen leaving and entering variables:
>>> from lemke.lemke import PrintingCallback
>>> problem = lcp.from_file("examples/lcp_2x2")
>>> result = runlemke(lcp=problem, callback=PrintingCallback())
verbose=False z0=False lexstats=False
n= 2
M=
1/2 1
-1 5/2
q=
-1 -1
d=
1 1
==================================
After filltableau:
Determinant: -1
var z0 z1 z2 RHS
scfa 1 2 2 1
w1 1 1 2 -1
w2 1 -2 5 -1
-----------------end of tableau-----------------
leaving: w2 entering: z0
leaving: w1 entering: z2
leaving: z0 entering: z1
Final tableau:
Determinant: 9
var w2 z0 w1 RHS
scfa 1 1 1 1
z2 -1 3 -2 3
z1 2 3 -5 3
-----------------end of tableau-----------------
Process finished successfully after 3 pivots.
Solution found:
basis= z1 z2
z= 0 2/3 2/3
w= 0 0
The options are:
verbose: print the full tableau after every pivot,z0: print the current value of \(z_0\) before each pivot,lexstats: print minimum-ratio test statistics at the end,stream: write to a file-like object instead ofsys.stdout.
Writing your own callback#
To react to the algorithm in your own way, subclass
LemkeCallback and override only the methods you need.
All of them do nothing by default.
For example, to log entering variables during each pivot,
you can override on_pivot_start():
>>> from lemke.lemke import LemkeCallback
>>> class LogEnteringVariables(LemkeCallback):
... def on_pivot_start(self, tableau, leave, enter):
... print(enter)
>>> result = runlemke(lcp=problem, callback=LogEnteringVariables())
z0
z2
z1
Note
No callbacks are invoked in the trivial case \(q \ge 0\).