[ Update: question solved, answer is at the end of the post]
I need to solve this problem
$$ \min \| x \|_{\infty} \; \text{s.t.} \; Ax = b $$
where $Ax = b$ is an underdetermined system, e.g., $A \in \mathbb{R}^{n \times m}$, with $n < m$. Unfortunately this problem seems to be not so treated in the research literature. After some googling, I found this thesis by Earle AC to be a clear reference. Many approaches are provided, unfortunately with no code at all, and the most pragmatic way to solve such a problem would be to cast this problem as a linear programming one, i.e.,
$$\min c^T [x |s] \\ A x = b \; , \; x_i \leq s , \; -x_i \leq s\\ $$ where $s$ is meant to be $\| x \|_{\infty}$, $c = [\textbf{0}|1]$, $\textbf{0} \in \mathbb{R}^m$ being the zero vector.
Such an approach is well described in page 36 of thesis by Earle AC
Unfortunately when I try to code this in Python, I obtain not feasible solutions: in particular when using the simplex method of Scipy, at some point in the tableau reduction the algorithm is no more able to find a pivot in the column.
An example is provided at the end of the post.
The questions hence are:
Are there any errors in the code below? Or the problem is effectively ill posed?
Are there any out of the box implemented solvers that deal with the same problem?
Thanks!
# divergence of linear programming for minimization of l_inf norm
# for overdetermined linear solution
import numpy as np
# p.36 of Infinity Norm Minimization EarleAC
n, m = 2,5
np.random.seed(10)
b = np.random.rand(n)
A = np.random.rand(n*m).reshape((n,m))
m = A.shape[1]
# check that system has solution
if (np.linalg.matrix_rank(np.column_stack((A,b))) - np.linalg.matrix_rank(A)):
raise ValueError("matrix rank problem")
# matrix costraints for equality Ax = b
A_eq = np.column_stack((A, np.zeros(n)))
b_eq = b
# matrix constraint for |x|<s
M_3 = np.column_stack((np.diag(np.ones(m)), np.ones(m)))
b_3 = np.zeros(m)
M_4 = np.column_stack(( - np.diag(np.ones(m)), np.ones(m)))
b_4 = np.zeros(m)
# contraint s > 0
M_5 = np.zeros(M_4.shape[1])
M_5[-1] = 1
b_5 = [0.]
A_ub = np.row_stack((M_3, M_4, M_5))
b_ub = np.concatenate((b_3, b_4, b_5))
# define target
from scipy.optimize import linprog
c = np.append(np.zeros(m), [1.])
sol_scipy = linprog(c = c,
A_ub = -A_ub, # minus is to consider <= instead of >=
b_ub = -b_ub, # minus is to consider <= instead of >=
A_eq = A_eq,
b_eq = b_eq,
method = 'simplex'
)
print sol_scipy["message"]
# >> Optimization failed. Unable to find a feasible starting point.
# however a feasible solution exists, given by
x_fea, residuals, rank, __ = np.linalg.lstsq(A, b)
inf_norm_fea = np.max(np.abs(x_fea))
tmp_sol = np.concatenate((x_fea, [inf_norm_fea]))
np.testing.assert_allclose(A_eq.dot(tmp_sol), b_eq)
assert np.all(A_ub.dot(tmp_sol) >= b_ub)
[UPDATE]
Thanks to the suggestion by Michael Grant, solving the problem with cvxpy (readily available out of the box through conda install -c omnia cvxpy) solved the problem in a straightforward manner. Thanks a lot to him (and its amazing job with cvx) and at jjjjjj for suggesting other options. Here is the code that solves the above problem.
import cvxpy as cvx
# Construct the problem.
x = cvx.Variable(c.shape[0])
objective = cvx.Minimize(c * x)
constraints = [A_eq*x == b_eq, A_ub*x >= b_ub ]
prob = cvx.Problem(objective, constraints)
# The optimal objective is returned by prob.solve().
result = prob.solve()
# The optimal value for x is stored in x.value.
print x.value
