对于 L1 norm 问题，我们通常先将其转化为线性规划，然后对线性规划问题进行求解。我们之前使用的是内点法，这里我们测试另一种方法——对偶单纯形法。

1. 问题说明

我们要求解一个 L1-norm regression problem:

\[\min _{\mathbf{x}}\|A \mathbf{x}-\mathbf{b}\|_{1}\]

这里的 $A \in \mathcal{R}^{m \times n}$ 和 $b \in \mathcal{R}^{m}$ 都是给定的，并且有 $m>n$.

我们的目标是寻找 $\mathbf{x} \in \mathcal{R}^{n}$ 使得 L1-norm

\[\|A \mathbf{x}-\mathbf{b}\|_{1}\]

最小. 对于向量 $\mathbf{y}=\left(y_{1}, \cdots, y_{m}\right)^{\top}$, L1-norm 定义为

\[\|\mathbf{y}\|_{1}=\sum_{j=1}^{m}\left|y_{j}\right|\]

2. 对偶问题

上面的优化问题可以转化成线性规划的形式:

\[\begin{aligned} &\min _{\mathbf{x}, \mathbf{u}, \mathbf{v}} \quad \mathbf{e}^{\top}(\mathbf{u}+\mathbf{v}) \\ &\text { s.t. } \quad A \mathbf{x}-\mathbf{b}=\mathbf{u}-\mathbf{v}, \quad \mathbf{u}, \mathbf{v} \geq \mathbf{0} \end{aligned}\]

这里 $\mathbf{e}=(1, \cdots, 1)^{\top} \in \mathcal{R}^{m}$. 接下来我们要考虑它的对偶问题, 我们先构造 Lagrange 函数:

\[\begin{aligned} L(\mathbf{x},\mathbf{u},\mathbf{v},\lambda,\sigma_1,\sigma_2)&= \mathbf{e}^{\top}(\mathbf{u}+\mathbf{v})+\lambda^{\top}(A \mathbf{x}-\mathbf{b}-\mathbf{u}+\mathbf{v})-\sigma_1^{\top}\mathbf{u} -\sigma_2^{\top}\mathbf{v}\\ &=\lambda^{\top}A\mathbf{x}+(\mathbf{e}-\lambda-\sigma_1)^{\top}\mathbf{u}+(\mathbf{e}+\lambda-\sigma_2)^{\top}\mathbf{v}-\lambda^{\top}\mathbf{b} \end{aligned}\]

其拉格朗日对偶函数为:

\[g(\lambda,\sigma_1,\sigma_2) = \inf_{\mathbf{x},\mathbf{u},\mathbf{v}} L(\mathbf{x},\mathbf{u},\mathbf{v},\lambda,\sigma_1,\sigma_2)= \begin{cases} -\mathbf{b}^{\top}\lambda, & A^{\top}\lambda=0,\sigma_1=\mathbf{e}-\lambda,\sigma_2=\mathbf{e}+\lambda\\ -\infty,& \text{else} \end{cases}\]

我们进行变量代换 $\mathbf{y}=-\lambda$ 后可得其对偶问题:

\[\begin{aligned} \max_{\mathbf{y}} &\quad \mathbf{b}^{\top} \mathbf{y} \\ \text { s.t. } &\quad A^{\top} \mathbf{y}=0 \\ &\quad -\mathbf{e} \leq \mathbf{y} \leq \mathbf{e} \end{aligned}\]

对于线性规划而言，原始问题与对偶问题是完全等价的，因此我们可以通过求解对偶问题获得原始问题的解.

3. 代码实现

clc,clear;
m = 10;
n = 5;
A = randn(m,n);
b = randn(m,1);
[data, info] = OneNormLP(A, b)


function [data, info] = OneNormLP(A, b)
    MAX_iter = 100;
    iter = 0;
    info.run = "Failure";
    [m,n] = size(A);
    Idx = 1:m;
    B = 1:n;
    M = inv(A(B,:));
    while(iter < MAX_iter)
        iter = iter+1;
        B_c = setdiff(Idx,B);
        A_ = A(B,:);
        b_ = b(B,:);
        if abs(det(A_))<1e-6
            info.msg = "Non-degeneracy assumptions are not met.";
            return
        end
        x = M*b_;
        data.obj = norm(A*x-b,1);
        data.x = x;
        data.loop = iter;
        h = A*x-b;
        A_c = A(B_c,:);
        h_c = h(B_c,:);
        y_c = sign(h_c);
        y_ = -A_'\(A_c'*y_c);
        [val, s] = max(abs(y_));
        if val<=1 || data.obj<1e-6
            info.run = "Success";
            return
        end
        j = B(s);
        e_s = zeros(n,1);
        e_s(s) = 1;
        t_c = -sign(y_(s)) .* y_c .* (A_c*(A_\e_s));
        ht = abs(h_c)./t_c;
        ht(ht<0) = -1e+9*ht(ht<0);
        [~, r] = min(ht);
        u = zeros(n,1);
        u(s) = 1;
        v = A(B_c(r),:)' - A(B(s),:)'; 
        M = M - (M*u*v'*M)/(1+v'*M*u);
        B(s) = B_c(r);
    end
    info.msg = "Maximum number of iterations exceeded.";
end