流形优化:从一个简单的例子开始

问题引入

先考虑一个经典的优化问题:瑞利商(Rayleigh quotient),如下:
\min_{x \in \real^n} \frac{x^TAx}{x^Tx}
由于瑞利商不随尺度变化而变化,因此,可以将x限制为单位范数:||x||=1,包含这样的向量的集合称为单位球:
S^{n-1}={x \in \real^n:||x||=1}
那么原问题可以转换为如下形式:
\min_{x\in S^{n-1}}x^TAx
或者是:
\min_{x \in \real^n} {x^TAx}\\ s.t. \space ||x||=1
这里选择一个具体的A作为示例
A = \begin{vmatrix} -1 & 0\\ 0 & 0 \end{vmatrix}
将(5)带入(4)得:
\min _{x_i\in\real}f(x)=-x_1^2\\ s.t. \space x_1^2+x_2^2=1


拉格朗乘数法

这个优化问题并不复杂,可以使用拉格朗乘数法求解。首先构造拉格朗日函数:
L(x_1,x_2,\lambda)=-x_1^2+\lambda(x_1^2+x_2^2-1)
L(x_1,x_2,\lambda)对三个参数的一阶偏导数等于零,即
image-20200607160851689
解上述方程可以得到目标函数的极值f_{max}=f(0,1)=0,f_{min}=f(1,0)=-1


流形优化

如上文所述,x的集合是一个单位球S^{n-1},这是一个很典型的流形结构,我们可以尝试使用流形的方法求解该问题。

类似于欧式空间上的最速下降法,在流形上,也可以求得流形上的梯度,从而进行最速下降。不同的是,为了保证每一步x都在流形上,需要进行拉回(retraction)操作,算法流程如下:
image-20200607160447664
使用几何图形对上述过程近似描述,如下图:

image-20200607141233294

回到上述问题,先计算目标函数在欧式空间的梯度
{\rm grad}\bar f(x)=\nabla f(x)=\begin{bmatrix}-2x_1\\0\end{bmatrix}
接下来,我们需要把改欧式空间梯度正交投影到切空间T_xM中,从而得到切空间内的梯度。欧式空间上的向量u向单位球的切空间的正交投影为
{\rm Proj}_x(u)=u-(x^Tu)x=(I_n-xx^T)u
将刚刚求得的目标函数在欧式空间中的梯度带入
{\rm grad}f(x)={\rm Proj}_x({\rm grad}\bar f(x))=\begin{bmatrix} -2x_1(1-x_1^2)\\ 2x_1x_2^2\end{bmatrix}
得到切空间中的梯度后,需要将其拉回到球面,球面的拉回函数如下:
R_x(v)=\frac{x+v}{||x+v||}
先在,公式和算法流程已经介绍完,我们尝试手动迭代解该问题。令初始位置x_0=[0.5, 0]^T,迭代步长\alpha_k \equiv 1

iteration 1
image-20200607160545952
iteration 2
image-20200607160600748
迭代到第二步时,便达到了终止条件,从而得到最优解f_{min}=f(1,0)=-1。该结果与使用拉格朗日乘数法得出的结果是一致的。


程序实现

使用编程语言实现流形优化很复杂,好在已经有现成的流形优化工具包供我们使用。本文使用的是manopt工具包。

安装

你可以在三中编程语言中使用manopt,分别是Matlab,Python和Julia,本文介绍前两种。

  • 【Matlab】下载manopt,在Matlab中进入manopt文件夹,执行importmanopt即可。

  • 【Python】执行下述指令即可。

    pip3 install --user pymanopt
    

使用manopt解决上述问题

在Matlab中,使用如下代码:

% (1) Instantiate a manifold
manifold = spherefactory(2);

% (2) Define the cost function and its Euclidean gradtitude
A = [1 0;0 0];
problem.cost = @(x) -x'*(A*x);
problem.egrad = @(x) -(2*A*x);
problem.M = manifold;

% (3) Choose a manopt solver to optimize
x0 = [0.5, 0]';
[x, xcost, info, options] = steepestdescent(problem,x0)

从info中查看迭代信息,发现和我们使用手工迭代的步骤是一致的。

image-20200607152606149

在Python中,使用如下代码:

import numpy as np

from pymanopt.manifolds import Sphere
from pymanopt import Problem
from pymanopt.solvers import SteepestDescent

# (1) Instantiate a manifold
manifold = Sphere(2)

A = np.array([[1,0],[0,0]])

# (2) Define the cost function
def cost(X): 
    xx = np.array(X)
    return  -np.dot(np.dot(xx.transpose(),A),xx)

problem = Problem(manifold=manifold, cost=cost)

# (3) Instantiate a Pymanopt solver
solver = SteepestDescent()

# let Pymanopt do the rest
Xopt = solver.solve(problem)
print(Xopt)

由于Pymanopt中设置起始点很麻烦,这里使用的随机起始点,得到的结果如下

image-20200607151713907

经过了19次迭代,最终逼近了最小值-1。