前些天看了Supervised Descent Method and Its Applications to Face Alignment这篇文章，非常喜欢。这篇文章提出了一种基于机器学习来解决复杂最小二乘问题(least squares problem)的方法(简称SDM方法)。该方法思路很简洁，从训练数据中学习梯度下降的方向并建立相应的回归模型，然后利用得到的模型来进行梯度方向估计。
牛顿法是迭代求解最小二乘问题的常用方法，但是由于直接求Hessian矩阵往往代价非常高，人们常常采用别的方法来近似Hessian矩阵。这些采用近似Hessian矩阵的方法成为拟牛顿法。只要估计的Hessian矩阵足够接近真实的Hessian矩阵，拟牛顿法就可以高效地解决最小二乘问题。即便如此，从最小二乘的代价函数估计梯度、Hessian矩阵对于某些复杂的最小二乘问题仍然是非常困难的(或者根本就不可能)，比如这篇文章中提到的关于SIFT特征的梯度以及Hessian矩阵。
那么能不能从别的角度来估计Hessian矩阵和梯度向量，或者，Hessian矩阵的逆与梯度向量的乘积呢？(其实对于牛顿法而言，梯度向量与Hessian逆矩阵的乘积直接决定了优化方向，所以如果可以直接估计他们的乘积，效果上是一样的。)
可是怎么估计呢？从大量的数据里面构造出这个关键的矩阵/向量来就好了嘛！假设要求解的最小二乘问题是 || f(x) – y ||^2 ，其中f是一个非常复杂的函数，它的计算方法已知，但是梯度和Hessian矩阵均无法直接求解。只要给定一个足够大的最优解集合 {(x*, y*) }， 就可以通过模拟牛顿法的过程来逐步构造出Hessian矩阵和梯度向量。不过直接构造Hessian矩阵和梯度向量需要比较多的训练数据，因为自由度往往很大。既然牛顿法里面最后只需要Hessian的逆和梯度向量的乘积(，估计这个向量显然要直接得多，也容易一些。SDM方法假设这个乘积是x的一个复杂函数，并且假设这个函数可以被一组函数的线性组合近似。具体而言，SDM用以下迭代公式
x_{k+1} = x_k + R_k * f(x_k) + b_k
取代了牛顿法的迭代公式
x_{k+1} = x_k – H(x)^{-1} f'(x)
其中R_k 和 b_k 是学习得到的用于近似Hessian矩阵和梯度的模型参数。通过模拟牛顿法的过程，可以通过不断得修正一组初始估计来逼近最优解集合，每一步迭代都可以计算出一对R_k 和 b_k，它们一起构成了最后的模型。
求解 R_k 和 b_k 的方法类似于贪心算法，即在每一步，，最小化迭代后的值与最优值的差:
|| x* – x_k + R_k * f(x_k) + b_k ||^2
这个其实就是最简单的线性最小二乘问题，可以很容易的求解。
当然啦，实际训练的时候，需要在所有训练样例计算上面的代价函数并求和，得到最终的代价函数，然后再求解R_k和b_k.
以下是原文中四个一维函数例子的实现。
close all;testCase = 3;switch testCasecase 1steps = [0.25, 0.025, 0.0025, 0.00025];case 2steps = [3, 1, 0.5];case 3steps = [3, 1, 0.5];case 4steps = [0.11, 0.055, 0.011];endterror = zeros(length(steps), 1);c = 0;for sid = 1:length(steps)step = steps(sid);switch testCasecase 1h = @(x) sin(x);g = @(x) asin(x);y0 = [-1:step:1]';x0 = g(y0);ystar = [-1:0.05:1]';case 2h = @(x) exp(x);g = @(x) log(x);y0 = [1:step:28]';x0 = g(y0);ystar = [1:0.5:28]';case 3h = @(x) x.^3;g = @(x) nthroot(x, 3);y0 = [-27:step:27]';x0 = g(y0);ystar = [-27:0.5:27]';case 4h = @(x) erf(x);g = @(x) erfinv(x);y0 = [-0.99:step:0.99]';x0 = g(y0);ystar = [-0.99:0.03:0.99]';endn = 10;% train modelR = cell(n,1);x = cell(n, 1);phi = cell(n, 1); error = zeros(n, 1);x{1} = repmat(c, length(x0), 1);error(1) = norm(x{1} – x0);for k=2:ndx = x0 – x{k-1};phi{k-1} = h(x{k-1});dphi = y0 – phi{k-1};R{k-1} = (dphi'*dphi)\(dphi'*dx);b{k-1} = mean(dx) – R{k-1} * mean(dphi);x{k} = x{k-1} + R{k-1} * dphi;error(k) = norm(x{k} – x0);endfigure; hold on;plot(x0, y0, '-gx');plot(x{n}, y0, '-ro');title('training');figure; plot(error);title('error');% test itxstar = cell(n, 1);xstar{1} = repmat(c, length(ystar), 1);for k=2:nxstar{k} = xstar{k-1} + R{k-1} * (ystar – h(xstar{k-1})) + b{k-1};enddiffx = xstar{n} – g(ystar);diffx = diffx(2:end-1);terror(sid) = sqrt((diffx'*diffx)/length(xstar))figure; hold on;plot(g(ystar), ystar, '-gx');plot(xstar{n}, ystar, '-ro');title('test');endfigure; plot(log(steps), terror, '-x');
在人生的大海中，我们虽然不能把握风的大小，却可以调整帆的方向。
