使用scipy.optimize.linprog求解Ax-b的L1范数最小化问题异常排查

问题背景 #

将最小化$|Ax - b|_1$的问题转化为线性规划模型后，使用scipy.optimize.linprog实现，但该实现仅在部分场景（如存在$Ax=b$的精确解时）有效；当无精确解时，得到的解的L1范数反而大于最小二乘解，且解的解读存在依赖场景的问题。

原实现代码 #

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：存在精确解（n<m且矩阵满秩） #

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.]

该场景下结果符合预期，$|Ax - b|_1=0$。

场景2：无精确解（n>m且矩阵满秩） #

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

该场景下，L1范数最小化的结果反而差于最小二乘，不符合预期。

问题原因 #

• 变量约束缺失：scipy.optimize.linprog默认所有变量的取值范围是$\geq 0$，但原代码未显式指定变量边界。L1范数最小化中，$x$可以是任意实数，而原实现强制$x$非负，导致仅当最优解$x$为非负向量时结果正确，其他场景下只能找到受限的局部最优解。
• 原线性规划模型的约束逻辑是正确的：通过引入非负变量$y$，将$|Ax - b|_1$的最小化转化为$\sum y_i$的最小化，约束为$|Ax - b| \leq y$，但变量边界的错误限制了$x$的取值空间。

修复后的代码 #

修改obtain_l1norm_minimizer函数，为$x$设置无界边界，$y$保持非负：

def obtain_l1norm_minimizer(A, b):
    """
    return op_result which contains parameters x which minimizes |Ax - b|_1
    """
    Nrow, Ncol = A.shape
    A_ub, b_ub, c = build_A_ub(A), build_b_ub(b), build_c(A)
    
    # 设置变量边界：前Ncol个变量(x)无界，后Nrow个变量(y)非负
    bounds = [(-np.inf, np.inf)] * Ncol + [(0, np.inf)] * Nrow
    
    return op.linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds)

修复结果验证 #

用场景2的代码重新测试，输出结果：

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

此时L1范数最小化的解的误差L1范数小于最小二乘解，符合预期。

内容的提问来源于stack exchange，提问作者Solarflare0