第43卷    计  算  机  学  报  Vol. 43 2020  论文在线号No.6  CHINESE JOURNAL OF COMPUTERS   2020Online No.6 ——————————————— 本课题得到国家自然科学基金  (61773410, 61673403)和中国博士后科学基金(2019M663234)资助.何笑雨，博士，主要研究领域为演化算法，hxyokokok@foxmail.com。周育人（通信作者），博士，主要研究领域为演化算法，zhouyuren@mail.sysu.edu.cn。陈泽丰，博士，主要研究领域为演化算法，zefeng.chen@ntu.edu.sg。  基于演化策略的矩阵流形黑盒优化方法 何笑雨1)  周育人1)    陈泽丰2) 1)中山大学数据科学与计算机学院  广州510006   2)南洋理工大学计算机科学与工程学院  新加坡  639798  摘  要  该文使用演化算法求解一类定义在黎曼流形之上的矩阵优化问题。传统演化算法大多针对欧式空间中的优化问题而设计，因此难以直接用于矩阵流形优化。该文将经典的协方差矩阵适应技术从欧式空间扩展到黎曼流形，提出了用于求解矩阵流形黑盒优化问题的流形搜索方向适应演化策略。该算法利用流形的几何结构，将原始的带约束问题转换为一系列切空间中的无约束问题。在每个切空间中，使用一个单位阵和若干个搜索方向构造协方差矩阵，并使用相应的多元高斯分布引导搜索。所得算法不依赖于全局坐标系、独立于特定坐标基的选择，所有演化算子均通过矩阵变换进行，达到较高的运行效率。在三类矩阵流形优化问题上进行的数值实验表明，所提出算法的性能显著优于或相当于主流的流形优化器。 关键词  演化策略；矩阵流形；搜索方向适应；黑盒优化；协方差矩阵适应 中图法分类号  TP18 An Evolution Strategy for Black-box Optimization on Matrix Manifold He Xiao-Yu1)  Zhou Yu-Ren1)    Chen Ze-Feng2) 1) School of Data and Computer Science, Sun Yat-sen University, Guangzhou 510006 2) School of Computer Science and Engineering, Nanyang Technological University, Singapore 639798  Abstract  This work concerns the black-box approach towards a class of matrix optimization problems defined on  Riemannian  manifolds.  Due  to  their  non-Euclidean  nature,  these  problems  are  intractable  for  traditional evolutionary  algorithms.  This  work  extends  the  covariance  matrix  adaptation  from  Euclidean  space  to  matrix manifolds, and proposes a manifold search direction adaptation evolution strategy (MSDA-ES). By exploiting the manifold structure, we transform an originally constrained problem into a sequence of unconstrained ones in the tangent spaces. In each tangent space, we use a multi-variate Gaussian distribution to guide the search where the  covariance  matrix  is  modeled  with  an  identity  matrix  and  multiple  search  directions.  The  proposed algorithm  is  independent  of  the  global  coordinates  as  well  as  the  choice  of  bases.  All  genetic  operators  are performed  using  matrix  transformations,  and  hence,  are  computationally  efficient.  The  proposed  algorithm demonstrates state-of-the-art performance on three matrix manifold optimization  problems. Key words  Evolution  strategy;  matrix  manifold;  search  direction  adaptation;  black-box  optimization; covariance matrix adaptation     网络出版时间：2020-01-20 17:15:08网络出版地址：http://kns.cnki.net/kcms/detail/11.1826.TP.20200120.1640.008.html2  计  算  机  学  报  2020年  1  引言 矩阵流形优化又称黎曼优化，主要关注于定义在黎曼流形上的矩阵优化问题。这类问题可定义为： min∈ℳ   其中，为矩阵形式的决策变量，为目标函数，ℳ为黎曼流形。通常，ℳ不是欧式空间甚至不是线性空间，但是在任意一点处都可以局部线性化为一个切空间，记为ℳ。在这一切空间中，可以定义关于的内积运算 ⋅,⋅ ，因此，ℳ为欧式空间。如上定义的问题广泛地存在于数据挖掘[1]–[3]、数值代数[4], [5]、信号处理[6], [7]、电子设计[8]等领域。例如，在使用谱聚类[9]进行无监督学习任务时，需要找到一组高维数据在低维子空间中的嵌入表示，该任务可归结为寻找使得  最大化的正交矩阵的优化问题，其中为数据点相似图的拉普拉斯矩阵，(⋅)表示矩阵的迹。由于决策变量需要满足正交约束，因此其可行域不是线性空间，而是欧式空间中的子流形。 虽然矩阵流形优化问题在现实中有着广泛的应用，但传统数值优化方法在求解这类问题时，通常面临着两方面困难：一是在优化过程中保持流形约束非常困难。传统的流形约束处理方法，例如使用罚函数、将流形参数化等，依赖于特定参数的选择或问题本身的结构。二是搜索区域的非凸性使得问题具有许多局部最优解，在大多数情形下传统算法无法保证得到全局最优解。 早在18 世纪关于微分流形的研究中，研究者就已经观察到矩阵流形具有的特殊几何结构可以用于构造流形优化算法。然而直到最近二十年，由于这类问题在学术研究和工程实践中的广泛存在，相关课题才在数值优化领域被系统性的提出。目前，流形优化算法主要基于两类框架：   1）基于测地线的框架。Edelman等人[10]最早在数值优化的范畴内研究了矩阵流形的几何结构，提出了矩阵流形优化的通用框架。其基本思想是沿着测地线（Geodesic）方向[11]迭代下降。测地线是直线在曲面中的扩展，其计算过程通常借助于微分流形的相关工具和技术。对于某些常用流形，测地线有解析式，因此所需计算量较少。此外，也可以使用伪测地线[12]代替测地线，进一步提升算法运行效率而不影响算法性能。然而，对于那些测地线不存在解析式的流形，测地线依然需要通过微分方程求解，因而运行效率低下，只能用于小规模问题。   2）基于回撤的框架。Manton[13]在微分几何的范畴之外重新设计了流形优化方法。其算法迭代地求解形如    的目标函数，其中:ℳ→ℳ为测地线的近似算子，因而将流形约束问题转换为了切空间中的无约束问题。近几年来，学者们提出了针对不同流形的不同近似算子[14]–[17]，并将它们统称为回撤（Retraction）。一个典型的回撤算子通常包含了在切空间中移动和投影回流形的两个基本步骤，其计算过程不依赖测地线，因而效率较高。此外，常见流形的回撤算子大多可以通过线性代数计算，因而适用面比基于测地线的方法更广。 目前常用的流形优化算法大多在上述框架中使用最速下降、共轭梯度、牛顿法等经典方法[5], [18]–[21]产生解。文献[14], [22], [23]对这些矩阵流形上的经典算法进行了细致地综述。由于矩阵流形局部同胚于欧式空间，因此这类算法大多能保证局部收敛性。同时，这些算法也存在一定缺陷：一方面，由于流形的非凸性，这些算法难以保证全局收敛性；另一方面，越来越多的现实问题被建模成非光滑或不可导的问题，因而经典方法往往难以直接应用。这些缺陷的存在对矩阵流形问题的黑盒优化方法提出了迫切需求。 演化算法（Evolutionary Algorithm, EA）[24]是黑盒优化的经典方法。EA不依赖于问题本身特性，具有较强的全局搜索能力，因此可以期望EA能在矩阵流形优化中有较好的表现。然而，现阶段的演化算法大多针对欧式空间中的优化问题而设计，当用于矩阵流形优化时，将导致以下两个方面问题： 1）大多数演化算法运行于向量空间而非矩阵空间，从而造成了演化流形优化的第一个核心困难：如何对矩阵进行有效编码。这一困难的本质在于矩阵流形往往不能占据整个矩阵线性空间，因而，算法在表示矩阵或某个搜索方向时所需的自由度小于矩阵大小。显然，若直接把矩阵编码为向量，则会引入冗余变量，造成计算资源的浪费；将矩阵参数化虽然可以避免这一问题，但往往依赖于问题本身特性，且计算量较大。 2）目前常用的演化算子大多依赖于决策空间的欧氏几何结构，引入了演化流形优化的第二个核心困难：如何设计演化算子。例如，差分演化（Differential Evolution, DE）[25]中的变异算子、演  何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  3  化策略（Evolution Strategy, ES）[26]中的重组算子在本质上都是线性变换，因而需要决策空间为线性空间。再比如，DE中的交叉算子以及遗传算法（Genetic Algorithm, GA）[27]中几乎所有的交叉和变异算子都隐含地使用到了坐标的概念，因此这些算子必须在欧式空间中才能运行。此外，协方差矩阵适应演化策略（Covariance  Matrix  Adaptation Evolution Strategy, CMA-ES）[28]中协方差矩阵的更新以及粒子群优化（Particle  Swarm  Optimization, PSO）[29]中粒子速度的更新，不仅需要欧式空间，还要求该欧式空间在算法的整个运行过程中具有一个全局坐标系从而保证更新的有效性。 目前，演化算法领域的研究者已对上述两个问题进行了初步研究。文献[30]在PSO的速度更新机制中引入了黎曼梯度，提升算法在求解低秩张量近似问题时的性能。文献[6],  [31]则对标准PSO进行了简单的扩展。在所提出的算法中，粒子的速度位于该粒子所决定的切空间中，每个粒子沿着相应的测地线方向前进。文献[31]将CMA-ES推广至三维球面空间，在球面上沿测地线方向移动候选解及相应的概率分布，求解核磁共振成像中的形状重建问题。该文献还在球面空间上实现了DE用作对比算法，其使用对数映射将可行解映射至切空间，在切空间中使用DE算子产生候选解后用指数映射投影回球面。文献[32]将文献[31]中的变异强度适应机制与标准ES相结合，并在球面空间中分析了算法收敛性能。文献[33]提出了变异强度自适应的标准ES，用于在三维欧式空间与三维球面空间组成的积流形上求解室内定位问题。 上述研究成果虽然在黎曼流形的黑盒优化任务上表现出了可观的应用前景，但仍然存在一些可改进的空间。例如，文献[6], [30], [31]中的PSO改进算法在移动粒子后，需要将粒子投影回切空间以更新对应的速度向量。这一过程需要计算测地线的逆映射，因而只在少数几类流形中可以实现，算法欠缺通用性。文献[31]将CMA-ES中的协方差矩阵视为黎曼流形上的双线性映射，在移动概率分布时需要进行矩阵分解，因此并未在本质上解决概率分布需要额外编码的问题。文献[31]中的改进DE算法还需要在切空间中进行显式的正交化过程，计算效率较低。 本文立足于解决现有演化算法在求解矩阵流形优化问题时存在的困境，在基于回撤的算法框架内，提出了流形搜索方向适应演化策略（Manifold Search  Direction  Adaptation  Evolution  Strategy, MSDA-ES）。该算法有以下主要特性和创新之处： 1）把经典的CMA-ES算法从欧式空间推广到黎曼流形之中，设计了在黎曼流形上维护和更新多元高斯分布的方法。 2）提出了先在切空间产生候选矩阵，后映射至流形进行评估的采样方案，所有解可使用矩阵直接表示，而演化操作则通过线性代数完成，算法不再需要额外的编码策略。 3）提出了使用少量切向量构造和更新协方差矩阵的方法，算法运行过程与特定的坐标系无关，不依赖于全局坐标系，也不需要显式地构造局部坐标系，运行效率较高。 本文的第二章首先介绍矩阵流形的基础知识；随后，在第三章给出MSDA-ES算法的具体实现；第四章在一组测试问题中检验所提算法的有效性；最后，在第五章总结全文并给出未来的研究方向。 2  研究基础   本节介绍矩阵流形优化问题的几个具体例子，同时，给出了流形几何结构相关的结论和方法，这些基础知识将在下一章中用于算法设计。 2.1  三类常见的矩阵流形优化问题 本研究考虑自变量为×矩阵的黎曼流形优化问题，所涉及的流形结构包括Stiefel 流形、Grassmann流形以及Oblique流形。需要注意的是，本研究所提出的流形演化策略采用通用的流形优化框架，因此其适用范围不局限于这三类流形。 1）Stiefel流形。首先关注以下定义在×正交矩阵集合上的函数  最小化问题： min∈×   s.t.= 其中，表示×单位阵。这一问题的可行域称为Stiefel流形[34]，记作 , 。 2）Grassmann流形。本文要研究的第二个问题非常类似于上述第一个问题，其主要区别是，该问题的目标函数值表现出旋转不变性。该问题可通过以下方式定义： min∈×  =   s.t.= 其中，为任意的×正交矩阵。以上定义可知，4  计  算  机  学  报  2020年  目标函数的值仅依赖于的列空间，而和基的选取无关。这一问题的可行域称为Grassmann流形，记为 , ，该问题与子空间追踪等多种实际应用有紧密的联系[14]。 3）Oblique流形。本章关注的最后一个例子具有以下形式 min∈×   s.t.  = 其中，符号  表示将方阵中所有非对角线元素设为0时所得的对角阵。这一问题的可行域常常被称为Oblique流形，记作 , 。许多二值整数规划问题经过松弛后都可以转化为此类问题[12]。 2.2  矩阵流形的基础知识 本节简要介绍有关上述三类流形的常用符号和结论。需要注意的是，本节中出现的所有图例仅考虑了Stiefel流形。相对而言，Grassmann流形和Oblique流形更加抽象，无法使用图片展示其几何结构。然而，本节中给出的维度计算方法、切空间投影、回撤算子等内容对上述三类流形均有效。至于Grassmann流形和Oblique流形的具体性质，读者可在文献[10], [14]中找到。 2.2.1  切空间与法空间   对于流形ℳ上的任意一点，都存在一线性空间以为原点、与ℳ相切于。这一空间称为（关于）的切空间，记为ℳ。本研究假设ℳ为黎曼流形，因此切空间同时为欧式空间。将切空间ℳ的正交补称为法空间，记为ℳ。图  1 给出了 3,1 流形的切空间和法空间的例子。该流形为嵌入在3维空间中的一个球面，其任意一点对应的切空间为球的2维切平面，而相应的法空间为经过这一点的法线。  2.2.2  流形的维度   在演化算法中，所涉及的诸多参数往往和可行域的维度有关。因此，在进一步设计流形演化算法前，需要了解所求流形优化问题的维度。一个流形的维度等于该流形所对应切空间的维度，亦即参数化切空间所需自由元素的个数。对于上述三类流形，其维度计算方式列于表  1中。 2.2.3  切空间投影与回撤运算   对于矩阵流形ℳ上的点，可以定义从矩阵空间到切空间的正交投影:×→ℳ。为了方便起见，在投影时将视为矩阵空间×的原点。对于上述三类矩阵流形，表  1给出了相应的切空间投影计算方法。 同样，对于切空间中的任意一点，都可以定义表  1  三类矩阵流形的基本性质 ℳ  特性  ,    ℳ=−12 +1      =−12 +      = +     ,    ℳ= −      =−     =，其中和分别为矩阵+的SVD分解中的左奇异矩阵和右奇异矩阵  ,    ℳ= −1    = :,1 , ⋯, :,  ，其中  :, =:,−:, :,:,    = :,1 , ⋯, :,  ，其中  :, =:, +:, |:, +:, |2 Ÿ   ⋅ 表示QR分解中得到的正交矩阵； Ÿ  表中使用的QR分解和SVD分解均为精简型（economy-sized）分解； Ÿ  |⋅|2表示向量的2范数； Ÿ  :,表示矩阵中的第列。 图  1   3,1 流形及相应的切空间和法空间 图  2  切空间投影与回撤算子   何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  5  一个将该点映射回流形的变换，称为回撤（Retraction），记为:ℳ→ℳ。回撤算子可定义为满足如下两个条件的连续可微映射[15]：1）对任意的∈ℳ满足 0 =，其中0表示ℳ的原点；2）对任意的∈ℳ满足   =0=。理论上，测地线也可视为一个回撤算子，但通常认为回撤算子是测试线函数的一个近似。此外，对于同一类流形可以定义多种不同的回撤算子，本章使用文献[14],  [35]描述的定义方法，其具体的表达式在表  1中给出。   图  2给出了一个2维流形中切空间投影与回撤算子的示意图。图中，蓝色直线表示切空间投影，该投影将位于×中的点映射至ℳ。棕色曲线表示回撤算子，其将ℳ中的点映射至ℳ。 2.2.4  CMA-ES及其局限性   CMA-ES是求解连续型黑盒优化问题性能最好的演化算法之一。CMA-ES维持了一个多元高斯分布引导搜索，并且通过学习优质解的变量间依赖不断更新协方差矩阵，从而表现出对决策空间旋转变换的性能不变性。由于使用了演化路径技术对群体的成功变异方向进行累积，只需要较小的群体规模即可实现较快的协方差矩阵更新速度。CMA-ES还采用了累积步长控制技术（Cumulative  Step-size Adaptation, CSA），使得变异强度的适应机制独立于特定坐标系。此外，由于CMA-ES算法的环境选择过程仅依赖于目标函数值的比较结果，因此算法具有保序变换不变性。这些优良特性使得CMA-ES在许多现实中有着广泛的应用[36]–[38]。 如前文所述，CMA-ES存在一个显著的局限性，即其只能运行于欧式空间。例如，算法维持的协方差矩阵仅在欧式空间中才有定义，而协方差矩阵的更新还依赖于全局坐标系。对于流形优化问题，一方面，决策变量间的线性关系仅在切空间中的局部区域有定义，另一方面，可行域并不存在一个全局坐标系。这些性质使得CMA-ES无法直接应用于矩阵流形优化。文献[31]通过在流形上平行移动双线性映射实现了协方差矩阵的更新，通过对切向量的平行移动实现了演化路径的更新。然而，该方法至少在两方面存在可改进的空间。首先，双线性映射的平行移动需要对矩阵进行分解，隐含地引入了对流形参数化的过程，因此运行效率较低，无法应用于大规模环境；其次，将演化路径累积机制直接推广至黎曼流形后的几何意义不明确，事实上，文献[32]已经从理论上证明，在黎曼流形上直接使用CSA技术无法确保算法收敛。 3  流形搜索方向适应演化策略 3.1  基本思想 MSDA-ES是CMA-ES在矩阵流形中的推广。它使用标准的 , -ES算法框架，并利用一个多元高斯分布引导群体搜索。在每一次迭代中，算法在某一切空间中产生出个解，利用回撤算子将解映射至流形从而完成函数评估，最后，使用其中最好的个解更新搜索分布。 假设矩阵流形ℳ的维度为，即=ℳ。MSDA-ES在其第代运行中维持了以下策略参数： Ÿ    ∈ℳ：矩阵流形上的一个可行解； Ÿ    ∈+：变异强度； Ÿ    ∈×：协方差矩阵。 上述协方差矩阵  用于描述切空间  ℳ中决策变量间的线性依赖。使用这些参数，构造如下的多元高斯分布：     2 0,    其中，0表示维0向量。由于  是切空间  ℳ的原点，因此服从上述分布的随机变量自然地分布于  周围。因此，虽然在上述分布中并没有指明均值，但  实际上起到了均值的作用。 在第代中，MSDA-ES的工作流程如下： 1.  从高斯分布    2 0,   中采样出个位于  ℳ中的候选矩阵1,⋯,。 2.  对候选矩阵排序，使得  1:  ≤⋯≤  :  。其中，:表示排名第的候选矩阵。 3.  对个最好的候选矩阵重组以得到新的矩阵。将使用回撤算子映射至流形作为下一代群体均值 +1 。 4.  使用更新搜索分布。 5.  将分布移动至新的切空间 +1 ℳ，同时计算σ +1 和 +1 。 上述步骤在以下几个方面区别于标准CMA-ES： Ÿ  CMA-ES在整个决策空间中维护概率分布，产生的候选解也分布于整个决策空间，而在MSDA-ES中，概率分布的维护以及候选解的6  计  算  机  学  报  2020年  采样均在切空间中完成。其原因在于流形与切空间是局部同胚的，因此在法空间中进行搜索只会浪费计算资源。 Ÿ  CMA-ES产生候选解后可直接进行函数评估，而MSDA-ES采样得到的候选矩阵并非可行解，因此需要在函数评估前使用回撤算子映射回流形（见上述第二步）。 Ÿ  CMA-ES在一个固定的全局坐标系中维护其概率分布，而在MSDA-ES中，群体的演化过程伴随着切空间的变化，因此，其概率分布在不同代可能位于不同的切空间中。对于均值 +1 ，需要使用回撤算子将其变换至矩阵流形以确保其可行性（见上述第三步）。对于协方差矩阵以及变异强度，还需要将其移动至新的切空间（见上述第五步）。 3.2  编码策略 对矩阵进行有效的编码是使用传统演化算法求解流形优化问题时面临的主要困难。本节介绍MSDA-ES中对解和概率分布的表示方法，可以有效克服上述困难。 3.2.1  协方差矩阵建模 标准CMA-ES使用一个规模为×的协方差矩阵存储所有自变量之间的两两线性关联，因此具有 2 的时间和空间复杂度。由于在矩阵流形优化中，=  常常为一个较大的值，因此上述表示协方差矩阵的方法不适用于矩阵流形优化。 为了降低存储和更新协方差矩阵的复杂度，MSDA-ES中的协方差矩阵具有以下形式：   = 1−� +�       =1# 1  其中，1  ,⋯,  ∈为个向量，�为权系数。在算法中，向量  描述了可能产生高质量解的方向，因此将其称为搜索方向。 公式 1 使用        =1 估计协方差矩阵，当≪时，维护协方差矩阵所需复杂度可以显著降低，从而在性能和效率之间达到折衷。同时，由于该估计量的秩不超过，为了保证协方差矩阵的正定性，将其与单位阵进行线性组合，并使用�作为组合系数。该线性组合使得公式 1 具有Shrinkage估计器[39]的形式，因此在保持原有估计量特征向量不变的情况下，改善了其正定性，有助于描述高维空间中的变量依赖关系。需要注意的是，上述模型只用于说明解的采样原理，模型本身并不真正的出现在算法实现中。 3.2.2  采样方法 除了具有较低的算法复杂度，使用公式 1 对协方差矩阵建模还有另外一个优点，即，采样过程与坐标系无关。 假设{∈×|  =1⋯}为  ℳ的一组基。为了从该组基所确定的坐标系中产生出一个候选矩阵，可首先对 0,   进行概率采样： = 1−�0+ �   �=1 该过程产生的向量中包含了候选矩阵的坐标值，其中，0∈为服从标准高斯分布 0, 的随机向量，�∈为服从标准高斯分布 0,1 的随机数， ∈{1,⋯,}。随后，该坐标向量所对应的候选矩阵可表示为 =   1−��0,+ � ,  �=1 =1 = 1−� �0,=1+ �   ,  =1 �=1 其中，�0,和,  分别表示0和  中的第个元素。 令  =  ,  =1∈×，则  是一个定义在  ℳ中、坐标为  的矩阵。因此，在上述协方差矩阵的建模中，可以使用矩阵  代替向量  从而在采样过程中避免使用坐标基。由于  和  发挥了相同的作用，本文的后续部分将该矩阵也同样称为搜索方向。 对于  �0,=1 ，从直观上看，其计算过程中不可避免地涉及到了坐标基。然而，可以观察到其本质上为  ℳ中服从标准高斯分布的随机变量。由于标准高斯分布与特定坐标基无关，因此可以将其替换成同一空间中呈标准高斯分布的其他随机变量。另一方面，记∈×为每一个元素均独立地服从 0,1 的矩阵，可知为矩阵空间×中服  何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  7  从标准高斯分布的随机变量；又由于  ℳ为×的子空间，而  为相应的正交投影，因此    的分布为的分布的边缘分布，亦即，    为  ℳ中服从标准高斯分布的随机变量。因此，可使用    代替  �0,=1 完成采样过程。 综上所述，当考虑变异强度时，MSDA-ES使用以下公式产生一个候选矩阵 =    1−�    + �   �=1 # 2  其中，∈×为每个元素都独立地服从 0,1 的矩阵，  ∈×为矩阵形式的搜索方向，�∈为服从 0,1 的随机数， ∈{1,⋯,}。 3.2.3  解的表示方法   从公式 2 可知，候选矩阵的产生过程与切空间中的坐标基无关，因此，每一个候选矩阵可以直接表示成矩阵形式而无需额外的参数化过程。需要注意的是，这一方法在表示一个矩阵时引入了冗余的参数，但由于在常见的流形优化问题中通常为一个较小的数，因而该编码策略中消耗的额外存储空间几乎可以忽略。另一方面，使用该方法时，矩阵的采样过程仅需要标准的线性代数运算，因此其计算效率较高。 3.3  协方差矩阵适应 标准CMA-ES通过协方差矩阵适应提高在整个欧式空间中重现成功变异方向的似然性，与此对应，MSDA-ES中协方差矩阵适应的目标是在切空间中提高重现成功变异方向的似然性。然而在MSDA-ES中，该过程伴随着切空间的改变，需要进行如下的特殊处理。 3.3.1  确定切空间 在第代中，当产生了全部个候选矩阵后，使用回撤算子将它们映射至流形ℳ。对映射回流形的候选矩阵进行评估和排序，使得    1:  ≤⋯≤    :  ，其中，:表示排名第的候选矩阵。 此后，对排名最好的个候选矩阵进行线性组合，得到新的分布均值： = �:=1 其中，�1,⋯,�为组合系数且满足  �=1=1,�1≥⋯≥�>0。上述操作与标准CMA-ES更新群体均值的过程一致。需要注意的是，由于所有候选矩阵都位于  ℳ中，而根据定义  ℳ为欧式空间，因此，此处进行的线性组合是合法运算。 最后，使用回撤算子将从切空间  ℳ映射至流形ℳ，得到以下新的矩阵：  +1 =     上述得到的新矩阵 +1 将在下一代中作为概率分布的均值，其相应的切空间则在下一代中作为搜索空间。换言之，上述过程完成了一次切空间的变换。切空间的变换在MSDA-ES中至关重要，如果算法在搜索过程中保持切空间固定不变，则新产生的候选矩阵将逐渐远离初始切空间的原点，由于流形与其切空间的同胚关系仅仅是定义在局部区域的，因此切空间中的候选矩阵与流形上的可行解之间的对应关系将越来越弱。相反，使用以上方法在每一次迭代中确定一个新的切空间作为当前代的搜索区域，使得每一代中产生的候选矩阵都分布在原点的周围，能够有效避免上述问题。 3.3.2  计算成功变异方向 演化策略中，成功变异方向是指群体均值的移动方向。MSDA-ES通过学习成功变异方向的决策变量间依赖关系以不断更新其概率分布。该方向由下式给出： =    其中，=1/  �2=1。可以验证，在不存在选择压力时，∼ 0,   。这表明的方差为协方差矩阵的一个无偏估计，因而在后续步骤中，将使用更新搜索方向。 3.3.3  更新搜索方向 搜索方向适应机制的核心思想是利用成功变异方向分别构造出个协方差矩阵的无偏估计，分别用于更新个搜索方向。 （1）更新1  ：令1=，则1的方差为  的一个无偏估计，因而第一个搜索方向按如下方式更新： 1′= 1− 1  +  2− 1 其中，∈ 0,1 为搜索方向的更新率。当不存在选择压力时，若1  ∼ 0,    ，则1′∼8  计  算  机  学  报  2020年   0,   ，亦即上式为无偏更新。需要注意的是，上式与标准CMA-ES中演化路径的更新过程具有相同形式，因此MSDA-ES中所使用的第一个搜索方向即为演化路径。此外，由于1=，因而其方差为（关于成功变异方向的）极大似然估计，即，以上更新得到的协方差矩阵增大了重现成功变异方向的似然性。   （2）更新2  ,⋯,  ：同样使用成功变异方向对余下的−1个搜索方向进行更新。然而，若使用与上述完全相同的更新方式，则这些搜索方向将变得越来越接近，导致概率分布的退化。为了保持搜索方向之间的独立性，为余下的每个搜索方向  构造一个新的成功变异方向；在完成相应的更新后，从中减去其在已更新搜索方向（记为′）上的投影，用以构造+1从而更新下一个搜索方向。换言之，算法使得每个+1都正交于′，因此，所得的+1′也将趋于与′保持独立，增强了在高维空间中估计协方差矩阵的鲁棒性。 具体而言，在更新第 个搜索方向后（ ∈{1,⋯,−1}），计算与其在′上的投影间的残差： −⋅ ′ 其中，= ,′   / ′,′   为投影的相对长度， ⋅,⋅   为流形在   的切空间中导出的内积运算。 上述投影残差已经与′正交，但由于其概率分布在运算中发生了改变，因此无法直接用于更新下一个搜索方向。为了避免这一问题、提升协方差矩阵更新的无偏性，对投影残差进行如下缩放： +1=1 1+2 −⋅ ′  可以验证，缩放系数1 1+2使得在满足下列三个条件时+1∼ 0,   ：1）不存在选择压力，2）∼ 0,   ，3）′∼ 0,   。即，+1的方差同样是  的无偏估计。利用这一性质，第+1个搜索方向可按如下方式更新： +1′= 1− +1  +  2− +1 可见，在更新个搜索方向时，算法在保证搜索方向无偏更新的同时减少了搜索方向间的依赖关系。 3.3.4  移动协方差矩阵 在上述过程中，更新后的搜索方向′位于当前代的切空间  ℳ中，为了能够使协方差矩阵用于下一代搜索过程，需要将搜索方向移动至由 +1 所决定的新的切空间。为此，可对其使用如下的正交投影：  +1 = +1  ′  上述所得的 +1 矩阵可用于在新的切空间中重构协方差矩阵。 图  3 给出了上述整个协方差矩阵适应过程的示意图。图中，  ℳ中的黑色箭头表示当前代搜索方向  ，使用当前代的成功变异方向对其进行更新后，产生了一个新的矩阵′（紫色箭头标注）。在矩阵空间×中将′的起点平行移动（使用蓝色虚线表示）到 +1 ，将其投影至 +1 ℳ从而得到新的搜索方向 +1 。 +1 将用于在下一代中重构协方差矩阵，图中可见，新的协方差矩阵的等密度轮廓线在更新后向着 +1 方向进一步延伸，因此提高了重现成功变异方向的似然性。 3.4  算法实现 本节给出MSDA-ES的具体实现细节。需要注意的是，虽然上述章节已经完整的描述了协方差矩阵适应过程，但还需要额外的机制对变异强度  进行适应。考虑到MSDA-ES使用了非标准的低秩模型对协方差矩阵进行建模，本节选择文献[40]中提出的基于群体的成功率准则（Population  Success Rule, PSP）完成变异强度适应。 图  3  协方差矩阵适应   何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  9  3.4.1  伪代码 算法I 给出了MSDA-ES算法的具体实现。在算法开始前，将矩阵  初始化至流形ℳ的任意位置，将搜索方向  初始化为0矩阵，将用于PSR的统计量  初始化为0，将变异强度  设置为1（第2-3行）。 在第代，算法首先利用公式 2 产生个候选矩阵，并将这些矩阵映射至ℳ后按目标函数值排序（第5-7行）。排序结果中较好的个候选矩阵通过线性组合得到新的群体均值和成功变异方向（第8-10行）。 随后，按照第3.3节中描述的方法进行协方差矩阵适应（第11-16行）。具体而言，首先使用成功变异方向对搜索方向  进行更新，再计算在更新后的搜索方向（即′）上的相对投影长度，最后对在′上的投影残差进行缩放，作为新的成功变异方向+1用于更新下一个搜索方向。 最后，在算法的第17-19 行，使用文献[40]中介绍的PSR对变异强度进行适应。该方法首先将当前代与上一代群体混合，对混合群体按目标函数值从小到大排序，分别计算两代群体在混合群体中的秩和。将秩和变化值归一化后与预设值�∗进行比较并进行累积，得到 +1 ，其中为累积率。上述步骤不断重复，直至满足终止条件。 3.4.2  参数设置 、、�等参数来源于标准CMA-ES，因此使用文献[41] 中 的 方 法 对 其 设 置：=4+ 3ln   ,= 0.5 ,�= +1 −   +1 −    =1。�为协方差矩阵的更新率，由于文献[42]指出，在信息几何优化框架中，协方差矩阵以 1/  的速度收敛至目标函数的Hessian矩阵的逆，因此将其设置为�=0.4/ 。为搜索方向的学习率，其合理的值应与成功变异方向的自由度成反比，因而将其设置为=0.25/ 。和�∗为变异强度适应机制PSR中引入的参数，因此按文献[40]中介绍的原则进行设置：=0.3,�∗=0.25。为使用的搜索方向的个数。通常，越大，算法能够捕获的有效搜索方向越多，因而性能越好；然而，使用过多的搜索方向也将增加算法复杂性。综合以上考虑，将其设置为=10。需要注意的是，上述部分参数虽然来源于其他算法，但其原始文献中的设置是为欧式空间中的优化而设计的，因此，当把这些参数应用于MSDA-ES时，这些参数中涉及到的问题维度应该设置为流形维度（即=ℳ）而非自变量矩阵的大小（即）。 3.4.3  时间复杂度 算法的时间复杂度由算法执行的回撤算子和切空间投影的次数决定。在算法的每一代中，产生每一个候选矩阵和更新每一个搜索方向各需执行一次切空间投影。评估每一个候选矩阵和产生新的群体均值各需执行一次回撤算子。因此，在使用第3.4.2节中推荐的参数时，算法平摊至每一次函数评估的时间复杂度，在Stiefel流形和Grassmann流形上为 2 ，在Oblique流形上为  。 算法I.    流形搜索方向适应演化策略MSDA-ES 输入：ℳ：矩阵流形；:×→：目标函数；:ℳ→ℳ：回撤算子；:×→ℳ：切空间投影。 输出：已找到的最优解。 1:    ←0   2:    ←ℳ中随机分布的矩阵,  σ  ←1 3:    ←0,    ←0×,   ∈{1, ⋯, } 4:  WHILE(未满足终止条件){ 5:  1, ⋯,  ←公式 2 采样 6:    ={    1:λ  , ⋯,     λ:λ  }   7:  将  按升序排列 8:  ←  �:=1   9:   +1 ←       10:  1= μ/σ     11:  FOR (←1:m){ 12:          ′= 1−   +  2−  13:          = , ′   / ′, ′    14: +1=1 1+2 −⋅′    15:   +1 = +1  ′    16:      } 17:       ,  −1 ←  ,  −1 在  ∪ −1 中的秩和 18:   +1 ← 1−σ   +σ  −1 −  λ2−�∗    19:  σ +1 ←σ  �  +1     20:      ←+1 21:  }  10  计  算  机  学  报  2020年  3.4.4  不变性分析 不变性是设计演化策略时需要重点考量的性质之一。不变性是指对目标函数或决策空间施加某些变换后，算法的运行行为或最终结果保持不变。当演化策略具有不变性时，其在某一个具体问题上展现出的性能可以推广至一大类问题，表明该算法具有普适性。 本文所提出的MSDA-ES与标准CMA-ES类似，具有目标函数保序变换不变性。具体而言，若表示任意严格单调递增函数，则MSDA-ES在求解min∈ℳ()和min∈ℳ    时将表现出完全相同的性能。其原因在于MSDA-ES并未直接利用目标函数值，而是利用了不同解之间的目标函数值的大小关系。由于函数是严格单调递增的，在群体中不改变目标函数值的序关系，因此算法性能与无关。 另一方面，当决策空间具有某些不变性时，MSDA-ES能够维持这种不变性。以Grassmann流形上的优化问题为例，该流形具有旋转不变性，因此当对某一优化问题的决策变量施加旋转变换后，虽然问题的形式发生了变换，但MSDA-ES依然将其视为同样的优化问题，展现出一致的算法性能。以下命题描述了MSDA-ES在求解Grassmann流形上的优化问题时所具有的决策空间旋转变换不变性。  命题1.  令�1和�2分别为算法I求解优化问题min∈(,)()和min∈(,)   =  时的两个算法实例，其中∈×为正交矩阵。令()=((),1(),…,(),σ(),()) 和  ()=( (), 1(),…, (),σ (), ())分别表示�1和�2运行在第 代 时 的 状 态 。 令 () =((),1(),…,(),σ(),()) 表示状态变换函数，其中∈×为任意正交矩阵。假设�1在第代使用随机变量  ,�() ( =1,..,)产生候选矩阵，则当 (0)= (0) 且�2在第代使用随机变量 ()=  ,� ()=�() ( =1,..,)产生候选矩阵时，有  ()= ()  对任意=1,2,…成立。  命题1的证明见补充材料。命题1表明，当对Grassmann流形优化问题的决策变量施加某些旋转变换时，只要算法的初始条件进行相应的变换（例如， (), ()分别与(),()具有相同列空间，而σ (), ()分别与σ(),()保持一致），则算法的后续行为将保持一致。需要指出的是，虽然命题1需要假设 ()=  ,� ()=�()  ( =1,..,)，但由于这些随机变量独立地采样自标准正态分布，因此即使在不满足上述假设时， ()与 () 的概率分布也是相一致的。此外，命题1描述的不变性的本质，在于旋转变化在Grassmann流形上为等价关系，因此命题1可直接推广至任意流形与给定等价关系导出的商流形，即，对任意商流形，只要算法初始状态是等价的，则对决策变量施加的任意等价变换不影响算法运行行为。 4  仿真实验 4.1  测试问题 4.1.1  基于割线的数据降维 基于割线的降维（Secant-based  Dimension Reduction,  SDR）用于在一个高维数据集中找到一个线性子空间，使得相应的正交投影尽可能地成为单射函数[43], [44]。 令∈×为样本容量为的数据矩阵。的单位割线集合为Σ={:,−:,|:,−:,|2|=1,⋯,}。令表示所求的子空间，那么当Σ与的正交补不存在交集时，对应的正交投影为单射。因此，该问题的优化任务是在Grassmann流形上找到一个正交矩阵，该矩阵张成且最小化以下函数： min∈ ,   =−min∈Σ||2   本实验中，数据矩阵通过=Λ产生。其中，∈×为一随机正交矩阵，Λ∈×为对角阵且对角线元素Λ,=�1−，∈×为随机矩阵且每个元素独立地服从 0,1 。参数�控制数据矩阵的条件数，在实验中固定为1.05。从集合{50,100,150,200,250,300}中选择，从{3,5,7,9,11}中选择，共计得到30组不同的测试用例。   何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  11  4.1.2  Thomson问题 Thomson问题为著名的3维装箱问题，其在数学、化学、物理等领域有着广泛的应用[45]。这一问题需要在一个球面上安放个电荷，使得总静电能最小。该问题的定义如下： min∈ ,   = |:,=1|2−1 其中，固定为3。Thomson问题的目标函数虽然可微，但存在着大量局部最优解。目前，仅在取几个特殊值时该问题存在有效解法，而在更一般的情况下则难以求解。在实验中，设置=25, ∈{1,2,⋯,20}，因此，共产生20组测试用例。 4.1.3  二次最小化问题 二次最小化问题（Quadratic Minimization, QM）是科学计算和数据挖掘中的常见问题[17]，其试图找到一个×正交矩阵从而最小化如下函数： min∈ ,   =12 � +   其中，�∈×为对称矩阵，∈×为任意矩阵， ⋅ 表示矩阵的迹。线性特征值问题为该问题在=0时的特例。该问题的目标函数不具有旋转不变性，因而其可行域为Stiefel流形而非Grassmann流形。 在该问题中，使用�=Λ随机产生矩阵�，其中∈×为随机正交矩阵，Λ∈×为对角阵且Λ,=�1−。矩阵由=得到，其中∈×为随机矩阵，其每一个元素独立地采样自 0,1 上的均匀分布，∈×为对角阵且,=ζ−1/|:,|2 。使用�=ζ=1.05，且从∈{50,100,150,200,250,300}以及∈{3,5,7,9,11}中产生出30组不同的测试用例。 4.2  实验设置 4.2.1  用于比较的算法 本研究选择最速下降法（Steepest  Descent Method, SD）、PSO、DE以及  CMA-ES作为对比算法。其中，SD和PSO为开源工具箱Manopt[35]中的内置求解器。DE和CMA-ES根据文献[31]中的方法实现，并进行了一定的修改以增强性能，读者可在附录中找到其伪代码。这些算法的设置如下： （1）SD：该方法在基于回撤的流形优化框架中使用经典的梯度下降法[14]。算法在每个切空间中使用Armijo准则[46]进行线搜索。由于本文内容针对黑盒优化，因此算法中所需的梯度通过有限微分法估计得到。具体而言，算法首先在矩阵空间×中使用前向微分估计目标函数的梯度，之后使用切空间投影将梯度映射至切空间。由于本文所考虑的三类流形均为嵌入子流形，因而所得的投影即为黎曼梯度。在上述过程中，估计一次梯度消耗+1次函数评估。该算法的所有参数按Manopt工具箱默认设置。 （2）PSO：该方法最早来源于文献[6]，其在标准PSO中引入测地线从而在3维旋转群中进行启发式搜索。Manopt工具箱对其进行了扩展，使其在测地线可逆时作为一般性的流形优化算法。在本实验中，由于Stiefel流形上的测地线不可逆，因而仅在Grassmann流形和Oblique流形上测试该算法。该算法的所有参数按Manopt工具箱默认设置。 （3）DE：按照文献[31]使用指数映射及其逆映射将传统DE算法推广至矩阵流形。与文献[31]采用DE/rand/1/bin 策略不同，本文采用了DE/rand/1/either-or策略，即，对每个解按概率执行变异操作或按概率1−执行算术交叉操作。该策略使所有演化算子均为线性变换，不再需要计算每个解的坐标值以进行传统的二进制交叉操作。依据文献[47]，设置=0.4，差分向量缩放因子=0.7，算术交叉系数=0.85。与PSO类似，该算法依赖指数映射的逆，因此不用于Stiefel流形上的测试问题。 （4）CMA-ES：根据文献[31]，使用平行移动实现矩阵流形上演化路径的更新。然而，考虑到文献[31]只涉及小规模问题，而本文实验使用的测试问题维度较高，为了避免平行移动双线性映射中的耗时操作，引入了文献[48]中提出的有限内存矩阵适应技术。该技术是协方差矩阵适应机制在大规模环境中的变体，使用少量几个演化路径重建协方差矩阵，因此与MSDA-ES有相同的算法复杂度，保证了实验的公平性。为了提升算法在大规模环境中的收敛速度，算法参数按文献[48]设置。 4.2.2  终止条件 本实验在函数评估次数超过�时终止所有算法。对于SDR问题，�=50，对于其他问题，�=100。此外，由于不同算法使用了不同的框架，因而为每个算法单独设置以下的终止条件：对SD，算法在梯度的2范数小于10−6或Armijo步长小于10−8时终止算法；对PSO和DE，最后100代中最优目标函数值的改进值小12  计  算  机  学  报  2020年  于10−8时终止算法；对于MSDA-ES和CMA-ES，变异强度小于10−6时终止算法。 4.2.3  性能指标与统计方法 对于每一个算法，记录其所得的目标函数值作为性能指标，目标函数值越小表明算法性能越好。所有算法均独立重复20 次，使用中值和四分位距衡量结果的平均水平和浮动范围。在每一个测试用例上，使用置信水平为0.05的Wilcoxon秩和检验比较MSDA-ES和相应算法之间的统计差异。使用多重检验比较所有算法在所有测试用例上的综合性能。该检验首先根据Friedman测试计算每个算法的平均排名，并使用Bonferroni后置过程进行校正，校正后的置信水平作为最终检验结果。 除了直接比较数值结果，还使用数据剖面图（Data profiles）对所有算法的综合性能做出可视化的对比。在数据剖面图中，横坐标（�）表示某个算法所分配的计算资源，纵坐标（ � ）表示在给定�值时算法能够求解的问题的百分比。在绘制剖面图时，对于一个算法，若其在一个测试用例上满足≤+10−3 �− ，则称该算法可在次函数评估内求解，其中为该算法在次函数评估内所获得的最优目标函数值，�为所有对比表  2  在SDR问题上MSDA-ES等算法最终结果的中值和四分位距。每个测试用例上的最优值和次优值分别用深色背景和浅色背景标注。 n  p  MSDA-ES  SD  PSO  DE  CMA-ES 50  3  -3.02E-1(3.2E-2)  -8.60E-2(2.7E-2)+  -1.21E-1(2.3E-2)+  -1.28E-1(1.7E-2)+  -2.18E-1(3.9E-2)+ 50  5  -5.52E-1(2.4E-2)  -1.82E-1(4.0E-2)+  -2.22E-1(2.1E-2)+  -2.37E-1(1.3E-1)+  -4.34E-1(2.9E-2)+ 50  7  -7.11E-1(3.6E-2)  -2.51E-1(4.4E-2)+  -2.93E-1(2.3E-2)+  -3.16E-1(1.6E-1)+  -5.56E-1(5.0E-2)+ 50  9  -8.05E-1(3.5E-2)  -3.25E-1(4.7E-2)+  -3.71E-1(2.2E-2)+  -3.83E-1(2.1E-1)+  -6.94E-1(4.2E-2)+ 50  11  -8.69E-1(1.1E-2)  -3.66E-1(8.0E-2)+  -4.18E-1(2.3E-2)+  -4.26E-1(2.1E-1)+  -7.62E-1(1.7E-2)+ 100  3  -3.19E-1(5.3E-2)  -5.81E-2(1.8E-2)+  -8.48E-2(1.7E-2)+  -9.40E-2(1.9E-2)+  -2.16E-1(2.0E-2)+ 100  5  -5.63E-1(3.8E-2)  -1.31E-1(3.2E-2)+  -1.61E-1(1.4E-2)+  -2.04E-1(1.1E-1)+  -4.40E-1(3.6E-2)+ 100  7  -7.15E-1(2.6E-2)  -1.95E-1(5.8E-2)+  -2.09E-1(1.8E-2)+  -2.48E-1(1.3E-1)+  -5.88E-1(3.0E-2)+ 100  9  -8.05E-1(2.0E-2)  -2.33E-1(4.7E-2)+  -2.55E-1(2.1E-2)+  -2.94E-1(2.2E-2)+  -6.86E-1(2.6E-2)+ 100  11  -8.75E-1(1.3E-2)  -2.69E-1(4.5E-2)+  -2.92E-1(1.6E-2)+  -3.35E-1(2.1E-1)+  -7.68E-1(3.5E-2)+ 150  3  -3.08E-1(3.1E-2)  -5.66E-2(1.7E-2)+  -7.22E-2(4.6E-3)+  -7.97E-2(1.2E-2)+  -2.17E-1(4.2E-2)+ 150  5  -5.56E-1(2.6E-2)  -1.20E-1(2.8E-2)+  -1.28E-1(9.7E-3)+  -1.62E-1(8.1E-2)+  -4.39E-1(4.2E-2)+ 150  7  -7.11E-1(2.4E-2)  -1.56E-1(3.8E-2)+  -1.72E-1(1.6E-2)+  -2.00E-1(1.2E-1)+  -5.99E-1(3.9E-2)+ 150  9  -8.03E-1(1.9E-2)  -2.04E-1(3.7E-2)+  -2.09E-1(2.0E-2)+  -2.83E-1(1.3E-1)+  -6.82E-1(2.5E-2)+ 150  11  -8.71E-1(1.3E-2)  -2.26E-1(2.9E-2)+  -2.39E-1(2.2E-2)+  -3.02E-1(1.6E-1)+  -7.47E-1(1.6E-2)+ 200  3  -3.20E-1(2.0E-2)  -5.38E-2(1.9E-2)+  -5.90E-2(1.5E-2)+  -7.18E-2(5.1E-2)+  -2.22E-1(2.8E-2)+ 200  5  -5.47E-1(3.2E-2)  -1.00E-1(2.6E-2)+  -1.12E-1(1.9E-2)+  -1.46E-1(8.2E-2)+  -4.47E-1(4.1E-2)+ 200  7  -7.16E-1(1.7E-2)  -1.42E-1(3.4E-2)+  -1.51E-1(1.5E-2)+  -1.91E-1(1.1E-1)+  -6.11E-1(1.7E-2)+ 200  9  -8.04E-1(2.4E-2)  -1.78E-1(4.1E-2)+  -1.84E-1(2.0E-2)+  -2.03E-1(3.6E-2)+  -6.97E-1(3.3E-2)+ 200  11  -8.66E-1(2.0E-2)  -2.06E-1(3.3E-2)+  -2.07E-1(8.0E-3)+  -3.03E-1(1.6E-1)+  -7.54E-1(3.3E-2)+ 250  3  -3.31E-1(1.9E-2)  -4.89E-2(1.6E-2)+  -5.47E-2(7.7E-3)+  -6.83E-2(5.4E-2)+  -2.21E-1(2.3E-2)+ 250  5  -5.67E-1(1.9E-2)  -8.41E-2(3.9E-2)+  -1.02E-1(1.5E-2)+  -1.20E-1(8.2E-2)+  -4.62E-1(3.3E-2)+ 250  7  -7.10E-1(3.2E-2)  -1.31E-1(3.4E-2)+  -1.33E-1(1.4E-2)+  -1.64E-1(9.9E-2)+  -6.09E-1(2.8E-2)+ 250  9  -8.03E-1(1.3E-2)  -1.62E-1(1.8E-2)+  -1.63E-1(1.3E-2)+  -1.98E-1(1.1E-1)+  -7.04E-1(4.4E-2)+ 250  11  -8.67E-1(1.1E-2)  -1.90E-1(3.0E-2)+  -1.85E-1(1.3E-2)+  -2.09E-1(3.2E-2)+  -7.54E-1(2.3E-2)+ 300  3  -3.27E-1(3.2E-2)  -4.63E-2(1.6E-2)+  -5.32E-2(7.1E-3)+  -6.85E-2(3.8E-2)+  -2.30E-1(2.7E-2)+ 300  5  -5.64E-1(3.0E-2)  -7.95E-2(1.7E-2)+  -9.33E-2(9.8E-3)+  -1.43E-1(6.9E-2)+  -4.51E-1(2.8E-2)+ 300  7  -7.22E-1(2.3E-2)  -1.17E-1(2.4E-2)+  -1.21E-1(1.4E-2)+  -2.24E-1(9.5E-2)+  -6.07E-1(3.0E-2)+ 300  9  -8.07E-1(1.9E-2)  -1.47E-1(2.2E-2)+  -1.49E-1(1.8E-2)+  -1.74E-1(6.1E-2)+  -6.97E-1(3.6E-2)+ 300  11  -8.66E-1(1.1E-2)  -1.77E-1(1.7E-2)+  -1.70E-1(1.5E-2)+  -2.02E-1(8.7E-2)+  -7.61E-1(2.2E-2)+  + \ - \ =  30 \ 0 \ 0  30 \ 0 \ 0  30 \ 0 \ 0  30 \ 0 \ 0   “+”表示MSDA-ES显著优于（基于显著性水平0.05 的Wilcoxon 秩和检验）对比算法；“-”表示相反，即对比算法显著优于MSDA-ES；若不存在显著差异，则用“=”表示。这些符号在本文其他表格中具有相同含义。   何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  13  算法初始解中的最大目标函数值，为次函数评估内所有算法得到的最优函数值。数据剖面（曲线）可定义为  � =1    ∈,+1≤�   其中， ⋅ 表示集合的元素个数，表示测试用例的集合，表示测试用例的大小。设置为矩阵大小，即=，因此+1可解释为一次梯度估计所需的评估次数。这一缩放因子使得数据剖面独立于问题规模。在数据剖面图中，可以在给定�值时通过数据剖面曲线的位置衡量算法优劣。同时，可以计算数据剖面的曲线下面积（Area  under  the  Curve, �=  � �）衡量算法的综合性能，算法的�值越大表明算法的综合性能越好。 4.3  数值结果 4.3.1  SDR问题上的最终结果 表  2给出了MSDA-ES、SD、PSO、DE以及CMA-ES在基于割线的数据降维问题上的最终实验结果。显然，MSDA-ES在所有问题上均显著优于其他对比算法。SDR问题的难度主要来源于其不可微性，因此SD使用的有限微分法难以估计出精确的下降方向。此外，该问题的目标函数非凸，因而确定性算法在该问题上难以有效收敛至全局最优解。由于这两方面原因，SD的性能表现最差。PSO、DE和CMA-ES由于采用了演化算法框架，不依赖于问题的可微性或凸性，因而在该问题上性能优于SD。然而，上述三个演化算法依然在所有测试用例上均弱于MSDA-ES。PSO和DE性能较弱的原因可能在于其无法有效求解大规模问题[49],  [50]。事实上，目前绝大多数PSO和DE的相关研究都仅针对于100维以下的问题，而本文中问题的维度普遍在150至2000之间。目前使用PSO和DE求解大规模问题的研究大多依赖于对全局坐标系的分解[51]–[53]，因此无法应用于流形优化之中。相比之下，由于MSDA-ES和CMA-ES引入了专门设计的大规模协方差矩阵建模技术，因此在求解流形优化问题中有明显的优势。MSDA-ES相比于CMA-ES的性能优势主要来源于其搜索方向适应机制。正如文献[31],  [32]指出的，传统的演化路径累积技术在黎曼流形中的几何意义不明确，因而间接影响了协方差矩阵适应和变异强度适应机制，造成算法性能的降低。由此可见，MSDA-ES中专门设计的流形表  3  在Thomson问题上MSDA-ES等算法最终结果的中值和四分位距。每个测试用例上的最优值和次优值分别用深色背景和浅色背景标注。 n  p  MSDA-ES  SD  PSO  DE  CMA-ES 3  25  2.44E+2(8.0E-4)  2.44E+2(1.0E-3)=  2.46E+2(1.1E+0)+  2.57E+2(2.2E+0)+  2.44E+2(2.9E-2)+ 3  50  1.06E+3(8.7E-2)  1.06E+3(3.9E-2)+  1.08E+3(4.3E+0)+  1.11E+3(7.7E+0)+  1.06E+3(1.1E-1)+ 3  75  2.45E+3(9.4E-2)  2.45E+3(1.3E-1)=  2.52E+3(9.3E+0)+  2.59E+3(1.4E+1)+  2.45E+3(1.8E-1)+ 3  100  4.45E+3(2.5E-1)  4.45E+3(2.7E-1)=  4.57E+3(1.4E+1)+  4.68E+3(2.0E+1)+  4.45E+3(2.2E-1)+ 3  125  7.04E+3(1.6E-1)  7.04E+3(2.4E-1)+  7.24E+3(3.1E+1)+  7.38E+3(2.5E+1)+  7.04E+3(2.8E-1)+ 3  150  1.02E+4(2.5E-1)  1.02E+4(4.0E-1)+  1.05E+4(2.7E+1)+  1.07E+4(4.9E+1)+  1.02E+4(3.2E-1)+ 3  175  1.40E+4(2.4E-1)  1.40E+4(4.7E-1)+  1.45E+4(5.9E+1)+  1.47E+4(3.9E+1)+  1.40E+4(4.7E-1)+ 3  200  1.84E+4(3.7E-1)  1.84E+4(4.5E-1)+  1.90E+4(6.3E+1)+  1.93E+4(4.0E+1)+  1.84E+4(4.7E-1)+ 3  225  2.35E+4(4.0E-1)  2.35E+4(7.5E-1)+  2.42E+4(1.1E+2)+  2.45E+4(8.7E+1)+  2.35E+4(6.2E-1)+ 3  250  2.91E+4(4.0E-1)  2.91E+4(4.8E-1)+  2.99E+4(8.6E+1)+  3.03E+4(5.2E+1)+  2.91E+4(4.9E-1)+ 3  275  3.53E+4(3.3E-1)  3.53E+4(6.7E-1)+  3.63E+4(1.5E+2)+  3.67E+4(1.0E+2)+  3.53E+4(7.5E-1)+ 3  300  4.21E+4(4.2E-1)  4.21E+4(9.4E-1)+  4.33E+4(1.3E+2)+  4.38E+4(7.0E+1)+  4.21E+4(3.8E-1)+ 3  325  4.96E+4(5.1E-1)  4.96E+4(4.2E-1)+  5.09E+4(1.6E+2)+  5.15E+4(8.5E+1)+  4.96E+4(5.9E-1)+ 3  350  5.76E+4(6.6E-1)  5.76E+4(9.7E-1)+  5.92E+4(1.5E+2)+  5.99E+4(9.5E+1)+  5.76E+4(1.1E+0)+ 3  375  6.63E+4(6.2E-1)  6.63E+4(1.1E+0)+  6.80E+4(1.7E+2)+  6.88E+4(1.1E+2)+  6.63E+4(1.1E+0)+ 3  400  7.56E+4(5.3E-1)  7.56E+4(1.2E+0)+  7.76E+4(1.6E+2)+  7.84E+4(9.6E+1)+  7.56E+4(8.1E-1)+ 3  425  8.55E+4(1.1E+0)  8.55E+4(1.8E+0)+  8.78E+4(2.5E+2)+  8.85E+4(1.0E+2)+  8.55E+4(7.8E-1)+ 3  450  9.60E+4(1.1E+0)  9.60E+4(1.4E+0)+  9.85E+4(2.0E+2)+  9.94E+4(9.2E+1)+  9.60E+4(1.3E+0)+ 3  475  1.07E+5(7.1E-1)  1.07E+5(1.7E+0)+  1.10E+5(2.4E+2)+  1.11E+5(1.4E+2)+  1.07E+5(1.1E+0)+ 3  500  1.19E+5(9.4E-1)  1.19E+5(1.5E+0)+  1.22E+5(2.1E+2)+  1.23E+5(1.5E+2)+  1.19E+5(1.2E+0)+    + \ - \ =  17 \ 0 \ 3  20 \ 0 \ 0  20 \ 0 \ 0  20 \ 0 \ 0  14  计  算  机  学  报  2020年  搜索方向适应机制对于求解矩阵流形上的优化问题是必要的。 4.3.2  Thomson问题上的最终结果 由于Thomson问题包含了大量局部最优解，因此可以预期，使用演化策略框架的MSDA-ES能够在该问题上表现出良好的性能。表  3  给出的实验结果很大程度上验证了该观点。MSDA-ES在所有测试用例上均取得最优结果，且在几乎所有测试用例上显著优于所有对比算法。对于SD，由于目标函数具有可微性，因而其收敛速度较快，能够在所有测试用例中得到仅次于MSDA-ES的结果。PSO和DE虽然采用了演化算法框架，在避免陷入局部最优解方面存在一定优势，但由于收敛能力不足，因此最终结果相对较差。此外，我们发现MSDA-ES在大多数问题上对应的IQR 值较小，表明其在每一次独立运行中都能够产生出类似的解，说明了MSDA-ES采用的演化策略框架能够在保持收敛性能的前提下避免算法陷入局部最优解。CMA-ES的表现类似于MSDA-ES，但IQR值相对更大。其原因如前所述，可能是由于在黎曼流形上使用了针对欧式空间设计的演化路径累积技术。事实上，使用切空间投影在矩阵流形上移动演化路径时，不可避免地会缩短演化路径长度，导致变异强度过度缩小，造成算法提前收敛。MSDA-ES中的变异强度适应机制不依赖于演化路径或搜索方向，仅和群体的目标函数值有关，因此未发生类似的问题。 4.3.3  QM问题上的最终结果 QM问题为定义在Stiefel流形上的矩阵优化问题，由于在该流形中测地线不可逆，因此本组实验不对PSO和DE进行测试。相比于以上两类问题，QM问题中的目标函数是光滑的，因此使用梯度的算法能够有效求解该问题，同时，在满足某些假设时能够达到全局最优解[17]。因此，该组实验的主要目的并非是比较算法的优劣，而是验证MSDA-ES作为通用流形优化器的可能性。表  4给出了实验结果，可以观察到，MSDA-ES和CMA-ES性能类似，在2个测试用例上取得了显著优势，而在另外28 个测试用例上不存在显著差异。相比于SD，MSDA-ES在30个问题中的18个上获得了更好的最终结果，但同样在28 个测试用例上不存在显著的统计差异。因此，可以从该组实验中得出如下结论：在相对简单的问题上，MSDA-ES的性能不弱于具有理论支撑的算法。 4.3.4  收敛行为研究 上节中的实验仅反映了各个算法在终止后的最终结果，其实验结果可能受到终止条件的影响。因此，本节将使用收敛曲线图描述算法在特定测试用例上的目标函数值随评估次数的变化趋势。对于每一个测试问题，分别选择规模最大和最小的两个测试用例。用于绘图的数据来源于20 次独立运行中结果最接近于中值的那一次。 图  4、图  5、图  6给出了各算法的收敛曲线，有如下观察结果： Ÿ  MSDA-ES的收敛速度很快，在大多数测试问题上，能够在约20%∼30%的总评估次数内接近最终解。在图  4(b)和图  5(b)所示的高维问题表  4  在QM问题上MSDA-ES等算法最终结果的中值和四分位距。每个测试用例上的最优值和次优值分别用深色背景和浅色背景标注。 n  p  MSDA-ES  SD  CMA-ES 50  3  -1.98E+0(5.9E-2)  -1.98E+0(6.7E-2)=  -1.99E+0(7.0E-2)= 50  5  -3.36E+0(1.7E-1)  -3.32E+0(7.9E-2)=  -3.37E+0(9.3E-2)= 50  7  -4.78E+0(1.2E-1)  -4.75E+0(1.2E-1)=  -4.83E+0(1.4E-1)= 50  9  -6.42E+0(1.8E-1)  -6.41E+0(1.7E-1)=  -6.37E+0(1.3E-1)= 50  11  -8.12E+0(1.9E-1)  -8.01E+0(3.0E-1)=  -8.13E+0(1.1E-1)= 100  3  -2.16E+0(4.8E-2)  -2.14E+0(5.9E-2)=  -2.17E+0(3.4E-2)= 100  5  -3.69E+0(6.0E-2)  -3.69E+0(6.7E-2)=  -3.67E+0(6.5E-2)= 100  7  -5.31E+0(8.1E-2)  -5.34E+0(6.7E-2)=  -5.30E+0(1.1E-1)= 100  9  -7.07E+0(1.2E-1)  -7.06E+0(8.3E-2)=  -7.03E+0(1.8E-1)= 100  11  -8.93E+0(1.8E-1)  -8.98E+0(1.1E-1)=  -8.96E+0(1.3E-1)= 150  3  -2.21E+0(3.5E-2)  -2.23E+0(3.0E-2)-  -2.22E+0(4.1E-2)= 150  5  -3.80E+0(4.4E-2)  -3.82E+0(5.4E-2)=  -3.80E+0(3.9E-2)= 150  7  -5.48E+0(1.1E-1)  -5.48E+0(8.5E-2)=  -5.47E+0(6.3E-2)= 150  9  -7.36E+0(1.5E-1)  -7.32E+0(6.0E-2)=  -7.32E+0(7.2E-2)= 150  11  -9.27E+0(1.0E-1)  -9.25E+0(1.3E-1)=  -9.26E+0(6.2E-2)= 200  3  -2.25E+0(3.5E-2)  -2.25E+0(3.2E-2)=  -2.24E+0(2.7E-2)= 200  5  -3.85E+0(5.1E-2)  -3.86E+0(4.8E-2)=  -3.86E+0(5.2E-2)= 200  7  -5.58E+0(4.3E-2)  -5.56E+0(4.5E-2)+  -5.57E+0(5.8E-2)= 200  9  -7.43E+0(5.9E-2)  -7.41E+0(1.0E-1)=  -7.38E+0(7.8E-2)+ 200  11  -9.40E+0(1.1E-1)  -9.41E+0(1.1E-1)=  -9.40E+0(1.1E-1)= 250  3  -2.26E+0(1.7E-2)  -2.27E+0(3.3E-2)=  -2.27E+0(3.2E-2)= 250  5  -3.89E+0(4.9E-2)  -3.89E+0(4.9E-2)=  -3.88E+0(4.4E-2)= 250  7  -5.61E+0(4.8E-2)  -5.61E+0(3.7E-2)=  -5.61E+0(4.1E-2)= 250  9  -7.49E+0(4.6E-2)  -7.48E+0(3.6E-2)=  -7.46E+0(5.4E-2)+ 250  11  -9.51E+0(8.1E-2)  -9.49E+0(4.7E-2)=  -9.49E+0(8.7E-2)= 300  3  -2.28E+0(2.3E-2)  -2.28E+0(2.5E-2)=  -2.27E+0(3.0E-2)= 300  5  -3.90E+0(3.4E-2)  -3.91E+0(3.3E-2)=  -3.89E+0(3.6E-2)= 300  7  -5.65E+0(3.8E-2)  -5.63E+0(5.4E-2)=  -5.65E+0(5.4E-2)= 300  9  -7.52E+0(9.0E-2)  -7.52E+0(7.0E-2)=  -7.52E+0(7.0E-2)= 300  11  -9.55E+0(7.9E-2)  -9.54E+0(6.1E-2)=  -9.55E+0(8.1E-2)=    + \ - \ =  1 \ 1 \ 28  2 \ 0 \ 28    何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  15   (a) =50,=3  (b)  =300,=11 图  4 SDR问题上的收敛曲线 中，其收敛速度甚至超过SD。CMA-ES具有类似的收敛行为，而它们的收敛速度都比PSO和DE更快。这一观察结果与欧氏空间优化的相关研究基本一致[54]。 Ÿ  MSDA-ES在所有测试用例中具有相似的收敛行为。与之相反，SD对目标函数的可微性高度敏感，其在不可微的SDR问题上几乎在初始化后立刻就陷入了停滞状态，而在其他问题上的收敛速度则相对较快。DE和PSO的收敛行为受到问题类型和问题规模的影响，例如，由图  5可知，PSO和DE在问题规模变大时更易提前收敛。 Ÿ  在相对较简单的低维QM问题上，MSDA-ES和CMA-ES的收敛速度较SD更慢，但SD的速度优势在问题维度升高后则不再显著。其可能的原因在于SD中梯度估计时所产生的误差在运行中不断累积，导致SD在接近最优解时难以达到较高的精度。      4.3.5  综合性能 (a) =3,=25  (b)  =3,=500 图  5 Thomson问题上的收敛曲线  (a) SDR  (b) Thomson  (c) QM 图  7  数据剖面图 16  计  算  机  学  报  2020年  上述实验仅反映了算法在特定测试用例上的性能，本节使用数据剖面图和多重检验测试所提出的算法在所有测试用例上的综合性能。 （1）数据剖面图。使用数据剖面图（见图  7）展示各个算法在给定问题能够成功求解的测试用例的比例。由于PSO和DE无法应用于QM问题，而在SDR和Thomson问题上也无法成功求解任何一个测试用例，因此它们并未绘制在数据剖面图中。 SDR 问题（见图  7(a)）具有不可微的目标函数，因此在三个测试问题中求解难度最大。除 MSDA-ES  以外的所有对比算法（包括同样使用演化策略框架的CMA-ES）均无法求解任何一个测试用例。因此，仅有  MSDA-ES的数据剖面曲线反映在该图中。 对于图  7(b)所示的Thomson问题，MSDA-ES和SD的数据剖面曲线相互重叠，说明它们在运行的不同阶段有各自的优势。例如，当50<�<60时五个算法中仅有SD 能够求解出约10%的测试用例，而当�>60时，MSDA-ES的数据剖面曲线总体上位于SD的数据剖面曲线的左上方，表明当允许使用更多的计算资源时，MSDA-ES能够求解出更多的问题。其原因可能是，尽管Thomson 问题具有可微的目标函数，但其并非光滑，因此 MSDA-ES  比  SD 具有性能优势。CMA-ES的整个数据剖面曲线都位于MSDA-ES的右侧，虽然其最终可以求解所有问题，但达到相同比例时大约比MSDA-ES多使用一倍的计算资源，表明其收敛速度较慢。 在图  7(c)所示的QM问题上可以观察到，SD仅需要使用极少的计算资源即可求解相当数量的问题，而MSDA-ES在获得充分计算资源时往往可以求解更多的问题。在�=20时，CMA-ES 比MSDA-ES能够求解稍多的问题，而在此之后所能求解的问题反而最少。考虑到QM问题相对简单， MSDA-ES在运行初期收敛速度较慢的原因可能在于更新协方差矩阵需要较多的计算资源，因此，提升协方差矩阵适应机制的效率可作为未来的一个研究方向。 （2）多重检验。本节使用Friedman 测试和Bonferroni后置过程对所有算法在给定测试问题的所有测试用例上进行多重检验。相比于前文中直接统计Wilcoxon秩和检验在每个测试用例上的结果，多重检验可降低族内误差，体现算法的整体性能。 表  5 给出了在三个不同测试问题上的多重检验结果。根据Friedman测试，MSDA-ES在三个不同的测试问题上均排名第一。根据Bonferroni后置过程，当置信水平为0.05时，MSDA-ES在SDR问题上显著优于SD、PSO以及DE，在Thomson问题上显著优于PSO、DE和CMA-ES。若将置信水平设置为0.1，则MSDA-ES在SDR问题上还能够显著优于CMA-ES。 4.3.6  问题规模对算法性能的影响 本节研究MSDA-ES关于流形维度的可扩展性。对于每个测试问题，将所有测试用例的流形维度从小到大排序，等分为小规模、中规模和大规模三组。在每一分组中，使用4.3.5节中使用的多重检验对所有算法的性能进行分析，所得结果呈现在表  6中。 从表  6可知，在SDR和Thomson问题上，MSDA-ES在所有分组中均取得了最好的平均排名，其他四个对比算法的平均排名结果在三个不同分组中保持一致，且与表  5中的结果基本相符。这说明在上述两个问题上，所有算法对问题规模都不敏感。产生这一结果的原因是由于上述两个问题的困难来源于目标函数本身的非光滑性而非决策空间的复杂性。因此，算法性能主要由算法处理复杂目标函数的能力决定，受到流形维度的影响较小。 相比于SDR和Thomson问题，QM问题的目标函数非常简单，且本质上等价于凸优化问题。在该问题上，算法性能受到目标函数性质的影响较小，而受到流形维度的影响相对较大，因此，在该测试问题上呈现出了与在SDR和Thomson问题上完全不同的实验结果。如表  6所示，MSDA-ES在小规模分组中排名最差，但在中规模和大规模分组表  5  不同测试问题上的多重检验结果。根据Friedman检验结果，排名第一的算法使用粗体标注。  SDR  Thomson  QM  Rank  p-Value  Rank  p-Value  Rank  p-Value MSDA-ES  10.52    14.37    29.63  SD  81.93  0.00  29.16  1.00  30.72  1.00 PSO  74.22  0.00  70.50  0.00  N/A  N/A DE  55.35  0.00  90.50  0.00  N/A  N/A CMA-ES  30.48  0.08  47.97  0.00  31.15  1.00 “Rank”表示在给定测试问题上，根据Friedman 检验得到的平均排名。 “p-Value“表示使用Bonferroni后置过程进行矫正后，MSDA-ES与对比算法在所有测试用例上性能差异的显著性水平。由于PSO和DE无法应用于QM问题，相关数据使用N/A 标准。这些符号在本文其他表格中具有相同含义。 表  6  不同规模的测试问题上的多重检验结果。根据Friedman检验结果，排名第一的算法使用粗体标注。    小规模    中规模  大规模    Rank  p-Value  Rank  p-Value  Rank  p-Value SDR MSDA-ES  10.55     10.53     10.50   SD  85.05   0.00   81.60   0.00   79.15   0.00  PSO  68.87   0.00   75.88   0.00   77.91   0.00  DE  57.57   0.00   54.03   0.01   54.44   0.01  CMA-ES  30.48   1.00   30.48   1.00   30.50   1.00  Thomson MSDA-ES  18.22     13.01     11.46   SD  26.36   1.00   29.98   1.00   31.47   1.00  PSO  70.50   0.01   70.50   0.00   70.50   0.00  DE  90.50   0.00   90.50   0.00   90.50   0.00  CMA-ES  46.91   0.64   48.51   0.22   48.58   0.27  QM MSDA-ES  30.80     29.25     28.85   SD  30.46   1.00   30.98   1.00   30.72   1.00  CMA-ES  30.25   1.00       31.27   1.00       31.94   1.00     何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  17  中排名最好，可见，MSDA-ES算法性能受到流形维度的影响，且当维度较高时具有一定的性能优势。MSDA-ES的这一性能优势主要来源于其使用的流形搜索方向适应机制。该机制仅需使用少量的几个切向量即可描述切空间中的线性依赖关系，因此可在高维流形中捕获有效的搜索方向，实现较快的搜索。 5  结论 本文将经典的CMA-ES算法从欧式空间推广至黎曼流形之中，提出了一个用于求解矩阵流形优化问题的搜索方向适应演化策略MSDA-ES。MSDA-ES使用回撤算子处理流形约束，将原始的流形约束问题转换成为一系列切空间中的无约束问题。在每个切空间中，MSDA-ES使用一个单位阵和若干个搜索方向构造协方差矩阵，并使用相应的多元高斯分布引导搜索。本文在涉及三种矩阵流形的经典问题中测试了MSDA-ES的性能，实验结果表明，MSDA-ES比已有算法在黑盒优化任务中具有显著优势。下一阶段，我们计划将MSDA-ES推广至幂流形、积流形等更为复杂的拓扑结构。同时，求解多目标矩阵流形优化也是未来研究任务之一。    致  谢 本项研究选自于国家自然科学基金项目(61773410, 61673403)  以及中国博士后科学基金项目(2019M663234)。  参 考 文 献 [1]  B.  Vandereycken,  “Low-Rank  Matrix Completion  by  Riemannian  Optimization,” SIAM  J.  Optim.,  vol.  23,  no.  2,  pp.  1214–1236, Jan. 2013. [2]  L. Cambier and P. Absil, “Robust Low-Rank Matrix  Completion  by  Riemannian Optimization,” SIAM  J.  Sci.  Comput.,  vol.  38, no. 5, pp. S440–S460, Jan. 2016. [3]  R.  H.  Keshavan,  A.  Montanari,  and  S.  Oh, “Matrix Completion From a Few Entries,” IEEE Transactions  on  Information  Theory,  vol.  56, no. 6, pp. 2980–2998, Jun. 2010. [4]  W. W. Bradbury and R. Fletcher, “New iterative methods  for  solution  of  the  eigenproblem,” Numer.  Math.,  vol.  9,  no.  3,  pp.  259–267,  Dec. 1966. [5]  Z.  Zhao,  Z.-J.  Bai,  and  X.-Q.  Jin,  “A Riemannian  Newton  Algorithm  for  Nonlinear Eigenvalue Problems,” SIAM Journal on Matrix Analysis  and  Applications,  vol.  36,  no.  2,  pp. 752–774, Jan. 2015. [6]  P.  B.  Borckmans  and  P. -a Absil, “Oriented Bounding  Box  Computation  Using  Particle Swarm  Optimization,”  in  18th  European Symposium  on  Artificial  Neural  Networks, Bruges, Belgium, 2010. [7]  L. I. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithms,” Physica  D:  Nonlinear  Phenomena,  vol.  60,  no. 1–4, pp. 259–268, Nov. 1992. [8]  C.  Yang,  J.  C.  Meza,  and  L.-W.  Wang,  “A constrained  optimization  algorithm  for  total energy  minimization  in  electronic  structure calculations,”  Journal  of  Computational Physics,  vol.  217,  no.  2,  pp.  709–721,  Sep. 2006. [9]  U.  von  Luxburg,  “A  tutorial  on  spectral clustering,” Statistics  and  Computing,  vol.  17, no. 4, pp. 395–416, Dec. 2007. [10]  A. Edelman, T. A. Arias, and S. T. Smith, “The Geometry  of  Algorithms  with  Orthogonality Constraints,” SIAM  Journal  on  Matrix  Analysis and  Applications,  vol.  20,  no.  2,  pp.  303–353, Jan. 1998. [11]  M.  Spivak, A  Comprehensive  Introduction  to Differential  Geometry,  Vol.  1,  3rd  Edition,  3rd edition. Houston, Tex: Publish or Perish, 1999. [12]  Z. Wen and W. Yin, “A feasible method for optimization  with  orthogonality  constraints,” Mathematical  Programming,  vol.  142,  no.  1–2, pp. 397–434, Dec. 2013. [13]  J.  H.  Manton,  “Optimization  algorithms exploiting  unitary  constraints,”  IEEE 18  计  算  机  学  报  2020年  Transactions  on  Signal  Processing,  vol.  50,  no. 3, pp. 635–650, Mar. 2002. [14]  P.-A.  Absil,  R.  Mahony,  and  R.  Sepulchre, Optimization  algorithms  on  matrix  manifolds. Princeton,  N.J. ;  Woodstock:  Princeton University Press, 2008. [15]  P.  Absil  and  J.  Malick,  “Projection-like Retractions  on  Matrix  Manifolds,” SIAM  J. Optim., vol. 22, no. 1, pp. 135–158, Jan. 2012. [16]  B.  Jiang  and  Y.-H.  Dai,  “A  framework  of constraint  preserving  update  schemes  for optimization  on  Stiefel  manifold,”  Math. Program.,  vol.  153,  no.  2,  pp.  535–575,  Nov. 2015. [17]  B. Gao, X. Liu, X. Chen, and Y. Yuan, “A New First-Order  Algorithmic  Framework  for Optimization  Problems  with  Orthogonality Constraints,” SIAM J. Optim., vol. 28, no. 1, pp. 302–332, Jan. 2018. [18]  W. Ring and B. Wirth, “Optimization Methods on Riemannian Manifolds and Their Application to Shape Space,” SIAM J. Optim., vol. 22, no. 2, pp. 596–627, Jan. 2012. [19]  W.  Huang,  P.-A. Absil, and K. A. Gallivan, “A Riemannian  symmetric  rank-one  trust-region method,” Math.  Program.,  vol.  150,  no.  2,  pp. 179–216, May 2015. [20]  S.  Hosseini,  W.  Huang,  and  R.  Yousefpour, “Line Search Algorithms  for  Locally  Lipschitz Functions on Riemannian Manifolds,” SIAM  J. Optim., vol. 28, no. 1, pp. 596–619, Jan. 2018. [21]  W.  Huang,  P.  Absil,  and  K.  Gallivan,  “A Riemannian  BFGS  Method  Without Differentiated  Retraction  for  Nonconvex Optimization Problems,” SIAM  J.  Optim.,  vol. 28, no. 1, pp. 470–495, Jan. 2018. [22]  P.-A.  Absil,  R.  Mahony,  and  R.  Sepulchre, “Optimization  On  Manifolds:  Methods  and Applications,”  in  Recent  Advances  in Optimization  and  its  Applications  in Engineering, 2010, pp. 125–144. [23]  W.  Huang,  “Optimization  Algorithms  on Riemannian  Manifolds  with  Applications,” Florida  State  University,  Tallahassee,  Florida, 2013. [24]  I.  Boussaïd,  J.  Lepagnot,  and  P.  Siarry,  “A survey  on  optimization  metaheuristics,” Information Sciences, vol. 237, pp. 82–117, Jul. 2013. [25]  R. Storn and K. Price, “Differential Evolution – A  Simple  and  Efficient  Heuristic  for  global Optimization over Continuous Spaces,” Journal of  Global  Optimization,  vol.  11,  no.  4,  pp. 341–359, Dec. 1997. [26]  H.-G.  Beyer  and  H.-P.  Schwefel,  “Evolution strategies–A  comprehensive  introduction,” Natural  computing,  vol.  1,  no.  1,  pp.  3–52, 2002. [27]  D.  E.  Goldberg, Genetic  algorithms  in  search, optimization,  and  machine  learning.  Reading, Mass: Addison-Wesley Pub. Co, 1989. [28]  N.  Hansen  and  A.  Ostermeier,  “Completely Derandomized  Self-Adaptation  in  Evolution Strategies,” Evolutionary  Computation,  vol.  9, no. 2, pp. 159–195, Jun. 2001. [29]  A. Banks, J. Vincent, and C. Anyakoha, “A review  of  particle  swarm  optimization.  Part  I: background  and  development,”  Natural Computing, vol. 6, no. 4, pp. 467–484, 2007. [30]  P.  B.  Borckmans,  M.  Ishteva,  and  P.-A.  Absil, “A  Modified  Particle  Swarm  Optimization Algorithm  for  the  Best  Low  Multilinear  Rank Approximation  of  Higher-order  Tensors,”  in Proceedings of the 7th International Conference on  Swarm  Intelligence,  Berlin,  Heidelberg, 2010, pp. 13–23. [31]  S.  Colutto,  F.  Fruhauf,  M.  Fuchs,  and  O. Scherzer,  “The  CMA-ES  on  Riemannian Manifolds  to  Reconstruct  Shapes  in  3-D  Voxel Images,”  IEEE  Transactions  on  Evolutionary Computation,  vol.  14,  no.  2,  pp.  227–245,  Apr. 2010. [32]  D.  V.  Arnold,  “On  the  Use  of  Evolution Strategies  for  Optimization  on  Spherical Manifolds,” in Parallel  Problem  Solving  from Nature – PPSN XIII, 2014, pp. 882–891. [33]  D.  V.  Arnold  and  A.  Lu,  “An  evolutionary algorithm  for  depth  image  based  camera  pose   何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  19  estimation  in  indoor  environments,”  in 2016 IEEE  Congress  on  Evolutionary  Computation (CEC),  Vancouver,  BC,  Canada,  2016,  pp. 3871–3878. [34]  E.  Stiefel,  “Richtungsfelder  und Fernparallelismus  in  n-dimensionalen Mannigfaltigkeiten,” ETH Zurich, 1935. [35]  N.  Boumal,  B.  Mishra,  P.  A.  Absil,  and  R. Sepulchre,  “Manopt,  a  matlab  toolbox  for optimization on manifolds,” Journal of Machine Learning  Research,  vol.  15,  no.  1,  pp. 1455–1459, 2013. [36]  N.  Hansen,  A.  S.  P.  Niederberger,  L.  Guzzella, and P. Koumoutsakos, “A Method for Handling Uncertainty  in  Evolutionary  Optimization  With an  Application  to  Feedback  Control  of Combustion,”  IEEE  Transactions  on Evolutionary  Computation,  vol.  13,  no. 1,  pp. 180–197, Feb. 2009. [37]  H.-G. Beyer and S. Finck, “On the Design of Constraint  Covariance  Matrix  Self-Adaptation Evolution  Strategies  Including  a  Cardinality Constraint,” IEEE Transactions on Evolutionary Computation, vol. 16, no. 4, pp. 578–596, Aug. 2012. [38]  T. Suttorp, N. Hansen, and C. Igel, “Efficient covariance  matrix  update  for  variable  metric evolution strategies,” Machine  Learning,  vol. 75, no. 2, pp. 167–197, May 2009. [39]  J.  Schäfer  and  K.  Strimmer,  “A  Shrinkage Approach  to  Large-Scale  Covariance  Matrix Estimation  and  Implications  for  Functional Genomics,” Statistical  Applications  in  Genetics and Molecular Biology, vol. 4, no. 1, Jan. 2005. [40]  I.  Loshchilov, “LM-CMA:  An  Alternative  to L-BFGS  for  Large-Scale  Black  Box Optimization,” Evolutionary  Computation,  vol. 25, no. 1, pp. 143–171, Mar. 2017. [41]  Nikolaus Hansen, Sibylle D. Müller, and Petros Koumoutsakos, “Reducing the Time Complexity of  the  Derandomized  Evolution  Strategy  with Covariance  Matrix  Adaptation  (CMA-ES),” Evolutionary  Computation,  vol.  11,  no.  1,  pp. 1–18, Mar. 2003. [42]  H.-G.  Beyer,  “Convergence  Analysis  of Evolutionary Algorithms That Are Based on the Paradigm  of  Information  Geometry,” Evolutionary  Computation,  vol.  22,  no.  4,  pp. 679–709, Dec. 2014. [43]  D.  S.  Broomhead  and  M.  J.  Kirby, “Dimensionality Reduction Using Secant-Based Projection  Methods:  The  Induced  Dynamics  in Projected Systems,” Nonlinear Dyn, vol. 41, no. 1–3, pp. 47–67, Aug. 2005. [44]  P.-A. Absil and S. Hosseini, “A collection of nonsmooth Riemannian optimization problems,” ICTEAM  Institute,  UCLouvain,  B-1348 Louvain-la-Neuve,  Belgium, UCL-INMA-2017.08-v1, Sep. 2017. [45]  J. J. Thomson, “XXIV. On the structure of the atom:  an  investigation  of  the  stability  and periods of oscillation of a number of corpuscles arranged  at  equal  intervals  around  the circumference of a circle; with application of the results to the theory of atomic structure,” The London,  Edinburgh,  and  Dublin  Philosophical Magazine and Journal of Science, vol. 7, no. 39, pp. 237–265, Mar. 1904. [46]  E.  Polak,  Optimization:  Algorithms  and Consistent  Approximations,  1997  edition.  New York: Springer, 1997. [47]  S.  Das  and  P.  N.  Suganthan,  “Differential Evolution:  A  Survey  of  the  State-of-the-Art,” IEEE  Transactions  on  Evolutionary Computation,  vol.  15,  no.  1,  pp.  4–31,  Feb. 2011. [48]  I. Loshchilov, T. Glasmachers, and H. G. Beyer, “Large  Scale  Black-box  Optimization  by Limited-Memory  Matrix  Adaptation,” IEEE Transactions on Evolutionary Computation, vol. 23, no. 2, pp. 353–358, Apr. 2019. [49]  F. van den Bergh and A. P. Engelbrecht, “A Cooperative  Approach  to  Particle  Swarm Optimization,” Trans. Evol. Comp, vol. 8, no. 3, pp. 225–239, Jun. 2004. [50]  M. Sepesy Maučec and J. Brest, “A review of the  recent  use  of  Differential  Evolution  for Large-Scale  Global  Optimization:  An  analysis 20  计  算  机  学  报  2020年  of  selected  algorithms  on  the  CEC  2013  LSGO benchmark  suite,”  Swarm  and  Evolutionary Computation, Aug. 2018. [51]  J. Vesterstrom and R. Thomsen, “A comparative study  of  differential  evolution,  particle  swarm optimization,  and  evolutionary  algorithms  on numerical  benchmark  problems,”  2004,  pp. 1980–1987. [52]  X. Li and X. Yao, “Tackling high dimensional nonseparable  optimization  problems  by cooperatively  coevolving  particle  swarms,” 2009, pp. 1546–1553. [53]  X. Li and X. Yao, “Cooperatively Coevolving Particle Swarms for Large Scale Optimization,” IEEE  Transactions  on  Evolutionary Computation,  vol.  16,  no.  2,  pp.  210–224,  Apr. 2012. [54]  D. Molina, A. LaTorre, and F. Herrera, “An Insight  into  Bio-inspired  and  Evolutionary Algorithms  for  Global  Optimization:  Review, Analysis,  and  Lessons  Learnt  over  a  Decade  of Competitions,” Cogn.  Comput.,  vol.  10,  no.  4, pp. 517–544, Aug. 2018.       何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  21  附录1. 命题1的证明 假设   = () ，使用归纳法证明  +1 =  +1  。  (1)首先，证明  +1 = +1 。 根据假设 ()=(), � ()=�()，有   ()=σ ()  1−� ()  () + �   ()� () =1  =σ()  1−�() () + � ()�()=1  =σ()  1−�() () + � ()�()=1  =() 其中，上式第二行利用了Grassmann流形上切空间投影的性质() =()。  根据回撤算子的定义，可知 ()   () =() () =()  () ，进而有   ()   ()  = ()  ()  = ()  ()   因此，�1和�2对应的候选矩阵具有相同的列空间，被映射至流形上的同一位置，具有相同的目标函数值。  由于�1和�2的群体具有相同的目标值和排序结果，因此在重组时，有   = � :()=1= �:()=1= ,   +1 = ()   =()   = +1   (2)其次，证明  +1 = +1   ( =1, …, )。 在更新第一个搜索方向时，    1= μ  ()= μ σ()=1，可知�1和�2中成功变异方向具有相同列空间，因此后续更新得到的搜索方向同样具有相同列空间。具体而言，对�2中的第一个搜索方向更新可得，  1′= 1−  1()+  2−  1 = 1− 1()+  2− 1=1′ 进而有  1=  1,  1′  ()  1′,  1′  ()= 1, 1′  () 1′, 1′  ()= 1, 1′  () 1′, 1′  ()= 1, 1′ () 1′, 1′ ()=1 即，成功变异方向在更新后的第一个搜索方向上的相对投影长度相同。同时，由  2=1 1+ 12  1− 1 1 +1  =2 可知用于更新第二个搜索方向的成功变异方向同样具有相同的列空间，同理可得，  ′=′,  =1, …,  因此，对 =1, …, ，有   +1 =  +1   ′ = +1  ′ = +1  ′ = +1 。  (3)最后，由于�1和�2中的群体具有相同的目标函数值，因此连续两代群体目标函数值的秩和差值相同。由于变异强度适应机制仅依赖于该秩和差值，因此有   +1 = +1 ,   +1 = +1 。 证毕. 22  计  算  机  学  报  2020年  附录2. 用于矩阵流形优化的DE算法                        附录3. 用于矩阵流形优化的CMA-ES算法     算法II.    DE 输入：ℳ：矩阵流形；:×→：目标函数；�:ℳ→ℳ：指数映射；:ℳ→ℳ：对数映射；：群体规模；：差分向量缩放因子；：算术交叉系数；：执行变异的概率。 输出：已找到的最优解。 1: FOR ( ←1:)  ←ℳ中随机分布的矩阵 2: WHILE(未满足终止条件){ 3:    FOR ( ←1:) { 4:    1, 2, 3←{1,2, …, }中不同的三个随机数 5:  IF(<) 6:      = 1 +  2 − 3   7:  ELSE 8:      = 1 +  2 + 3 −                        2 1     9:  IF((�  )<())    ←�   10:  } 11: } 指 0,1 上均匀分布的随机数。 算法III.    CMA-ES 输入：ℳ：矩阵流形；:×→：目标函数；:ℳ→ℳ：回撤算子；:×→ℳ：切空间投影。 输出：已找到的最优解。 1:   ←0   2:   ←ℳ中随机分布的矩阵,  σ  ←1 3:   ←0×,    ←0×,   ∈{1, ⋯, } 4: WHILE(未满足终止条件){ 5:  FOR ( ←1:λ) { 6:       ←     7:       ← 8:      FOR (←1:)  9:           ← 1−, +,       10:  } 11:       ={      1  , ⋯,       λ  }   12:  将  按升序排列 13:  �←   �:=1   14:     +1 ←     �  :=1    15:  (+1)= 1− ()+  2− �   16:  FOR (←1:m) 17:     +1 = +1   1−,   +                                                  , 2−, �  18: σ +1 ←σ  � 2 | +1  2ℳ−1     19:  ←+1   20: }    何笑雨、周育人、陈泽丰：基于演化策略的矩阵流形黑盒优化方法  23     He  Xiao-Yu,  born  in  1988,  Ph.D., postdoctoral  fellow.  His  research interests  focus  on  evolutionary computation.     Zhou  Yu-Ren,  born  in  1965,  Ph.D.,  professor.  His research interests focus on evolutionary computation.    Chen  Ze-Feng,  born  in  1990,  Ph.D., postdoctoral fellow.  His research  interests  focus  on  evolutionary computation.                   Background This  work  was  supported  by  National  Natural Science  Foundation  of  China  (61773410,  61673403) and  China  Postdoctoral  Science  Foundation (2019M663234).   This  work  focuses  on  the  manifold  constrained matrix  optimization  problems  which  abound  in  data mining,  numerical  linear  algebra,  signal  processing, and  electronic  design.  Most  existing  approaches are white-box,  thereby  the  efficiency  and  applicability strongly  dependent  on the  problem structure  and the inherent  mathematic  properties.  The  algorithm presented  here,  to  the  best  of  our  knowledge,  is  the first retraction-based evolutionary manifold optimizer. This algorithm is a general-purpose solver that can be applied  to  most  kinds  of  manifolds whenever  the retraction  operator  and  tangent  space  projection  are well  defined. This  work  gives  some  inspirations  on black-box approaches to complicated matrix manifold problems  which  may  not  be  handled  easily  with traditional algorithms.