线性方程的非负解 - 卤蛋的自我修养

本文最后更新于：2021年12月15日中午

写在前面

最近半年没有写博客，最近笔记都用Notion做的，后面陆续会转移到博客里来。本篇博客记录大学同学的一个问题『求线性方程组非负解』，算法流程参考文献^[1]。

问题回顾

求方程组\(A_{m\times n}x_{n\times1}=b_{m\times1}\)，满足解向量\(x\)元素非负，即\(\forall i,x_i\geq0\)。

这类问题可转为非负最小二乘曲线拟合问题的形式 \[ \min_x\|Ax-b\|_2^2,\quad\text{s.t.,}\quad x\geq0 \] 这种拟合问题已经有成熟的工具包来求解，例如 MATLAB 的lsqnonneg^[2]和 Python 的scipy.optimize.nnls^[3]。判断解的存在性也很简单，因为最小二乘拟合的结果带入原方程得到残差

存在非负解：\(\|A\hat x-b\|_2=0\)
不存在非负解：\(\|A\hat x-b\|_2>0\)

注意：计算过程不可避免地会产生截断误差，所以得到了一个非常小但不严格为\(0\)的结果。

另外，求非负解问题也可转为线性规划问题 \[ \min_x \sum_i x_i,\quad\text{s.t.,}\quad Ax=b, x\geq0 \] 估计结果\(\hat x\)的分量\(\sum_i \hat x_i\)也可以作为是否存在非负解的判断条件

存在非负解：\(\sum_i \hat x_i>0\)
不存在非负解：\(\sum_i \hat x_i=0\)

然而这些方法都只能在成熟的数值计算软件中完成，对小白不友好，下面给一个初级的方法。

线性方程组求非负解

step 1: 得到增广矩阵\(M=(A_{m\times n}|b_{m\times 1})\in\mathbb R_{m\times (n+1)}\)
step 2: 对增广矩阵使用 Gauss 消去法，简化行阶梯形矩阵得到\(M^0\)，并判断方程组是否有解。
step 3: 进入循环(\(k=0,1,\ldots,k_{\max}\))，处理矩阵\(M^k\)的最后一列\(M_{(:,n+1)}^k\)
- 如果最后一列\(M_{(:,n+1)}^k\)存在负元素，不妨假设第\(i\)个元素\(M_{(i,n+1)}^k<0\)，判断矩阵\(M^k\)的第\(i\)行前\(n\)列是否存在负元素，即\(M_{(i,j)}^k<0,\forall j=1,2,\ldots,n\)
  - 如果对所有\(j=1,2,\ldots,n\)，都有\(M_{(i,j)}^k>0\)，即 \[ M_{(i,1)}^k y_1+\cdots+M_{(i,n)}^k y_n = M_{(i,n+1)}^k<0 \] 根据非负解的非负组合必非负，此方程无非负解。停止循环，输出原方程无非负解。
  - 如果第\(i\)行前\(n\)列存在负元素，不妨记第\(j\)个元素小于\(0\)，即\(M_{(i,j)}^k<0\)，注意\(1\leq j\leq n\)。对矩阵\(M^k\)的第\(i\)行第\(j\)列使用 Gauss 消去法，保证矩阵\(M^k\)的第\(i\)行第\(j\)列为\(1\)，第\(j\)列的其他行均为\(0\)，即 \[ M_{(k,j)}^k=\begin{cases} 1& k=i\\ 0& k\neq i\\ \end{cases} \] 如果第\(i\)行前\(n\)列存在多个负元素，取其中任意一个为基准进行简化行阶梯形矩阵。
- 如果最后一列\(M_{(:,n+1)}^k\)所有元素都非负，说明方程存在非负解，拆分增广矩阵\(M^k=(\hat A|\hat b)\)，其中向量\(\hat b_i\geq0,\forall i\)。找到矩阵\(\hat A\)的负元素对应的列，这些列迫使解向量对应元素为为\(0\)，这些分量称为自由变量。带入原方程后得到规模更小的方程，利用\(\hat b\)反解非自由变量，从而得到非负解。

下面给一个例子来解释自由变量和非自由变量。 \[ \begin{pmatrix} -1&0&1\\ 1&1&0\\ \end{pmatrix} \begin{pmatrix} x_1\\x_2\\x_3 \end{pmatrix}= \begin{pmatrix} 1\\0 \end{pmatrix} \] \(\hat A\)非负元素仅出现在第一列，对应第一个分量\(x_1=0\)为自由变量。带入原方程得 \[ \begin{pmatrix} 0&1\\ 1&0\\ \end{pmatrix} \begin{pmatrix} x_2\\x_3 \end{pmatrix}= \begin{pmatrix} 1\\0 \end{pmatrix} \] 解得非自由变量\(x_2=0,x_3=1\)。

注记：

自由变量如果不取\(0\)，而取正数，例如取\(x_1=0.1\)，带入得 \[ \begin{pmatrix} -1&0&1\\ 1&1&0\\ \end{pmatrix} \begin{pmatrix} 0.1\\x_2\\x_3 \end{pmatrix}= \begin{pmatrix} 1\\0 \end{pmatrix} \] 出现\(x_2=-0.1,x_3=1.1\)，即无法保证非自由变量非负性，即使\(x_1\)取得再小。
非自由变量如果出现某一列均为0的方程，例如 \[ \begin{pmatrix} 0&1\\ 0&0\\ \end{pmatrix} \begin{pmatrix} x_2\\x_3 \end{pmatrix}= \begin{pmatrix} 1\\0 \end{pmatrix} \] 则对应方程无数解，即\(x_2\geq0，x_3=1\)。
步骤2中判断方程是否存在解其实是有必要的，或者说方程是否存在矛盾方程，例如 \[ \begin{pmatrix} 1&0\\ 0&0\\ \end{pmatrix} \begin{pmatrix} x_2\\x_3 \end{pmatrix}= \begin{pmatrix} 1\\1 \end{pmatrix} \] 这个可以从阶梯型矩阵很容易判断出来，方程无解必然无法求非负解。
循环迭代数\(k_{\max}\)如何设置，迭代是否存在无意义的循环，这些都可以深入思考，正如文献^[1]里提到：如果\(M^k\)第\(i\)行前\(n\)列存在多个负元素的情况：

If this happens, just make sure that pick your pivot entries so that you’re not getting into a “loop” – i.e. continually repeating the same set of three pivots again and again and again.

MATLAB 实现

clear; clc;

% 参考文献
% http://web.math.ucsb.edu/~padraic/caltech/math1b_2011/math1b_2011_recitation_wk2.pdf
%% case 1
% 原例子
A = [1,1,1,1,1,1,1,1;-1,1,2,2,1,-1,-2,-2;2,2,1,-1,-2,-2,-1,1];
b = [1,-3,0]';
disp('====== case 1 ======');

x_sol = nonsol(rref([A,b]));
disp(['error_nonsol = ', num2str(norm(A*x_sol-b))]);

% MATLAB 标准解方程
x_reg = A\b;
disp(['error_stand = ', num2str(norm(A*x_reg-b))]);

% MATLAB 方程非负解
x_lsq = lsqnonneg(A,b);
disp(['error_lsqnon = ', num2str(norm(A*x_lsq-b))]);

%% case 2
% 文献例子1
A = [1,1/3,4/3,1/3;0,-1,-6,2];
b = [5/3,-12]';
disp('====== case 2 ======');

x_sol = nonsol(rref([A,b]));
disp(['error_nonsol = ', num2str(norm(A*x_sol-b))]);

% MATLAB 标准解方程
x_reg = A\b;
disp(['error_stand = ', num2str(norm(A*x_reg-b))]);

% MATLAB 方程非负解
x_lsq = lsqnonneg(A,b);
disp(['error_lsqnon = ', num2str(norm(A*x_lsq-b))]);
%% case 3
% 文献例子2
A=[1,2,3;0,2,6;0,-1,-3];
b = [2,4,-2]';
disp('====== case 3 ======');

x_sol = nonsol(rref([A,b]));
disp(['error_nonsol = ', num2str(norm(A*x_sol-b))]);

% MATLAB 标准解方程
x_reg = A\b;
disp(['error_stand = ', num2str(norm(A*x_reg-b))]);

% MATLAB 方程非负解
x_lsq = lsqnonneg(A,b);
disp(['error_lsqnon = ', num2str(norm(A*x_lsq-b))]);
%% case 4
% 随机矩阵A和向量b；
% A = rand(4,4);
% b = rand(4,1);
% case 4
% 设置矩阵A和非负解，计算b = A*x，利用A和b反解x；
A = rand(4,4);
x = rand(4,1);
b = A*x;
disp('====== case 4 ======');

x_sol = nonsol(rref([A,b]));
disp(['error_nonsol = ', num2str(norm(A*x_sol-b))]);

% MATLAB 标准解方程
x_reg = A\b;
disp(['error_stand = ', num2str(norm(A*x_reg-b))]);

% MATLAB 方程非负解
x_lsq = lsqnonneg(A,b);
disp(['error_lsqnon = ', num2str(norm(A*x_lsq-b))]);

%% nonsol function
function x_sol = nonsol(Ab)
itermax = 50;
AA = cell(1,itermax);

for iter = 1 : itermax
    AA{iter} = Ab;
    ind_raw = find(Ab(:,end)<0);
    if isempty(ind_raw)
        break;
    elseif length(ind_raw) > 1
        ind_raw = ind_raw(1);
    end
    ind_col = find(Ab(ind_raw,:)<0);
    if ind_col(1) == size(Ab,2)
        disp('no nonnegative solution is possible');
        break;
    else
        ind_col = ind_col(1);
    end
    Ab(ind_raw,:) = Ab(ind_raw,:)/Ab(ind_raw,ind_col);
    k = 1;
    for k = 1:size(Ab,1)
        if k == ind_raw
            continue;
        end
        Ab(k,:) = Ab(k,:)-Ab(k,ind_col)/Ab(ind_raw,ind_col)*Ab(ind_raw,:);
    end
    if iter == itermax
        disp('Reach the maximum number of iterations');
    end
end

x_sol = zeros(size(Ab,2)-1,1);
ind = 1:size(Ab,1);
ind = ind(sum(Ab < 0,2)==0);
x_sol(ind) = Ab(:,ind)\Ab(:,end);

end

输出结果

====== case 1 ======
no nonnegative solution is possible
error_nonsol = 2.7331
error_stand = 4.9651e-16
error_lsqnon = 0.44721
====== case 2 ======
no nonnegative solution is possible
error_nonsol = 12.1152
error_stand = 2.2204e-15
error_lsqnon = 0.97619
====== case 3 ======
error_nonsol = 0
警告: 矩阵为奇异工作精度。  
error_stand = NaN
error_lsqnon = 0
====== case 4 ======
error_nonsol = 2.2204e-16
error_stand = 2.2888e-16
error_lsqnon = 2.2888e-16

迭代得到的结果与lsqnonneg基本一致，验证方法的有效性。另外，这个简单的想法很容易用其他语言重新编译一遍。

小结

这个初级想法^[1]的核心是通过初等行变换使得增广矩阵最后一列均非负，自由变量设为\(0\)即可消除方程负元素的影响，从而保证解也是非负的。

References

线性代数

解方程非负解 MATLAB

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议，转载请注明出处！

🔒在线主成分分析下一篇