Supervised Descent Method 简单实现

前些天看了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.

以下是原文中四个一维函数例子的实现。

[plain]  view plain  copy

1. close all;  
2.   
3. testCase = 3;  
4.   
5. switch testCase  
6.     case 1  
7.         steps = [0.25, 0.025, 0.0025, 0.00025];  
8.     case 2  
9.         steps = [3, 1, 0.5];  
10.     case 3  
11.         steps = [3, 1, 0.5];  
12.     case 4  
13.         steps = [0.11, 0.055, 0.011];  
14. end  
15.   
16. terror = zeros(length(steps), 1);  
17. c = 0;  
18. for sid = 1:length(steps)  
19.     step = steps(sid);  
20.       
21.     switch testCase  
22.         case 1  
23.             h = @(x) sin(x);  
24.             g = @(x) asin(x);  
25.             y0 = [-1:step:1]';  
26.             x0 = g(y0);  
27.             ystar = [-1:0.05:1]';  
28.         case 2  
29.             h = @(x) exp(x);  
30.             g = @(x) log(x);      
31.             y0 = [1:step:28]';  
32.             x0 = g(y0);  
33.             ystar = [1:0.5:28]';  
34.         case 3  
35.             h = @(x) x.^3;  
36.             g = @(x) nthroot(x, 3);  
37.             y0 = [-27:step:27]';  
38.             x0 = g(y0);  
39.             ystar = [-27:0.5:27]';  
40.         case 4  
41.             h = @(x) erf(x);  
42.             g = @(x) erfinv(x);  
43.             y0 = [-0.99:step:0.99]';  
44.             x0 = g(y0);  
45.             ystar = [-0.99:0.03:0.99]';  
46.     end  
47.       
48.     n = 10;  
49.       
50.     % train model  
51.     R = cell(n,1);      x = cell(n, 1);  
52.     phi = cell(n, 1);   error = zeros(n, 1);  
53.     x{1} = repmat(c, length(x0), 1);  
54.     error(1) = norm(x{1} - x0);  
55.     for k=2:n  
56.         dx = x0 - x{k-1};  
57.         phi{k-1} = h(x{k-1});  
58.         dphi = y0 - phi{k-1};  
59.         R{k-1} = (dphi'*dphi)\(dphi'*dx);  
60.         b{k-1} = mean(dx) - R{k-1} * mean(dphi);  
61.         x{k} = x{k-1} + R{k-1} * dphi;  
62.         error(k) = norm(x{k} - x0);  
63.     end  
64.       
65.     figure; hold on;  
66.     plot(x0, y0, '-gx');  
67.     plot(x{n}, y0, '-ro');  
68.     title('training');  
69.       
70.     figure; plot(error);  
71.     title('error');  
72.       
73.     % test it  
74.     xstar = cell(n, 1);  
75.     xstar{1} = repmat(c, length(ystar), 1);  
76.     for k=2:n  
77.         xstar{k} = xstar{k-1} + R{k-1} * (ystar - h(xstar{k-1})) + b{k-1};  
78.     end  
79.     diffx = xstar{n} - g(ystar);  
80.     diffx = diffx(2:end-1);  
81.     terror(sid) = sqrt((diffx'*diffx)/length(xstar))  
82.     figure; hold on;  
83.     plot(g(ystar), ystar, '-gx');  
84.     plot(xstar{n}, ystar, '-ro');  
85.     title('test');  
86. end  
87.   
88. figure; plot(log(steps), terror, '-x');