Sub optimal algorithm for linear programming

Question

Consider the linear programming problem \begin{align} f^* = \max_{x}&~p^Tx~\\~&A^Tx\leq b~,~0\leq x_i\leq 1 \end{align}where $c$ is a $n\times 1 $ vector, $A$ is a $n\times c$ matrix and $b$ is a $c \times 1$vector. Here $x=[x_1,\dots,x_n]$ is the $n \times 1$ vector to be found. Thus, in addition to the box constraint $0\leq x_i\leq 1$, we have $c$ constraints. It is also known that $p,A,b$ are all element-wise positive. Assume that this problem is feasible and let $(x^*,\lambda^*)$ be the optimal solutions for the primal and dual respectively. Note that $\lambda^*$ is a $c\times 1$ vector.

Consider the function \begin{align} d(\lambda)=\max_{x}(p-A\lambda)^Tx+\lambda^Tb~,~0\leq x_i \leq 1 \end{align} where $\lambda$ is a $c \times 1$ vector. It is not hard to see that $d(\lambda^*)=f^*$. Define the vector $r$ such that $$r_i=[p-A\lambda^*]_i$$ where $[]_i$ denotes the $i^{th}$ entry.

Thus \begin{align} d(\lambda^*)=\max_{x}r^Tx+\lambda^Tb~,~0\leq x_i \leq 1 \end{align} Now, consider the following strategy constructing the $n\times 1$ vector $y$ \begin{align} y_i = \begin{cases}1 ~,& r_i>0 \\0~,& r_i \leq 0\end{cases} \end{align} You can see the motivation of defining $y$ from $d(\lambda^*)$. Given this,

• Can $y$ violate the constraints?
• Can we comment on the gap $$f^*-p^Ty$$ in terms of $p,A,b,n,c$?

if anyone is interested, background: This is from an engineering problem. Typically $c$ won't be more than 2, thus that many constraints. $n$ will run into millions of variables. $p,A,b,n,c$ comes from our data. In some sense, finding $\lambda^*$ is an easier problem for us. And also from an engineering perspective not related to linear programming, $y$ is much easier to implement in our system rather than solving a linear programming system. Often, empirically, for our data (for smaller sampled sets), we have seen that it does well also. Is there any justification?

Answer 1

Yes, $y$ can violate the constraints. The one about box constraints ($0 \leq y \leq 1$) is OK by how $y$ is defined, the problem is the other one ($A^T y \leq b$). First I tried getting a feasibility proof unsuccessfully, then I decided to build a reproducible counterexample using Python:

import cvxpy as cp
import numpy as np

# Generate a random instance according to the conditions stated by OP (small c, larger n).
n = 150
c = 2

# Build p, A and b, which are elementwise positive.
np.random.seed(1)
p = np.random.uniform(0, 1, n)
A = np.random.uniform(0,1, (n,c))
b = np.random.uniform(0,1, c)

# Define and solve the CVXPY problem. 
# Note that we'll solve the primal problem only to get the dual variable lambda^*
x = cp.Variable(n)
prob = cp.Problem(cp.Maximize(p.T@x),
                 [A.T @ x <= b, 
                  x <= 1,
                  x >= 0 ])
prob.solve()

print("A dual solution is")
print(prob.constraints[0].dual_value)

# Construct the suggested heuristic solution, starting from the dual solution:
lambda_opt = prob.constraints[0].dual_value
r = p - np.matmul(A, lambda_opt)
y = (r > 0).astype(int)

# Verify feasibility for y:
print("Check nonpositive:", (A.T @ y) -b )

This piece of code outputs the array [0.138349, 0.11109537] at the last line, therefore, $A^T y - b \leq 0$ doesn't hold and $y$ is unfeasible. As a conclusion, without further information or hypotheses about $A, b$ and $p$, the proposed approach won't work in general.

Regarding the second part of your question, it reminds me about Approximation Algorithms. I cite from the Wikipedia article:

In computer science and operations research, approximation algorithms are efficient algorithms that find approximate solutions to optimization problems (in particular NP-hard problems) with provable guarantees on the distance of the returned solution to the optimal one.

Assuming the problem of finding $\lambda^*$ easier as you mention, there are two scenarios:

1) If you can find $\lambda^*$ by other means and not solving the linear problem nor its dual, then you could use complementary slackness conditions to try find $x^*$. This somewhat relates to what you are building with the binding/nonbinding dual constraints (i.e. $r_i$ values). There are several examples available on how to procede, like here, here or here.

2) If the dual problem is "easier" to solve, your approach sounds similar to using a primal-dual method. You can read an explanation on how it works at Section 5.1 of this document, while the original source where it was proposed for solving linear programs is several decades old.¹

The general idea will be as follows:

1. Find a feasible dual solution $\lambda$.
2. Given $\lambda$, find some $x$ that minimizes the violation of complementary slackness in the primal problem. (This is the step that reminds me of what you are trying to approach when constructing $y$).
3. If complementary slackness holds, $y$ is optimal, and the algorithm terminates.
4. Otherwise, change $\lambda$ so as to improve the dual objective, and go to 2.

---

¹ Dantzig, George Bernard, Lester Randolph Ford Jr, and Delbert Ray Fulkerson. A PRIMAL--DUAL ALGORITHM. No. P-778. RAND CORP SANTA MONICA CA, 1956.