In this post we don't even consider optimization, but restrict to elementary linear algebra. Suppose we want to find a solution to a system of linear equations:
\begin{array}{llllllllllll} & 2x_1 & - & x_2 & + & x_3 & - & 3x_4 & + & x_5 & = & 1, \\ & & & x_2 & -& 2x_3 & & & +& x_5 & = & -1,\\ & 3x_1 & & & - & x_3 & - & x_4 & & & = & 2, \\ & 2x_1 & + & x_2 & - & 3x_3 & - & 3x_4 & + & 3x_5 & = & -0.5. \end{array}
This system of linear equations is in fact infeasible (has no solution). One way to see it is to multiply the equations by coefficients given on the right-hand side and add them up:
\begin{array}{lllllllllllll} & 2x_1 & - & x_2 & + & x_3 & - & 3x_4 & + & x_5 & = & 1 & / \cdot (-1) \\ & & & x_2 & -& 2x_3 & & & +& x_5 & = & -1 & / \cdot (-2) \\ & 3x_1 & & & - & x_3 & - & x_4 & & & = & 2 & / \cdot 0 \\ & 2x_1 & + & x_2 & - & 3x_3 & - & 3x_4 & + & 3x_5 & = & -0.5 & / \cdot 1 \\ \hline & & & & & & & & & 0 & = & 0.5. & \end{array}
We get an obvious contradiction, which proves the system has no solution. The vector
y = [-1, -2, 0, 1]
of weights used above is therefore a proof (certificate) of infeasibility. MOSEK produces such a certificate automatically. Here is a simple script in MATLAB that computes precisely that vector:
As output we get:
Here are a few comments:
In the example above y_3=0 which means that the third equation does not matter: infeasibility is caused already by some combination of the 1st, 2nd and 4th equation. In many practical situations the infeasibility certificate y will have very few nonzeros, and those nonzeros determine a subproblem (subset of equations) which alone cause infeasibility. The user can configure MOSEK to print an infeasibility report, which in our example will look like:
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
prob = []; | |
prob.a = [ 2 -1 1 -3 1; | |
0 1 -2 0 1; | |
3 0 -1 -1 0; | |
2 1 -3 -3 3 ]; | |
prob.buc = [ 1 -1 2 -0.5]'; | |
prob.blc = [ 1 -1 2 -0.5]'; | |
[r, res] = mosekopt('minimize', prob); | |
res.sol.bas.y |
As output we get:
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
ans = | |
-1 | |
-2 | |
0 | |
1 |
Here are a few comments:
- General theory guarantees that the dual variable y is a convenient place that can store the certificate.
- When a system of equations has no solution then an appropriate certificate is guaranteed to exist. This is a basic variant of Farkas' lemma, but it should be intuitively clear: your favorite method of solving linear systems (for instance Gaussian elimination) boils down to taking successive linear combinations of equations, a process which ends up either with a solution or with an "obvious" contradiction of the form 0=1.
Using the certificate to debug a model
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
MOSEK PRIMAL INFEASIBILITY REPORT. | |
Problem status: The problem is primal infeasible | |
The following constraints are involved in the primal infeasibility. | |
Index Name Lower bound Upper bound Dual lower Dual upper | |
0 1.000000e+00 1.000000e+00 0.000000e+00 1.000000e+00 | |
1 -1.000000e+00 -1.000000e+00 0.000000e+00 2.000000e+00 | |
3 -5.000000e-01 -5.000000e-01 1.000000e+00 0.000000e+00 |
This report is nothing else than a listing of equations with nonzero entries in y, and y is the difference between Dual lower and Dual upper. Analyzing this set (which hopefully is much smaller than the full problem) can help locate a possible modeling error which makes the problem infeasible.
To conclude, let us now phrase the above discussion in matrix notation. A linear equation system
Ax=b
is infeasible if and only if there is a vector y such that
A^Ty = 0 \quad \mbox{and}\quad b^Ty \neq 0.
The situation is analogous for linear problems with inequality constraints. We will look at an example of that in another blog post.