线性方程的非负解

本文最后更新于: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


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