Minimizing the L1 norm of Ax - b using scipy.optimize.linprog

Question

It is known that the problem of minimizing the 1-norm of Ax - b over x can be written as a linear programming problem. I would like to use scipy.optimize.linprog to solve this minimization problem, but I am having trouble interpreting the output array of linprog. Below is my code implementation of the method to solve Ax - b via linear programming:

import numpy as np
import scipy.optimize as op
from scipy.linalg import lstsq

def build_A_ub(A):
    Nrow, Ncol = A.shape
    I          = np.identity(Nrow)

    return np.block([
        [A , -I],
        [-A, -I]
    ])

def build_b_ub(b):
    return np.concatenate([b, -b])

def build_c(A):
    Nrow, Ncol = A.shape
    return np.concatenate([np.zeros(Ncol), np.ones(Nrow)])

def obtain_l1norm_minimizer(A, b):
    """
    return op_result which contains parameters x which minimizes |Ax - b|_1
    """
    A_ub, b_ub, c = build_A_ub(A), build_b_ub(b), build_c(A)
    
    return op.linprog(c, A_ub=A_ub, b_ub=b_ub)

This sometimes works. For example, in the following case where the minimal 1-norm is 0, we have:

#
#    example implementation of L1 norm minimization of Ax - b for n x m matrix 
#    with n < m and rank(n) (i.e: exact solution to Ax - b)
#

from numpy.random import uniform

np.random.seed(2)
Nrow = 3
Ncol = 4
A    = uniform(-1, 1, [Nrow,Ncol])
b    = uniform(-1, 1, [Nrow])

assert np.linalg.matrix_rank(A) == min([Nrow, Ncol])

op_result  = obtain_l1norm_minimizer(A, b)
difference = A @ op_result.x[:Ncol] - b

print("Solution vector (x,y):")
print(np.round(op_result.x, 6))
print("Ax - b:")
print(np.round(difference, 6))

produces output

Solution vector (x,y):
[1.366741 0.373728 0.       1.557993 0.       0.       0.      ]
Ax - b:
[ 0. -0.  0.]

However, in the following problem where there is no exact solution, it does not produce a vector with minimal 1-norm:

np.random.seed(2)
Nrow = 6
Ncol = 3
A    = uniform(-1, 1, [Nrow,Ncol])
b    = uniform(-1, 1, [Nrow])

assert np.linalg.matrix_rank(A) == min([Nrow, Ncol])

op_result        = obtain_l1norm_minimizer(A, b)
x_lstsq, _, _, _ = lstsq(A, b)

# compare 1-norm of original, solution obtained from 1-norm minimization, and least squares solution

print(f'original 1-norm : {np.linalg.norm(b, 1)}')
print(f'1-norm from 1-norm minimization: {np.linalg.norm(A @ op_result.x[:Ncol] - b, 1)}')
print(f'1-norm from least-squares minimization : {np.linalg.norm(A @ x_lstsq - b, 1)}')

with output

original 1-norm : 3.364444701650462
1-norm from 1-norm minimization: 3.2226922553095174
1-norm from least-squares minimization : 2.374926794178961

Therefore, my implementation only sometimes give the right answer. Empirically, for systems in which there is an exact solution, I find that the answer is correct if and only if all the non-zero entries of the output op_result.x are in the first Ncol entries. When this isn't true, there doesn't seem to be any method to obtain a solution from op_result.x; I tried, for example, permuting all the entries and seeing if the permuted parameters generated a solution. My question is: why is my implementation only sometimes working and what can I modify to make it always work?

Answer 1

Alternatively, you can remove the bounds from your decision variables. By default, they must be non-negative. Presumably you intended for them to be unbounded, though, so:

op.linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=(None, None))

I haven't checked your code carefully since there is a better way to solve such a problem, but this seems to work.