非负二次规划的乘性更新法则

本文原文为《Multiplication Updates for Nonnegative Quadratic Programming》，第一作者为Fei Sha，他是大牛Michael I. Jordan的学生。

由于本文有三十多页且算法收敛性的证明篇幅较长、公式较多在此只简单介绍算法部分。

一、介绍

在神经计算与统计学习问题中，往往会出现非负的约束。比如最大间隔分类器SVM、贝叶斯网络中的密度估计、非负矩阵分解中的降维以及声学中回声的消除技术。在这些优化求解中是很难有直接的一个解析式能一步就算出最优值的，而这时候就会出现很多迭代算法，一步步地让迭代出来的值靠近最真实的解。

比如对于最小化目标函数F(v)的求解中，最简单一般的算法便是加性迭代梯度下降（GD，Gradient Descent）：

（1）

然而上式会出现不满足非负限制的情况，这时候我们可以将迭代法则改为：

（2）

也就是说，只要一出现为负的情况，就全部取为0。注意到上俩式中的学习速率要启发式地设置为正数，而其大小会让求解过程的速度变慢（当太小）甚至会影响到更新能收敛还是在解附近震荡（当太大）。所以我们称其为启发式参数设置，就是说在各种问题中没有一个统一标准，要根据实际问题设置其大小。

对于非负的约束，有一个更为合适的更新法则叫做指数梯度（EG，Exponentiated Gradient ）为：

（3）

上式我们可以看到乘子是为非负的，只要初始值为非负，则可以保证整个迭代更新过程的非负性。若将两边同时取对数。则更新变为：

（4）

值得注意的是，如果优化问题的解是稀疏的，那么乘性更新求解的速度比加性更新要快很多。

在本文中，提出了一种最优解是被限制在非负的轴对齐区域的二次规划凸问题的乘性更新。这个更新规则虽然简单，但是能使目标函数值下降并单调收敛至全局最优点。

现在给出本文问题的形式，带非负约束的二次规划问题：

（5）

此处假设矩阵A为对称以及严格正定的，所以公式（5）是有界且其最优化是凸问题。特别地，它只有一个全局最优解而没有局部的最小解。针对此问题，已有简单的更新法则为：

（6）

此项更新有一个默认条件就是A中的元素必须是为非负的。然而实际中A不一定是非负的，公式（6）的分母在为负的情况下对于v的更新便会违反了约束规则。本文将公式（6）推广至范围更广的非负二次规划问题的求解中，提出新的更新规则。同时证明了该更新法则能单调收敛于全局最优点。事实上，很多算法收敛性的证明都是依赖于对辅助函数的构造，如期望最大化算法（EM）以及非负矩阵分解算法（NMF）等，该文也是如此。但是单纯的辅助函数构造还只能证明更新收敛于一个局部的静态点。于是该文还运用了一个特殊的结构来证明其收敛于全局最优点。下面会介绍其算法。

二、算法

本文只假设A为半正定矩阵。对于A中可能同时存在正负元素的情况，我们将其分解和为两部分：

和 （7）

于是有，目标函数也被分解为三个部分：

（8）

其中：

，， （9）

我们记：

，， （10）

由此可以得到的一个乘性更新法则如下：

（11）

说到这里就可以直接上程序。

clear all
clc

M = 30;
NumoIter = 100;

%create A that is symmetric and semipositive

 a = 10*(rand(M,M) -0.5.*ones(M,M));
 A = a'*a;
% eig(A)  if A is symmetric and semipositive 
% all the elements of eig(A) should be positive

% elementes of A are bounded in [-0.5,0.5] both negative and positive
index_positive = (A>0);
index_negative = (A<0);
Apositive = A.*index_positive;
Anegative = -A.*index_negative;

% rand('state',100); 
b = 10*(rand(1,M) - 0.5*ones(1,M));
% elementes of b are bounded in [-0.5,0.5]both negative and positive
% rand('state',1001); 
v = abs(rand(M,1));

F_v=[];
for iter = 1:NumoIter
   a = Apositive*v;
   c = Anegative*v;
    for i = 1:M
    v(i) = v(i)*(-b(i)+sqrt(b(i)^2+4*a(i)*c(i)))/(2*a(i));
    end
    
    F_v(iter) = (1/2)*v'*A*v + b*v;
end

% figure;
plot(F_v);
grid on;
xlabel('Number of iterations');
ylabel('Values of cost function');

其中，有需要注意的是矩阵A的生成。之前不懂什么叫A是半正定的，以为是A中的元素要大于等于0就可以（太年轻了）。生成的A矩阵去算总是看着代价函数乱跑。查了一下才知道半正定不是这个意思。关于正定与半正定可参见维基百科。其中提到“Inlinear algebra, asymmetricn×nrealmatrixMis said to bepositive definiteifz^TMzis positive for any non-zero columnvectorzofnreal numbers. Herez^Tdenotes thetransposeofz.”“For any real matrix,is a positive semi-definite matrix.” 程序中对于A的生成就是采用此方法。