高效的Matlab：mex & profile

04/17/2011 in matlab | Tags: , , , | 2 comments

刚一条句子让wordpress bug了，拆成两篇算了

for的终极取代：多线程混合编程

我们也常会遇到并不能简答并行化的情形，例如针对数组元素的复杂的操作。这时的一个选择就是混合编程: mex, 它可以实现matlab与C/C++/Fortan之间的相互调用。matlab自带的user guide里有详细介绍。

在这里我们讨论两个小trick. 首先注意到向量化的高斯核代码中我们用到了exp(A)，A可能是一个非常大的矩阵。而matlab却会不谙风情的先将这个大矩阵在内存里复制一次，然后再逐一算exp. 有时我们明明算好内存能存下，但matlab却报out of memory了。

我们可以用C来对内存进行直接操作。例如for(i...){A(i)=exp(A(i));}. 不过还有另外一个问题：现在我们计算机大多是多核多线程，matlab也很早就支持了多线程加速，而简单的for只会单线程，所以mex后的代码可能会慢很多。这时的我们的一个选择就是openmp，它为共享内存的多线程提供了非常方便的实现。

下面附上一计算exp(A)的多线程代码，它无需额外内存:

#include "mex.h"
#include
#include
#include
using namespace std;

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray*prhs[] )
{
  double *A, *pp;
  mwSize n, m, p, i, j;

  A = mxGetPr(prhs[0]);
  n = mxGetN(prhs[0]); // # of rows in C
  m = mxGetM(prhs[0]); // # of columns in C
  p = omp_get_num_procs();

#pragma omp parallel for private (j)    \
    default(shared) num_threads(p)
  for (i = 0; i < n; i ++) {
    for (j = 0; j < m; j ++)
      A[i*m+j] = exp(A[i*m+j]);
  }
  return;
}

编译时需要稍修改下编译参数来开启openmp的支持。例如linux下gcc需要加上flag-fopenmp，详见此说明，Windows下的visual studio的相应的flag是/openmp.

但用这段代码时需要很小心。matlab使用了lazy copy的技术，既当我们运行命令A=B时，matlab不会立即将B复制一份给A, 只有当A或者B被修改的时候才执行复制。如果我们用上面代码对B做指数运算，例如mex_exp(B), 这时matlab不知道B变了，所以复制还是不会执行，那么A就会和B一起被修改了，这也许会给后面的代码造成麻烦。

代码究竟慢在哪里

我们可以方便使用tic,toc来知道每段代码的执行时间。但要知道每条语句的详细运行情况的话，非profile莫属。

先运行profile on，再运行需要测试的代码，然后使用profile viewer来查看报告。例如，上面两段不同的计算高斯核的实现的profile报告如下（为了报告更有信息，将其中的一条长语句拆开成了三句）。

根据报告我们就知道每条语句的执行情况，从而可以对其中时间耗费最多的部分进行优化。

无止尽的优化

大多人类活动就是一优化过程。无论是金钱，利益，地位，还是算法，代码，编译，硬件。但除了最优化外，我们还追求一些其他。例如我们用matlab是为了实现和调试上的方便，虽然优化能让计算机做得更快，但是却花去了我们的工作时间，而且使得代码难读懂。糟糕的实现是大忌，过度的优化也不可取。所以，我们需要在计算时间，开发时间，代码的可读性之间做一个权衡。

高效的Matlab：向量化

04/17/2011 in matlab | Tags: , , | 4 comments

前面一直在谈理论，为了避免成为纯理论blog, 这次聊聊最实验的部分：coding.

Matlab是数值计算和仿真中使用最为广泛的软件。它便捷，高效，通常可以直接将伪代码算法转换成matlab代码。但在方便的同时，matlab里面也有很多性能的陷阱。例如令人诟病的for语句。

在matlab里，至少有两种方法将一个数组每个元素加1. 一种来自一般的程序语言for i=1:length(a), a(i)=a(i)+1; end；另一种则借鉴数学记法a=a+1。

但在matlab中，这两实现的性能却有天壤之别。虽然matlab在不断优化for的执行，不过即使是matlab 2010b, 在我的笔记本上（Core 2 Duo 2.4 GHz, 4GB, Win7），对于一个10万元素的一维数组，第一种实现需5微秒，而第二种只用0.4微妙。后者快上10倍。

事实上，用过matlab就会发现，matlab自带的函数通常都有很高的性能。例如刚提到的a=a+1，或是排序 sort(a)，以及令人称道的矩阵乘法，SVD分解等等。事实上，对于两矩阵a和b相乘，如果简单用三层for循环来实现（用C），性能肯定会比matlab差很多。

虽然matlab的界面是用java实现以方便跨平台，在性能上不尽人意，不过内核是用C实现的。可以试试命令matlab -nojvm。对于矩阵运算，matlab是调用高度优化的BLAS/LAPACK库。例如BLAS level3提供的一般矩阵乘法，它将矩阵拆成很多小块，充分利用内存的连续性，cpu的cache架构，多线程等来优化执行。matlab做的只是提供一个方便的人机接口。当然，其他函数例如sort matlab则提供自己优化过的实现。

但为什么用for时就会慢呢。主要原因是matlab是一种脚本语言，通过解释执行。每次我们运行一段代码时，matlab先预解析一遍，生成p-code, 然后一句一句执行。预解析当然是会有代价。在调用for时，matlab会每循环一次就会重新解析下for循环内部的代码，所以自然就慢了。

提高性能的办法是尽量使用matlab提供的函数，而非自己写。例如all（是否元素全非0）, cumsum（累加）, diff（相邻元素的差）,prod（累乘）。另外一些小技巧包括预先通过zeros等操作分配内存。在这里我们着重讨论for的优化，通常称其为向量化（vectorizing）。

用bsxfun来取代for

bsxfun是一个matlab自版本R2007a来就提供的一个函数，作用是”applies an element-by-element binary operation to arrays a and b, with singleton expansion enabled.”

举个例子。假设我们有一列向量和一行向量。

a = randn(3,1), b = randn(1,3)

a =
   -0.2453
   -0.2766
   -0.1913

b =
    0.6062    0.5655    0.9057

我们可以很简单的使用matlab的外乘c=a*b来得到c(i,j)=a(i)*b(j) . 但如果我们想用”外加”呢？既c(i,j)=a(i)+b(j) . 这时我们可以用c=bsxfun(@plus,a,b)来实现。

bsxfun的执行是这样的，如果a和b的大小相同，那么c=a+b. 但如果有某维不同，且a或b必须有一个在这一维的维数为1, 那么bsxfun就将少的这个虚拟的复制一些来使与多的维数一样。在我们这里，b的第一维只有1（只一行），所以bsxfun将b复制3次形成一个3×3的矩阵，同样也将a复制成3×3的矩阵。这个等价于c=repmat(a,1,3)+repmat(b,3,1)。这里

repmat(a,1,3)

ans =
   -0.2453   -0.2453   -0.2453
   -0.2766   -0.2766   -0.2766
   -0.1913   -0.1913   -0.1913

repmat是显式的复制，当然带来内存的消耗。而bsxfun是虚拟的复制，实际上通过for来实现，等效于for(i=1:3),for(j=1:3),c(i,j)=a(i)+b(j);end,end。但bsxfun不会有使用matlab的for所带来额外时间。实际验证下这三种方式

>> c = bsxfun(@plus,a,b)
c =
    0.3609    0.3202    0.6604
    0.3296    0.2889    0.6291
    0.4149    0.3742    0.7144

>> c = repmat(a,1,3)+repmat(b,3,1)
c =
    0.3609    0.3202    0.6604
    0.3296    0.2889    0.6291
    0.4149    0.3742    0.7144

>> for(i=1:3),for(j=1:3),c(i,j)=a(i)+b(j);end,end,c
c =
    0.3609    0.3202    0.6604
    0.3296    0.2889    0.6291
    0.4149    0.3742    0.7144

从计算时间上来说前两种实现差不多，远高于for的实现。但如果数据很大，第二种实现可能会有内存上的问题。所以bsxfun最好。

这里@plus是加法的函数数柄，相应的有减法@minus, 乘法@times, 左右除等，具体可见 doc bsxfun.

下面看一个更为实际的情况。假设我们有数据A和B, 每行是一个样本，每列是一个特征。我们要计算高斯核，既K(i,j)=\exp(-\|A(i,:)-B(j,:)\|^2/\beta . 当然可以用双重for实现（如果第一直觉是用三重for的话…）。

K1 = zeros(size(A,1),size(B,1));
for i = 1 : size(A,1)
    for j = 1 : size(B,1)
        K1(i,j) = exp(-sum((A(i,:)-B(j,:)).^2)/beta);
    end
end

使用2,000×1,000大小的A和B, 运行时间为88秒。

考虑下面向量化后的版本：

sA = (sum(A.^2, 2));
sB = (sum(B.^2, 2));
K2 = exp(bsxfun(@minus,bsxfun(@minus,2*A*B', sA), sB')/beta);

使用同样数据，运行时间仅0.85秒，加速超过100倍。

如要判断两者结果是不是一样，可以如下

assert(all(all(abs(K1-K2)<1e-12)))

在高斯核中，我们通常可以取\beta 是样本之间的平均距离，既\beta=\sum_{i,j=1}^n \frac{\|X(i,:)-X(j,:)\|^2}{(n-1)n} , 此时我们可以使用下面的代码：

n = size(X,1);
S = sum(X);
beta = 2*norm(X,'fro')^2/(n-1) - (S*S')*2/n/(n-1);

各位有兴趣可以比比与for实现的性能差异，呵呵。

漫谈在线学习：阴阳论☯

04/10/2011 in machine learning, online learning | 21 comments

Primal-dual的观点

《红楼梦》第三一回云：天地间都赋阴阳二气所生。世间有阴便有阳，在优化界也是如此。我们将需要最小化的目标函数称之为primal problem. 其对应就有dual problem. 例如常见的Lagrange Duality. 如果primal是凸的，那么dual便是凹的。且在常见情况下，最小化primal等价于最大化dual. 有时候直接求解primal problem很麻烦，但dual却方便很多。或是通过考察dual problem的性质能得到最优解的一些性质。一个经典的案例就是SVM.

上节我们通过子空间投影\Pi_{\mathcal{H}} 来处理罚，这里我们直接考虑最小化损失+罚的形式，既\min_w f(w) + \sum_{t=1}^T\ell_t(w) （罚前的参数\lambda 在这里写进了f(w) 里了）。注意到罚f(w) 的存在，使得目标函数不能如前面那样自然的分成T 块，从而不能直接对每块做梯度下降得到online gradient descent. 于是我们转向dual problem.

Shai Shalev-Shwartz提出可以用Fenchel Duality来研究其对偶式。Shai最近在online/stochastic界相当活跃。虽然和他不熟，不过因为他姓比较长，所以下面亲切的用Shai来称呼他。Shai 毕业于耶路撒冷希伯来大学，PHD论文便是online learning. 然后他去了online learning重地TTIC, 更是习得满身武艺，现在又回了耶路撒冷希伯来大学任教。他的数学直觉敏锐，善于将一些工具巧妙的应用到一些问题上。

两个概念

对一个凸函数，如下图所示，我们可以用两种不同的描述方式。一种就是函数上的所有点，函数曲线变由这些点构成；另外一种是每一个点所对应的切线，此时函数由所有切线形成的区域的上边缘确定。

函数上每一个点可由(w,f(w)) 表示，每一条切线可由斜率\theta 及与y轴交点，记为-f^\star(\theta) ，来确定。前一种表示是函数的primal form, 而后者则被称之为Fenchel conjudate, 是一种dual form, 其定义为

\displaystyle f^\star(\theta) = \max_w \langle w, \theta \rangle - f(w).

既然我们称它们为primal-dual, 当然意味着可以相互转换。实际上，对于一个闭凸函数，有f^{\star\star}=f .

另外一个概念是关于函数凸的程度的。一个函数被称为\sigma -强凸（strongly convex），那意味着对任意的w 和u , 有

\displaystyle f(w)-f(u)-\langle \nabla f(u), w-u\rangle\ge \frac{\sigma^2}{2}\|u-w\|^2.

示意图如下。直观上来说就是这个凸函数要凸得足够弯。这个定义可以拓展到一般的范数和次导数上来处理诸如稀疏罚\|\cdot\|_1 之类，此乃后话了。

Regret Bound分析

使用这两个概念可以得到一个重要定理。假设f 是\sigma -强凸的，且f^\star(0)=0 . 对任意的向量序列v_1,\ldots,v_T , 记w_t=\nabla f^\star\big(\sum_{j<t}v_j\big) 。那么对任意的向量u , 有

\displaystyle \sum_t \langle u, v_t \rangle - f(u) \le f^\star \big(\sum_t v_t\big) \le \sum_t \big( \langle w_t, v_t\rangle + \frac{1}{2\sigma}\|v_t\|^2\big).

此定理可以为Boosting, learning theory里的Rademecher Bound提供简单的证明，不过我们先关注它在online上的应用。不等式两边左右移项得\sum_t \langle u - w_t, v_t\rangle \le f(u) + \frac{1}{2\sigma} \sum_t \|v_t\|^2 . 注意到w_t 只和比t 小v_j 相关，所以我们可以取v_t=-\nabla \ell_t(w_t) ，立即由f 是凸函数有f(w_t)-f(u)\le \langle u - w_t, v_t \rangle . 同前记w^* 是离线最优解，且假设\|\nabla \ell_t(w_t)\|\le L , 联立上面两不等式有，R(T)\le f(w^*) + \frac{TL^2}{2\sigma} .

如果我们用\ell_2 罚，那么f=\lambda \|\cdot\|^2 , 且\sigma=2\lambda . 如果\lambda=\frac{L\sqrt{T}}{2\|w^*\|} ，则有R(T)\le LD\sqrt{T} , 这里D=2\|w^*\| 是对子空间\mathcal{H} 直径的一个估计。与上节Online Gradient Descent的regret bound比较，他们一样！

一般化的online算法

不小心先把理论讲完了，下面来解释下算法。接回一开始我们谈到的要考虑dual problem. 先观察到primal problem可以写成

\displaystyle \min_{w_0,\ldots,w_T}\left( f(w_0) + \sum_{t} \ell_t(w_t) \right)\quad \textrm{s.t.}\quad \forall t,\ w_t=w_0.

引入T 个拉格朗日乘子u_t , 那么拉格朗日形式就是

\displaystyle \mathcal{L}(w_0,\ldots,u_1,\ldots) = f(w_0) + \sum_t\ell_t(w_t) + \sum_t\langle u_t,w_0-w_t\rangle.

这样dual problem就是最大化函数\min_{w_t} \mathcal{L}(w_0,\ldots,u_1,\ldots) ，记为\mathcal{D}(\cdot) . 套进Fenchel conjugate的定义，就有

\displaystyle \mathcal{D}(u_1,\ldots,u_T)=-f^\star\big(-\textstyle \sum_t u_t\big)-\displaystyle \sum_t \ell^\star_t(u_t).

注意到与primal problem不同，这时罚里面的变量也被拆成了T 份，所以在t 时刻我们可以假设比t 大的u_j 都是0, 从而只更新u_t 来增大\mathcal{D} .

总结下算法。在t 时刻，我们先计算w_t = \nabla f^\star\big(-\sum_{i<t} u_i\big) ，在收到数据的label后，选择u_t 来保证足够大的增加dual:

\displaystyle \mathcal{D}(u_1,\ldots,u_{t-1},u_t,0,\ldots) \ge \mathcal{D}(u_1,\ldots,u_{t-1},\nabla \ell_t(w_t), 0, \ldots).

下面是一个示意图。在每个时间里我们根据新来的数据不断的增大dual, 同时使得primal变小。

一个值得讨论的是，在这里我们可以控制算法的aggressive的程度。例如，如果我们每次只取u_t=\ell_t(w_t) ，就是意味着维持最小但使得理论分析成立的步伐长度。我们也可以使用更加激进的步伐，u_t = \arg\max_u\mathcal{D}(\ldots,u_{t-1},u,0,\ldots) ，来最大限度的增加dual以期望收敛更快。不过这样通常使得每一步的计算开销更大。所以是仁者见仁智者见智了。

这个算法框架可以解释一系列算法，例如前面提到的Halving, 它使用激进步伐，或是WM, 使用平缓的步伐，以及Adaboost. 当然，我们可以从这个算法框架里推出新的有no-regret保证的新的online算法。

到这里，我们介绍的三个算法（前面是Online Gradient Descent, Weighted Majority）的regret bound都是\mathcal{O}(\sqrt{T}) 。事实上，对于强凸的目标函数，目前已知最优的regret bound是\mathcal{O}(\log T) 。记得在OGD中我们使用固定的\mathcal{O}(1/\sqrt{T}) 的学习率，如果我们将其改成变长的\mathcal{O}(1/t) ，那么变可以达到这个最优上界。对于这节介绍的方法，我们可以通常将罚f(w) 里面的\lambda 分成T 份给每个loss \ell_t , 这样前t 个loss的和就是t\lambda/T -强凸，从而达到使用变长学习率的效果。具体可见Shai的NIPS ’08.

漫谈在线学习：在线梯度下降

04/08/2011 in machine learning, online learning | 11 comments

在线凸优化

回顾下上次聊的专家问题，在t 时刻专家i 的损失是\ell_t(e^i) ，于是这个时刻Weighted Majority（WM）损失的期望是\sum_{i=1}^m w_t^i\ell_t(e^i) ，是关于这m 个专家的损失的一个线性组合（因为权重w_t^i 关于i 的和为1，所以实际上是在一个simplex上）。将专家在t 时刻的损失看成是这个时候进来的数据点，于是我们便在这里使用了一个线性的损失函数。

虽然WM的理论在上个世纪完成了[Littlestone 94, Freund 99]，将其理论拓展到一般的凸的函数还是在03年由Zinkevich完成的。当时Zinkevich还是CMU的学生，现在在Yahoo! Research。话说Yahoo! Research的learning相当强大，Alex Smola（kernel大佬），John Langford（有个非常著名的blog），这些年在large scale learning上工作很出色。

回到正题。我们提到在online learning中learner遭受的累计损失被记成\sum_{t=1}^T\ell_t(h_t) ，如果挑选h_t 的策略集\mathcal{H} 是凸的，而且损失函数\ell_t(h_t) 关于h_t 是凸的，那么我们称这个问题为Online Convex Optimization（OCP）。通常我们将h_t 表示成一个向量w_t ，例如WM中的维护的m 专家信任度，所以这时\mathcal{H}\subseteq\mathbb{R}^m 。

在线梯度下降

Zinkevich提出的算法很简单，在时刻t 做两步操作，首先利用当前得到数据对w_t 进行一次梯度下降得到w_{t+1} ，如果新的w_{t+1} 不在\mathcal{H} 中，那么将其投影进来：

\displaystyle w_{t+1}=\Pi_{\mathcal{H}}(w_t-\eta_t\nabla\ell_t(w_t)),

这里\nabla\ell_t(w_t) 是\ell_t(w_t) 关于w_t 的导数（如果导数不唯一，就用次导数），\eta_t 是学习率，\Pi_{\mathcal{H}}(w) 是投影子，其将不在\mathcal{H} 中的向量w 投影成一个与w 最近的但在\mathcal{H} 中的向量（如果w 已经在\mathcal{H} 中了，那就不做任何事），用公式表达就是\Pi_{\mathcal{H}}(w)=\arg\min_{u\in\mathcal{H}}\|w-u\| 。此算法通常被称之为 Online Gradient Descent。

先来啰嗦几句其与离线梯度下降的区别。下图是一个区别示意图。在离线的情况下，我们知道所有数据，所以能计算得到整个目标函数的梯度，从而能朝最优解迈出坚实的一步。而在online设定下，我们只根据当前的数据来计算一个梯度，其很可能与真实目标函数的梯度有一定的偏差。我们只是减少\ell_t(w_t) ，而对别的项是否也是减少就难说了。当然，我们一直在朝目标前进，只是可能要走点弯路。

那online的优势在哪里呢？其关键是每走一步只需要看一下当前的一个数据，所以代价很小。而offline的算法每走一个要看下所有数据来算一个真实梯度，所以代价很大。假定有100个数据，offline走10步就到最优，而online要100步才能到。但这样offline需要看1000个数据，而online只要看100个数据，所以还是online代价小。

于是，我们很容易想到，offline的时候能不能每一步只随机抽几个数据点来算一个梯度呢？这样每一步代价就会很少，而且可能总代价会和online一样的少。当然可以！这被称之为Stochastic Gradient Descent，非常高效，下回再谈把。

在这里，\mathcal{H} 的作用是什么呢？记得在learning中的目标函数通常是损失+罚（\ell(w)+\lambda \Omega(w) ）的形式。例如ridge regression就是平方误差+\ell_2 罚，lasso是平方误差+\ell_1 罚，SVM是hinge loss+\ell_2 罚。最小化这个目标函数可以等价于在\Omega(w)\le\delta 的限制下最小化\ell(w) 。\lambda 和\delta 是一一对应的关系。实际上\Omega(w)\le\delta 就是定义了一个凸子空间，例如使用\ell_2 罚，既\|w\|^2 ，时就是一个半径为\delta 的球。所以，Online Gradient Descent可以online的解这一类目标函数，只是对于不同的罚要选择不同的投影子。

Regret Bound分析

记投影前的\tilde w_{t+1} = w_t-\eta_t\nabla\ell_t(w_t) ，以及offline最优解w^*=\arg\min_{w\in\mathcal{H}}\sum_{t=1}^T\ell_t(w) 。因为\mathcal{H} 是凸的且w^* 在其中，所以对\tilde w_{t+1} 投影只会减少其与w^* 的距离，既\|w_{t+1}-w^*\|\le\|\tilde w_{t+1}-w^*\| 。简记\nabla_t=\nabla \ell_t(w_t) ，注意到

\displaystyle \|\tilde w_{t+1}-w^*\|^2=\|w_t-w^*\|^2+\eta_t^2\|\nabla_t\|^2-2\eta_t\langle\nabla_t,w_t-w^*\rangle.

由于\ell_t 是凸的，所以有

\displaystyle \ell_t(w_t)-\ell_t(w^*)\le \langle\nabla_t,w_t-w^*\rangle \le \frac{1}{2\eta_t}\big(\|w_t-w^*\|^2 - \|w_{t+1}-w^*\|^2\big) + \frac{\eta_t}{2}\|\nabla_t\|^2.

取固定的\eta_t=\eta ，对t 进行累加就有R(T)\le \frac{1}{2\eta}\|w_1-w^*\|^2+\frac{\eta}{2}\sum_{t=1}^T\|\nabla_t\|^2 。记\mathcal{H} 的直径为D ，且对所有t 有\|\nabla_t\|\le L 成立（既Lipschitz常数为L ），再取\eta=\frac{D}{L\sqrt{T}} ，那么

\displaystyle R(T)\le LD\sqrt{T}.

这个bound可以通过设置变动的学习率\eta_t 加强，具体下次再谈咯。

漫谈在线学习：专家问题

03/25/2011 in machine learning, online learning | 31 comments

两个故事

从前，网络的一头住着个挨踢男，另一头住着一小编。每天小编写一封垃圾邮件给挨踢男。苦命的挨踢男日日分析邮件设计过滤器来过滤垃圾邮件。但聪明的小编每天检查邮件有没有被挨踢男收到，不断增加更多的甜言蜜语来忽悠挨踢男。较量一直进行下去，挨踢男是否能摆脱小编的骚扰呢？

这个故事都属于博弈论里的重复游戏（repeated game），它是对在线学习（online learning）最贴切的刻画：数据不断前来，我们需根据当前所能得到的来调整自己的最优策略。

熟悉机器学习的稳拿可能注意到了在线学习与在机器学习（统计学习）里所讨论的（离线）学习的区别。前者认为数据的分布是可以任意的，甚至是为了破坏我们的策略而精心设计的，而后者则通常假定数据是服从独立同分布。这两种不同的假设带来不一样的设计理念和理论模型。

统计学习考虑算法所求得到的模型与真实模型的差距。数据由真实模型产生，如果能有无限数据、并在包含有真实模型的空间里求解，我们大概率能求得真实模型。但实际上我们只有有限数据，且只能考虑某类模型，因此很可能我们得不到真实解。所以，好的算法是能够用不多的数据来得到一个效果不错的模型。

在线学习的一个主要限制是当前只能看到当前的和过去的数据，未来是未知，有可能完全颠覆现在的认知。因此，对在线学习而言，它追求的是知道所有数据时所能设计的最优策略。同这个最优策略的差异称之为后悔（regret）：后悔没能一开始就选对最优策略。我们的希望是，时间一久，数据一多，这个差异就会变得很小。因为不对数据做任何假设，最优策略是否完美我们不关心（例如回答正确所有问题）。我们追求的是，没有后悔（no-regret）。

如果与统计学习一样对数据做独立同分布假设，那么在线学习里的最优策略就是在得知所有数据的离线学习所能得到最优解。因此在线学习可以看成是一种优化方法：随着对数据的不断访问来逐步逼近这个离线最优解。

很早以前优化界都在追求收敛很快的优化算法，例如牛顿迭代法。而最近learning的人却发现，好的优化算法很可能是坏的learning算法。这类算法虽然迭代一些次就能算出一个精度很高的解，但每一步都很贵，而且数据一大根本算不动。而向来被优化界抛弃的在线学习、随机优化算法（例如stochastic gradient descent），虽然收敛不快，不过迭代代价很低，更考虑到learning算法的解的精度要求不高，所以在实际应用中这些方法通常比传统的优化算法快很多，而且可以处理非常大的数据。

随便看看这些年ICML，NIPS，JMLR上online, stochastic关键字文章的数量就知道这类算法的大红大紫了。一个通常的规律是learning界火过的topic在computer vision，data mining等偏应用的方向上会接着火几年，所以肯定CVPR, ICCV, KDD之类会议上还会有茫茫多这类paper。因为最近一直在做这个topic，所以想乘自己知识还新，就自己的理解来谈谈它。希望能用通俗的语言为大家介绍这个topic，为国内learning界做做贡献。

当然，相对于老地主online learning，stochastic绝对是新贵。我会接下来谈这类算法，以及他们动辄数页的convergence rate的证明。我的另外一个topic是大规模矩阵的低秩近似，也算大规模机器里面的一个重要问题，我想有时间也来浅出一下。

下面是有公式的正文。

准确点的定义

我们把挨踢男都称之为learner，并从learner的角度来看问题。在t 时刻，learner收到问题x_t ，这可以是一封邮件，或者一个样本。learner然后从策略集\mathcal{H} 中选择一个策略h_t ，例如一个过滤器规则，或是一个线性分类器，并且对问题做出判断h_t(x_t) 。learner然后收到正确答案y_t 。

用损失函数\ell(h_t(x_t),y_t) 来衡量learner在t 时刻做出不正确判断时所遭受的损失，例如最常用的平方误差损失是(h_t(x_t)-y_t)^2 。为了方便起见，简记为\ell_t(h_t) 。那么在前T 轮learner遭受的总损失就是\sum_{t=1}^T\ell_t(h_t) 。记得我们说过用regret来衡量learner对一开始没能采用最优策略的后悔，记为R(T) ，有

\displaystyle R(T)=\sum_{t=1}^T\ell_t(h_t)-\min_{h\in\mathcal{H}}\sum_{t=1}^T\ell_t(h).

learner的目标就是每次挑不错的h_t 来使得R(T) 最小。一个自然的问题是，R(T) 可能小于0吗？好吧，是可以的。毕竟最优策略是要固定h ，而如果我们每次都能选很好的h_t 来适应问题(x_t,y_t) ，可能总损失会更小。但这非常难。因为我们是要定好策略h_t ，才能知道正确答案y_t 。如果这个问题和以前的很不一样，答案也匪夷所思，那么猜对它而且一直都猜对的概率很小。而对于最优策略来说，它能事先知道所有问题和答案，所以它选择的最优策略的总损失较小的可能性更大。所以，我们一般只是关心平均regretR(T)/T 是不是随着T 变大而变小。我们称一个online算法是不是no-regret，或者说online learnable的，意味着

\displaystyle \lim_{T\rightarrow\infty}\frac{R(T)}{T}\rightarrow 0.

这个定义不是很准确，因为问题x_t 是随机给的，所以这里我们要考虑最坏情况（小编写最狠的邮件）。为了符号的简单，就免去精确的数学定义了。（事实上，除非是统计人写的文章，优化或者机器学习的文章很多采用这种简单的写法）。

听专家的话

no-regret是不是很难达到呢？事实证明很多简单的算法都是no-regret。举个最经典的例子，假设我们有m 个专家，他们在每轮问题都会给出的答案。简单起见假设答案就两种，0和1。平方损失在这里就是0-1损失函数，答案正确损失为0，否则为1.

先考虑最简单的情况：假设至少有一个完美专家，她永远给出正确的答案。那么什么才是learner的好策略呢？最简单的很work：看多数专家怎么说罗。learner在时刻t 先挑出前t-1 轮一直是给正确答案的专家们，然后采用他们中多数人给出的那个答案。

此方法被称做Halving Algorithm。记C_t 是前t-1 轮一直是给正确答案的专家们的集合，|C_t| 是这些专家的个数。如果learner在时刻t 给出了错误的答案，那么意味着C_t 中大部分专家都给了错误答案，所以下一时刻learner参考的专家就会少一半，既|C_{t+1}|\le |C_t|/2 。因为完美专家的存在，所以|C_t| 不可能无限变小使得小于1，所以learner不能无限犯错。事实上，|C_t| 至多有\log_2(m) 次变小机会，所以learner最多有\log_2(m) 的损失。另一方面，最优策略当然是一直选中一位完美专家，总损失为0，所以regret的上界是个常数

\displaystyle R(T)\le \log_2(m).

现在考虑并没有传说中的完美专家存在的情况。这时维护C_t 的做法就不行了，我们使用一个稍复杂些的策略。记第i 个专家为e^i ，对其维护一个信任度w^i\in[0,1] ，且使得满足\sum_{i=1}^m w^i=1 。在t 时刻learner以w^i 的概率挑选专家e^i ，并用她的答案来做作为t 时刻的答案。

关键是在于如何调整信任度。直观上来说，对于在某一轮预测不正确的专家，我们需降低对她们的信任度，而对预测正确的，我们则可以更加的相信她们。一种常见的调整策略是基于指数的。我们先维护一个没有归一化（和不为1）的信任度u 。一开始大家都为1，u_0^i=1 . 在t 时刻的一开始，我们根据i 专家在上一轮的损失，记为\ell_{t-1}(e^i) ，来做如下调整：

\displaystyle u_{t}^i=u_{t-1}^i\exp(-\eta\ell_{t-1}(e^i)),

这里\eta 是学习率，越大则每次的调整幅度越大。然后再将u_t 归一化来得到我们要的分布：w_t^i=u_t^i/\sum_i u_t^i 。此算法被称之为Exponential Weights或者Weighted Majority。

基于指数的调整的好处在于快速的收敛率（它是强凸的），以及分析的简单性。我们有如下结论：如果选择学习率\eta=\sqrt{8\ln m/T} ，那么

\displaystyle R(T)\le \sqrt{\frac{T\ln m}{2}}.

证明先分析\sum_i u^i_{t}/\sum_i u^i_{t-1} ，其上界是关于最优策略的损失，下界是关于learner的损失，移项就得regret的上界（online/stochastic算法的收敛性分析大都走这种形似）。具体的证明这里略过，后面后补上文献。

一个值得讨论的问题是，设定最优学习率\eta 需要知道数据的总个数T ，而在实际的应用中一般不能知道样本到底会有多少。所以如果设定work的参数很trick。在以后的算法中还会遇到这类需要上帝之手设定的参数，而如何使得算法对这类控制速率的参数不敏感，或者有意的调整参数来控制算法的灵敏度，都是当前研究的热点。这是后话了。

另外一个值得注意的是，同前面存在完美专家的情况相比，regret的上界由\log_2 m 增到了\sqrt{T\ln m/2} 。这两者在实际的应用中有着很大差别。例如我们的目标是要使得平均regret达到某个数值以下，假设前一种方法取1,000个样本（迭代1,000次）就能到了，那么后一种算法就可能需要1,000,000个样本和迭代。这样，在时间或样本的要求上，前者明显优于后者。类似的区别在后面还会更多的遇到，这类算法的一个主要研究热点就是如何降低regret，提高收敛速度。

这次谈的都是online的经典算法。下回我们介绍convex online optimization，这与机器学习更相关。