其实感觉 Derivate-free optimization(DFO)的方法当中，大多数还是需要依靠导数的，尤其是在确定算法的初值的时候。嘛，不过相比于其他的方法来说，倒确实少了很多。

DFO的方法往往用于处理无约束优化或者带有简单约束的优化(如边界等)。常见的做法往往是进行插值拟合，然后对拟合出的函数运用已有的优化算法。

model-basd method

根据 Taylor theorm，函数可以展开为下面的式子：
$m_k(x_k+p)=c+g^Tp+\frac{1}{2}p^TGp$
我们可以拟合上面的函数得到原函数的近似。插值点的个数 $q=\frac{1}{2}(n+1)(n+2)$ ，因为需要估算出这么多的变量(假设 $G$ 是对称的)。这种情况下，我们需要解 $q$ 个等式：
$$m_k(y^l) = f(y^l)\quad l=1,2,3...q \label{eq:sim}$$
其实本质上来说，这就是一个线性方程组。为了保证恒有解，我们挑选的 $y^l$ 必须保证系数非奇异。

确定了 $m_k$ 之后，就可以用其他方法去算更新量了。书中给的是 Trust-region，所以这里我们也以 Trust-region 为例进行说明。用 Trust-region 去计算更新量，也就是计算
$\min_p m_k(x_k+p) \quad s.t. \|p\|_ 2\le\Delta$
和 Trust-region 一样，我们以比值
$\rho=\dfrac{f(x_k)-f(x_k^+)}{m_k(x_k)-m_k(x_k^+)}$
表征本次迭代的效果($x_k^+$ 为下一次可能的$x_k$)。
但是和 Trust-region 不同的是，我们需要对 $m_k$ 做更新。

算法

先描述下整个算法的流程。对于 $m_k$ 更新的部分，我们在后面再说。算法的流程如下：

select $Y = \{y^1,y^2...y^q\}$ satisfy condition
select $x_0$ satisfy $f(x_0)\le f(y^i)$ for all $y^i\in Y$
given $\Delta_0,\ \eta\in(0,1)$
repeat
- form $m_k(x_k+p)$
- get $p$ by trust-region
- $x_k^+=x_k+p$
- compute $\rho$
- if $\rho > \eta$
  - replace an element in $Y$ by x
  - choose $\Delta_{k+1}\ge\Delta_{k}$
  - $x_{k+1}=x_k^+$
  - continue
- else if $Y$ don’t need improve:
  - choose $\Delta_{k+1}<\Delta_{k}$
  - $x_{k+1}=x_k$
  - continue
- update $Y$
- $\Delta_{k+1}=\Delta_{k}$
- $x=\arg_y\min_{y\in Y}f(y)$
- $x_k^+ = x$, compute $\rho$
- if $\rho \ge\eta$, $x_{k+1}=x_{k}^+$
- else $x_{k+1}=x_{k}$
  
  until convergence

算法首先会判断当前迭代的效果。和 Trust-region 一样，效果好的话，扩大 $\Delta$ 。效果不好则存在2种情况：用 $Y$ 解出的拟合效果不好；或者是 $\Delta$ 过大。 $Y$ 的拟合效果不好可以理解为 \ref{eq:sim} 中的解为一个低维子空间。如果 \ref{eq:sim} 中的系数矩阵 condition number 过大，就会出现这种问题(我也不知道为什么)。在这种情况下，对 $Y$ 进行调整即可；否则就调整 $\Delta$ 。

对于 $Y$ 的初值来说，选择单纯形每条边的中点即可。

可以看到每次迭代我们需要对 $m_k$ 进行拟合。对于二次函数来说，拟合的计算代价过大。我们可以将二次函数切换为线性函数，减小计算量。

插值模拟

对于插值拟合来说，我们需要解决
$m_k(y^l) = f(y^l)$
而
$m_k(x_k+p) = f(x_k) + g^Tp + \sum_{i<j}G_{ij}p_ip_j + \frac{1}{2}\sum_{i}G_{ii}p_i^2$
令
$$\begin{align} \hat{g} &= (g^T,\{G_{ij}\}_{i<j}，\{\frac{1}{\sqrt{2}}G_{ii}\})^T\\ \hat{p} &= (p^T,\{p_ip_j\}_{i<j}，\{\frac{1}{\sqrt{2}}p_{i}^2\})^T \end{align}$$
则 $m_k(x_k) = f(x_k)+\hat{g}^T\hat{p}$。直接求解等式即可。当然，我们也可以假设
$m_k(x) = \sum_{i=1}^q\alpha_i\phi_i(x)$
可以求解出 $\alpha$ 当
$\delta(Y) = det\begin{pmatrix}\phi_1(y^1)&\cdot\cdot\cdot & \phi_1(y^q)\\\cdot& & \cdot\\\cdot& & \cdot\\\cdot& & \cdot\\\phi_q(y^1)&\cdot\cdot\cdot & \phi_q(y^q)\end{pmatrix}$
非奇异。但是随着迭代次数的增加， $\delta(Y)$ 会越来越小，最后趋向于0。我们需要通过 $Y$ 的更新来解决这个问题。

$Y$ 的更新

为了解决上面的问题，我们需要保证在每轮迭代时 $\delta(Y)$ 不会下降。定义 $L(\cdot,y)$ 对于 $Y$ 中的所有点 $y, \bar{y}$ 满足如下条件：
$L(\bar{y}, y) = \begin{cases}0 &\bar{y}=y\\1&\bar{y}\neq y\end{cases}$
则有
$\|\delta(Y^+)\| \le \|L(y_+,y_-)\|\|\delta(Y)\|$
其中 $y_-$ 是要被替换的点， $y_+$ 是替换 $y_-$ 的点，$Y$ 是原来的集合， $Y^+$ 是替换后的集合。

再回到上面的算法。在第一个条件满足($\rho\ge\eta$)时，我们直接让 $y_+ = x^+$ ，然后去求 $y_-$：
$y_- = \arg\min_{y\in Y}|L(x^+,y)|$

该条件不满足的情况下，首先需要判断是不是 $Y$ 的问题。如果对于所有的 $y^i\in Y$ 的情况下都有 $\|x_k-y^i\|\le\Delta$ 。这样的情况下我们需要更新 $\Delta$ 。否则我们需要找到一个可以增加 $\delta(Y)$ 的 $y_+$ 。对于每一个 $y^i\in Y$ ，都可以找到对应的 $y^+_ i$ ：
$y^+_ i = \arg\max_{\|y-x_k\|\le\Delta} |L(y,y^i)|$
我们需要找到 $\max|L(y^+_ i, y^i)|$ 的 $y^i$ 作为 $y_-$。

基于 mininum-change 的更新

书中简单的提到了这种方法。在这种方法中，我们最小只需要维持 $n+1$ 个点，每轮迭代的优化目标为：
$$\min_{f,g,G} \|G-G_k\|_ F\\ s.t.\quad G^T=G\\ m(y^l) = f(y^l)\quad l=1,2,...\hat{q}$$
其中 $\|\cdot\|_ F$ 为 Frobenius 范式，即矩阵每个元素的平方和。为了保证 $G$ 和 $G_k$ 不一样，$\hat{q}$ 至少得大于 $n+1$；推荐为 $2n+1$ ，同时还需要保证 $y^i$ 不共面。

不过这种方法应该是最近提出的，所以介绍不多。

coordinate

坐标梯度下降，一种很简单的方法，每次下降的方向都是一个基坐标向量。在遍历完所有的下降方向之后，再从头开始。坐标梯度下降的主要问题在于收敛的判断条件。在这种情况下我们不可能用 $\|\nabla f\| <\epsilon$ 所以很难判断什么时候收敛。

Pattern-search

Pattern-search 类似于一种泛化的坐标梯度下降，每次迭代时从一个备选集中选出当前的迭代方向，在这个方向上下降。我们将这个备选集定义为 $D_k$ ，算法的具体过程如下：

given $\gamma_{tol},\ \theta_{max}$ , decrease function $\rho(t)$ that $\lim_{t\to 0} \dfrac{\rho(t)}{t}=0$
choose $x_0,\ \gamma_0>\gamma_{tol},\ D_0$
repeat
- if $\gamma_k<\gamma_{tol}$ break
- if $f(x_k+\gamma_kp_k)<f(x_k) - \rho(\gamma_k)$ for some $p_k\in D_k$
  - $x_{k+1} = x_k+\gamma_kp_k$
  - $\gamma_{k+1} = \phi_k\gamma_k$ for $\phi_k>1$
- else
  - $x_{k+1}=x_k$
  - $\gamma_{k+1} = \theta_k\gamma_{k},\ 0<\theta_k\le \theta_{max}$

$D$的选择

$D$ 在这个算法中的重要性自然不用说。所以 $D_k$ 中方向的选择就成了一个重要的问题。我们需要建立一些基础条件。在 Line search 中，每次迭代时都有 $\cos\theta=\dfrac{-\nabla f_k^Tp}{\|\nabla f_k\|\|p\|}\ge\delta$ 。所以对于 $D_k$ 中的任何一个方向 $p$ 来说，都要满足上面的条件。我们将上面的条件描述为：
$k(D_k) = \min_{v\in R^n}\max_{p\in D_k} \dfrac{v^Tp}{\|v\|\|p\|}\ge \delta$

我们还需要让 $D_k$ 中的每个向量的长度保持在一定范围内。这样的话，不同向量的选择不会导致步长的变化过大。也就是说
$\beta_{min} \le \|p\|\le\beta_{max}$
结合上面两个条件，可以得到
$-\nabla f_k^T p \ge k(D_k)\|\nabla f_k\|\|p\|\ge\delta\beta_{min}\|\nabla f_k\|$
有一些满足简单的满足条件的 $D$ 的例子：

$\{e_1,e_2,...e_n,-e_1,-e_2,...-e_n\}$
令 $e = (1,1,...,1)^T$，则 $p_i=\frac{1}{2n}e-e_i, i=1,2,..,n;\ p_{n+1} = \frac{1}{2n}e$也是一组满足条件集合

当然了，上面的例子往往只是 $D_k$ 的一个子集。其他的方向可以通过启发式算法得到，但是书中没有详细说明，所以这里也没法介绍(+ _ +)

$\rho(t)$

在上面的算法中我们引入了一个函数 $\rho(t)$ 书中给出的建议是 $\rho(t) = Mt^{\frac{3}{2}}$，$M$ 为正数

conjugate-direction method

这里的方法和 conjugate gradient 是一样的，只是这里在求向量的时候没有用到梯度(或者说很少用？)
和 conjugate gradient 一样，我们还是从二次函数说起。假设优化目标为
$f(x) = \frac{1}{2}x^TAx-b^Tx$
假设我们的优化方向为 $l_1(\alpha) = x_1+\alpha p,\ l_2(\alpha) = x_2+\alpha p$，且沿着两个方向优化的最终结果为 $x_1^*,\ x_2^*$ ，则 $x_1^*- x_2^*$ 与 $p$ 关于A 共轭。

将上面的结果推而广之，对于 $l$ 个方向的集合 $\{p_1,p_2,...,p_l\}$ 来说，
$$S_1=\{x_1+\sum^l_{i=1}\alpha_ip_i|\alpha_i\in R, i=1,2,...,l\}\\ S_2=\{x_2+\sum^l_{i=1}\alpha_ip_i|\alpha_i\in R, i=1,2,...,l\}$$
而 $x_1^*,\ x_2^*$ 为沿着两个方向计算出的最小值点，则 $x_1^*- x_2^*$ 与 $p$ 关于A共轭。

在这个结论的基础上，我们就可以构建算法了：

choose $x_0,\ p_i=e_i, i=1,2,..,n$
compute $x_1$ to min $f(x_0+\alpha p_n)$
repeat
- $z_1 = x_k$
- compute $\alpha_j = \arg\min f(z_j+\alpha_jp_j)$ and $z_{j+1} = z_j+\alpha_jp_j$ for j = 1,2,…,n
- $p_j = p_{j+1}$ for j = 1,2,…,n-1
- $p_n = z_{n+1}-z_1$
- compute $\alpha_n = \arg\min f(z_{n+1}+\alpha_np_n)$
- $x_{k+1} = z_{n+1}+\alpha_np_n$
  until convergence

简单解释一下上面的算法。假设 $n=3$ ，即我们一开始有3个向量 $e_1,e_2,e_3$。按照算法的步骤，我们计算出了 $x_1$ 以及 $z_4$ 。根据上面的结论 $p_1 = z_4-x_1$ 和 $e_1,e_2,e_3$ 共轭，而 $x_2$ 为 $f$ 在 $S_1=\{y+\alpha_1e_3+\alpha_2p_1\}$ 方向上的最小值。这种情况下我们得到了新的集合 $\{e_2,e_3,p_1\}$

非线性系统

上面的算法可以移植到任何的函数上，但是需要一些修改。我们不修改方向，只对长度做一些调整：
$\hat{p}_i = \dfrac{p_i}{\sqrt{p_i^TAp_i}}$
则如果我们最大化
$|det(\hat{p}_1,\hat{p}_2,...,\hat{p}_n)|$
那么 $p$ 就是关于A共轭的。我们可以调整上面的算法：

find $m\in \{1,2,...,n\}$ max $\psi_m=f(x_{m-1}) - f(x_m)$
$f_1=f(z_1), f_2=f(z_2),f_3=f(z_3)$
if $f_3\ge f_1$ or $(f_1-2f_2+f_3)(f_1-2f_2-\psi_m)\ge \frac{1}{2}\psi_m(f_1-f_2)^2$
- set $x_{k+1} = z_{n+1}$
else
- set $\hat{p}=z_{n+1}-z_1$, compute $\hat{\alpha}$ to min $f(z_{n+1}+\hat{\alpha}\hat{p})$
- $x_{k+1} = z_{n+1} + \alpha\hat{p}$
- remove $p_m$ , add $\hat{p}$

Nelder-mead method

The Nelder-mead simplex-reflection method，顾名思义和单纯形有关。在每次迭代时会所有的点的 convex hull 是一个单纯形。假设我们需要优化一个 $n$ 维函数 $f(x)$ 。易知， $n$ 维单纯形有 $n+1$ 个点。我们假设这 $n+1$ 个点为 $\{x_1,x_2,...,x_{n+1}\}$ ,其顺序满足下面的条件：
$f(x_1)\le f(x_2)\le...f(x_{n+1})$
令 $\bar{x} = \sum_{i=1}^nx_i, \bar{x}(t) = \bar{x}+t(x_{n+1}-\bar{x})$ 。我们的算法就是用 $\bar{x}(t)$ 替代最差的 $x_{n+1}$。算法的一次迭代如下：

compute $\bar{x}(-1), f_{-1}=f(\bar{x}(-1))$
if $f(x_1)\le f_{-1}\le f(x_n)$
- $x_{n+1}=\bar{x}(-1)$ and continue
else if $f_{-1}< f(x_1)$
- compute $f_{-2}$
- if $f_{-2} < f_{-1}$
  - $x_{n+1}=\bar{x}(-2)$ continue
  - $x_{n+1}=\bar{x}(-1)$ continue
else if $f_{-1}> f(x_n)$
- if $f(x_n)\le f_{-1} \le f(x_{n+1})$
  - compute $f_{-1/2}$
  - if $f_{-1/2} \le f_{-1}$
    - $x_{n+1}=x_{-1/2}$ continue
- else
  - compute $f_{1/2}$
  - if $f_{1/2} < f_{n+1}$
    - $x_{n+1}=x_{1/2}$ continue
$x_i = (1/2)(x_1+x_i), i=2,3,...,n+1$

在算法陷入停滞时，可以进行 restar。不过对于初值点的选择方面，似乎没有什么特别优秀的方案。

Implicit Filtering

一个非常奇怪的方法= =。算法如下：

choose $\epsilon_k\to 0$, $c,\rho \in (0,1), \alpha_{max}$
repeat
- increment_k = false
- repeat
  - comupte $f(x)，\nabla_{\epsilon_k}f(x)$
  - if $\|\nabla_{\epsilon_k}f(x)\|\le \epsilon_k$
    - increment_k = true
  - else
    - find $m \in (0,a_{max})$ that $f(x-\rho^m\nabla_{\epsilon_k}f(x))\le f(x)-\rho^m\|\nabla_{\epsilon_k}f(x)\|^2$
    - if $m$ not exist
      - increment_k = true
    - else
      - $x = x -\rho^m\nabla_{\epsilon_k}f(x)$
  until increment_k
- $x_k = x$
until convergence

当然，也可以用 $\nabla_{\epsilon_k}f$ 估计二阶导，然后用牛顿法解。