Минимизация нормы L1 Ax-b с помощью scipy.optimize.linprog

Известно, что задачу минимизации 1-нормы Ax - b над x можно записать как задачу линейного программирования. Я хотел бы использовать scipy.optimize.linprog для решения этой проблемы минимизации, но у меня возникли проблемы с интерпретацией выходного массива linprog. Ниже приведена реализация моего кода метода решения Ax - b с помощью линейного программирования:

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)

Иногда это работает. Например, в следующем случае, когда минимальная 1-норма равна 0, мы имеем:

#
#    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))

производит продукцию

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

Однако в следующей задаче, где нет точного решения, не создается вектор с минимальной 1-нормой:

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)}')

с выходом

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

Поэтому моя реализация лишь иногда дает правильный ответ. Эмпирически я обнаружил, что для систем, в которых существует точное решение, ответ правильный тогда и только тогда, когда все ненулевые записи выходных данных op_result.x находятся в первых Ncol записях. Когда это не так, кажется, что нет никакого способа получить решение от op_result.x; Я попробовал, например, переставить все записи и посмотреть, дадут ли переставленные параметры решение. Мой вопрос: почему моя реализация работает только иногда и что я могу изменить, чтобы она работала всегда?