使用 scipy.optimize.linprog 最小化 Ax - b 的 L1 范数

Question

众所周知，最小化

Ax - b

超过

x

的 1-范数问题可以写成线性规划问题。我想使用

scipy.optimize.linprog

来解决这个最小化问题，但我在解释

linprog

的输出数组时遇到困难。下面是我通过线性规划解决

Ax - b

的方法的代码实现：

import numpy as np
import scipy.optimize as op

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)

这有时有效。例如下面的代码：

#
#    example implementation of L1 norm minimization of Ax - b for n x m matrix 
#    with n > m and rank(m) (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.]

但是，使用随机种子

0

，我得到了

Solution vector (x,y):
[0.       0.097173 0.       1.050148 0.       0.       0.895963]
Ax - b:
[0.       0.       0.895963]

因此，我的实现有时只能给出正确的答案。根据经验，我发现当且仅当输出

op_result.x

的所有非零条目都在第一个

Ncol

条目中时，答案才是正确的。当这不是真的时，似乎没有任何方法可以从

op_result.x

获得解决方案；例如，我尝试排列所有条目并查看排列后的参数是否生成解决方案。我的问题是：为什么我的实现只是有时有效？我可以修改什么以使其始终有效？

Answer 1

由于方程少于未知数，因此您可以使用

scipy.linalg.lstsq

来修改它以使其始终有效。

import numpy as np
from numpy.random import uniform
from scipy.linalg import lstsq
np.random.seed(0)
Nrow = 3
Ncol = 4
A = uniform(-1, 1, [Nrow,Ncol])
b = uniform(-1, 1, [Nrow])
x, _, _, _ = lstsq(A, b)
print(A@x - b)
# [-2.77555756e-17  3.33066907e-16 -2.22044605e-16]