次模优化初探

前置概念

次模函数

次模性（Submodularity，又称子模、亚模等）是集合函数的一个属性，具有深层的理论意义和深远的应用前景。拥有着次模性的函数称为次模函数。

给定有限集$V$，通常称作基础集（ground set）。次模函数有两个等价的定义，第一个定义依赖于边际增益（marginal gain）的概念。

定义$f(\varnothing )=0$，即集合为空集时值为0。

定义 1（边际增益）对于集合函数$f:2^V\to \mathbb{R}$，子集$S\subseteq V$，集合中的元素$e\in V$，记 Δ f ( e ∣ S ) : = f ( S ∪ { e } ) − f ( S ) \Delta_f(e|S):=f(S\cup \{e\})-f(S) Δf​(e∣S):=f(S∪{e})−f(S)为元素$e$添加到集合$S$中产生的边际增益。由于$f$在上下文中是清晰的，故这里省去，记作$\Delta (e|S)$。

定义2（次模性）一个函数$f:2^V\to \mathbb{R}$被称作是次模的，当且仅当对于任意两个子集$A\subseteq B\subseteq V$，元素$e\in V\setminus B$，满足$\Delta (e|A)\ge \Delta (e|B)$。等价地，对于每两个子集$A,B\subseteq V$，满足$f(A\cap B)+f(A\cup B)\le f(A)+f(B)$时也满足次模性。

次模函数展示了一个自然的边际效用递减（diminishing returns）的性质。

在下面的例子¹中，$V$指的是传感器能放置的位置集合，$f(S)$指的是将传感器放置在位置$S$获得的效用（比如检测性能）。蓝色区域是已经存在的传感器（$S_1、S_2、S_3、S_4$）能探测的区域，红色区域是新增传感器（$S^{\prime }$）能探测的区域，如果一开始有更多的传感器被放置（$S_3$、$S_4$），那么新增传感器能增加探测的区域相对来说便小了，即 Δ ( S ′ ∣ { S 1 , S 2 } ) ≥ Δ ( S ′ ∣ { S 1 , S 2 , S 3 , S 4 } ) \Delta(S'|\{S1,S2\})\ge\Delta(S'|\{S1,S2,S3,S4\}) Δ(S′∣{S1,S2})≥Δ(S′∣{S1,S2,S3,S4})。

次模函数有很重要的一类，即单调次模函数，增大输入集不会使得函数值减小。

定义3（单调性）函数$f:2^V\to \mathbb{R}$被称作是单调的，当且仅当对于任意两个子集$A\subseteq B\subseteq V$，满足$f(A)\le f(B)$。

注意到$f$是单调函数时，其所有的边际增益都是非负的，即对于任意集合$A\subseteq V$，元素$e\in V$，满足$\Delta (e|A)\ge 0$。

次模优化

次模函数优化简称为次模优化。近年来，次模优化开始在机器学习和数据挖掘方面得到应用，包括自动摘要、多文档摘要、自动提取、主动学习、传感器放置等多个方面。次模优化相关问题主要有次模函数最大化、次模函数最小化等问题。这里主要讨论次模函数最大化，拥有以下形式：
$$\underset{S\subseteq V}{\max }f(S)subjecttosomeconstraintsonS.\text{\hspace{1em}(1)}$$
其中最简单的约束是基数约束（cardinality constraints），即$|S|\le k$，选取个数不超过$k$的子集使得次模函数的值尽可能大。其他的还有背包约束、拟阵约束等。对于无约束次模函数最大化问题，如果函数是单调的，该问题可以在多项式时间内精确求解，然而，大多数次模函数最大化问题是NP难的，主要都是用近似算法来求解。

流次模优化

在流模型中，数据以流（streaming）的形式呈现，算法只能“阅读”数据，即数据只有被访问到才能利用它们的信息。同时假设存储资源是有限的（远远小于数据规模），单位时间内，只能看到当前数据，并且要在下一个数据到达之前决定是否存储该数据点。该决定通常是不可更改的，即一旦该数据被丢弃，之后就无法再利用它的任何信息，目的是从数据流中选择满足某些性质的子集，使得其目标函数收益尽可能大。如果收益函数是次模函数，则称该模型为流次模函数最大化问题。

数据摘要

Data summarization（数据摘要，或许也可翻译成数据总结）已经成为甚至TB级别数据的有价值重大的工具。由于其引人注目的理论属性，次模函数成为了摘要算法的焦点。

简单的说，数据摘要就是从大量数据中选取有代表性的元素，即根据某些目标函数从数据集中提取子集。由于许多“代表性”的概念满足次模性，即收益递减的直观概念。因此上述问题能规约为最大化受基数约束的次模函数问题。

RBF核函数

在机器学习中，radial basis function（径向基函数）核，是一种常用的核函数。它是支持向量机分类中最为常用的核函数。关于两个样本$x$和$x^{\prime }$ 的RBF核可表示为某个“输入空间”的特征向量，它的定义如下所示：
$$k(x,x^{\prime })=\exp (-\frac{\Vert x-x^{\prime }{\Vert }_2^2}{2{\sigma }^2})\text{\hspace{1em}(2)}$$
$||x-x^{\prime }|{|}_2^2$可以看做两个特征向量间的平方欧式距离。$\sigma$是一个自由参数。一种等价但更为简单的定义是设一个新的参数$\gamma =\frac{1}{2{\sigma }^2}$，那么
$$k(x,x^{\prime })=exp(-\gamma ||x-x^{\prime }|{|}_2^2)\text{\hspace{1em}(3)}$$
由于RBF核的值随距离而减小，并且范围在0（极限）和1（当x = x’时）之间，因此它可以一种好的相似性度量表示法。

通过核函数可以构造核矩阵${\Sigma }_S=[k(e_i,e_j){]}_{ij}$，包含$S$中的所有相似度点对。

该矩阵在对角线上具有最大值，因为元素与自己本身是最相似的，而在非对角线上的值表示不同元素之间的相似性，因此通常较小。由于我们寻求全面的总结，我们更感兴趣的是0附近的非对角线上的值对。这种直觉已经在信息向量机（IVM）的上下文中得到了正式处理。

IVM

IVM（Informative Vector Machine）²是一个高斯过程（Gaussian Process，GP），贪心地选择数据点的子集并跟踪GPs的后部分布。基于多样性的论点，作者建议迭代选择能最好地涵盖了训练数据的点。这种直觉可以形式化为最大化核矩阵的对数行列式：
$$f(S)=\frac{1}{2}logdet(\mathcal{I}+\alpha {\Sigma }_S)\text{\hspace{1em}(4)}$$
其中$\alpha \in {\mathbb{R}}_{+}$是数值鲁棒性的缩放参数，$\mathcal{I}$是单位矩阵。

在（Seeger,2004）³中，该函数被证明是单调次模函数。它的函数值不依赖于数据集$V$，而仅取决于核函数$k$和摘要大小$K$。这使得其成为流模式中数据摘要的理想函数。

贪心算法求解

问题描述

在数据集中进行数据摘要，归约成了次模优化问题，即主要解决的问题是最大化下列函数，一个受基数约束的单调次模函数：
$$\underset{S\subseteq V}{\max }f(S)s.t.|S|\le k\text{\hspace{1em}(5)}$$
这里的目标函数采用IVM。由于贪心算法在流式数据中表现并不好，所以这里先不采用流式数据。

算法原理

由于该问题是NP难的，只能用近似算法求解，Nemhauser等⁴证明了简单的贪心算法是很有效的。从空集$S_0=\varnothing$开始，每一次在整个数据集的迭代过程$i$中，选择元素$e\in V$最大化(5)式，即：
S i = S i − 1 ∪ { a r g   max ⁡ e ∈ V Δ f ( e ∣ S i − 1 ) } (6) S_i = S_{i-1}\cup \{arg\ \max_{e\in V} \Delta_f(e|S_{i-1}) \} \tag{6} Si​=Si−1​∪{arg e∈Vmax​Δf​(e∣Si−1​)}(6)
用$S^g$表示贪心算法所求解，Nemhauser等证明了$f(S^g)\ge (1-1/e)OPT$，其中$e$是自然对数的底数，$OPT$是理想最优解，即贪心算法有着$(1-1/e)\approx 0.63$的近似比，对于若干种单调次模函数，$(1-1/e)$是最好的近似率了。

代码实现

代码链接 ：https://github.com/wangzhebufangqi/SubmodularOptimization

环境：Visual Studio 2019；C++17

其中

• DataTypeHandling.h文件中定义了一些数据类型
• Matrix.h文件中实现了一些对对称正定矩阵的操作如秩一校正、cholesky分解等
• Kernel.h中定义了Kernel类，表示核函数
• RBFKernel.h中定义了继承自Kernel类的RBFKernel类，表示RBF核函数
• SubmodularFunction.h中定义了SubmodularFunction类，表示次模函数
• IVM.h中定义了继承自SubmodularFunction类的IVM类，表示具体的次模函数IVM
• FastIVM.h中定义了继承自IVM类的FastIVM类，比起IVM类计算更快
• SubmodularOptimizer.h中定义了SubmodularOptimizer类，包含一个次模函数f、基数约束K等，表示一个次模优化器，解决一个次模优化问题
• Greedy.h中定义了继承自SubmodularOptimizer类的Greedy类，表示用贪心算法求解次模优化问题

主要代码

以下代码为Greedy.h文件中，fit函数部分。功能为在数据集上迭代K次，每次挑选具有最大边际增益的元素。

/*
*
 @brief 在整个数据集中挑选拥有最大边际增益的元素。一直重复直到选择完K个元素。
*        调用'get_solution'得到结果。
* @param X 整个数据集的常引用。
* @param iterations：其实没什么用，贪心算法在任何情况下在整个数据集上迭代K次。
*/
void fit(vector<vector<data_t>> const& X, vector<idx_t> const& ids,
    unsigned int iterations = 1) {
    vector<unsigned int> remaining(X.size());//数据集中剩下未被选择的元素序号
    iota(remaining.begin(), remaining.end(), 0);//0,1,2，...,X.size()-1
    data_t fcur = 0;

    while (solution.size() < K && remaining.size() > 0) {//K个元素未选择完、数据集中还剩余有元素
        vector<data_t> fvals;//每个剩余未被选择元素被添加进解后的函数值，每次选取一个最大值
        fvals.reserve(remaining.size());//共remaining.size()个，每次循环值会减小

        /*
        * 贪心算法挑选具有最大边际增益的元素。等价于挑选导致最大函数值的元素。
        * 没必要再显式地再计算这个增长值。
        * 由于pos = solution.size()，即pos>=solution.size()，
        * 故ftmp就是假设将X[i]添加至当前解后的函数值。保存至数组fvals中。
        */
        for (auto i : remaining) {
            data_t ftmp = f->peek(solution, X[i], solution.size());
            fvals.push_back(ftmp);
        }

        /*
        * max_ele为fvals中拥有最大函数值的元素的下标，范围为[0,remaining.size()-1]，（remaining.size()==fvals.size()）
        * fcur为当前的最大函数值，这个函数值导致了最大边际增益，即fcur-fcur(上一轮)=最大边际增益
        * max_idx为X中的元素序号，范围为[0,X.size()-1]，
        */
        unsigned int max_ele = distance(fvals.begin(), max_element(fvals.begin(), fvals.end()));
        fcur = fvals[max_ele];
        unsigned int max_idx = remaining[max_ele];

        /*
        * 得出了最大值，贪心选择该元素。故将序号为max_idx的元素添加进当前解，
        * 同时，次模函数进行更新。
        */
        f->update(solution, X[max_idx], solution.size());
        solution.push_back(X[max_idx]);

        /*
        * this->ids序列保存每一次被选择的元素序号
        */
        if (ids.size() >= max_idx) {
            this->ids.push_back(max_idx);
        }
        /*
        * 挑选完一个元素，将其从remaining中移去
        */
        remaining.erase(remaining.begin() + max_ele);

    }

    fval = fcur;//最后一步的函数值
    is_fitted = true;
}

结果测试

小规模测试

从含有8个元素的集合data中根据IVM从中选取4个最具代表性的元素。

#include<iostream>
#include<vector>
#include<numeric>

#include"DataTypeHandling.h"
#include"Greedy.h"
#include"FastIVM.h"
#include"RBFKernel.h"

int main() {
    vector<vector<data_t>> data = {
        {0.0, 0.0},
        {1.0, 1.0},
        {0.5, 1.0},
        {1.0, 0.5},

        {0.0, 0.5},
        {0.0, 1.5},
        {0.0, 1.0},
        {0.5, 0.5}
    };

    vector<idx_t> ids;
    ids.resize(data.size());
    iota(ids.begin(), ids.end(), 0);

    unsigned int K = 4;

    FastIVM fastIVM(K, RBFKernel(), 1.0);
    Greedy greedy(K, fastIVM);
    greedy.fit(data, ids);
    vector<vector<data_t>> solution = greedy.get_solution();
    vector<idx_t> solution_ids = greedy.get_ids();
    double fval = greedy.get_fval();

    cout << "Found a solution with fval = " << fval << endl;
    for (auto x : solution) {
        for (auto xi : x) {
            cout << xi << " ";
        }
        cout << endl;
    }

    cout << "ids = {";
    for (auto i : solution_ids) {
        cout << i << " ";
    }
    cout << "}" << endl;
    return 0;

}

运行结果：

KDDCup99数据集

数据集下载地址：https://www.dbs.ifi.lmu.de/research/outlier-evaluation/input/KDDCup99.tar.gz

以下代码读取数据后，从48113个元素中挑选出了50个最具代表性的元素。

#include<iostream>
#include<vector>
#include<numeric>
#include<fstream>
#include<sstream>
#include<tuple>

#include"DataTypeHandling.h"
#include"Greedy.h"
#include"FastIVM.h"
#include"RBFKernel.h"

vector<vector<data_t>> read_arff(string const& path) {
    vector<vector<data_t>> X;

    string line;
    ifstream file(path);

    if (file.is_open()) {
        while (getline(file, line)) {
            // Skip every meta information
            if (line.size() > 0 && line[0] != '@' && line != "\r") {
                vector<data_t> x;
                stringstream ss(line);
                string entry;
                // All entries are float, but the last one (the label, string) and the second to last(the id, integer). Skip both.
                while (getline(ss, entry, ',') && x.size() < 78) {
                    if (entry.size() > 0) { //&& entry[0] != '\''
                        x.push_back(static_cast<float>(atof(entry.c_str())));
                    }
                }
                if (X.size() > 0 && x.size() != X[0].size()) {
                    cout << "Size mismatch detected. Ignoring line." << std::endl;
                }
                else {
                    X.push_back(x);
                }
            }
        }
        file.close();
    }

    return X;
}

auto evaluate_optimizer(SubmodularOptimizer& opt, vector<vector<data_t>>& X, vector<idx_t> ids) {
    auto start = chrono::steady_clock::now();
    opt.fit(X,ids);
    auto end = chrono::steady_clock::now();
    chrono::duration<double> runtime_seconds = end - start;
    auto fval = opt.get_fval();
    cout << "Selected " << opt.get_solution().size() << endl;
    //最终的函数值；运行时间；存储元素个数；候选解集的个数
    return make_tuple(fval, runtime_seconds.count(), opt.get_num_elements_stored(), opt.get_num_candidate_solutions());
}

int main() {
    cout << "Reading data" << endl;
    auto data = read_arff("./KDDCup99/KDDCup99_withoutdupl_norm_1ofn.arff");
    cout << data.size() << endl;

    unsigned int K = 50;

    vector<idx_t> ids;
    ids.resize(data.size());
    iota(ids.begin(), ids.end(), 0);

    FastIVM fastIVM(K, RBFKernel(sqrt(data[0].size()), 1.0), 1.0);

    tuple<data_t, double, unsigned long, unsigned int> res;

    cout << "Selecting " << K << " representatives via fast IVM with Greedy" << endl;
    Greedy fastGreedy(K, fastIVM);
    res = evaluate_optimizer(fastGreedy, data, ids);
    cout << "\t fval:\t\t" << get<0>(res) << "\n\t runtime:\t" << get<1>(res) << "s\n\t memory:\t" << get<2>(res) << "\n\t num_sieves:\t" << get<3>(res) << "\n\n" << endl;

    cout << endl;
    vector<idx_t> solution_ids = fastGreedy.get_ids();
    for (auto x : solution_ids) {
        cout << x << ' ';
    }
    cout << endl;
}

运行结果：

后记

根据(Badanidiyuru et al.，2014)⁵做毕设，结合(Buschjäger et al.,2020)⁶、次模函数最大化的流算法综述⁷等文献，总算对次模优化有了一点点了解，在此开始记录学习过程。知识浅薄，错漏之处在所难免，欢迎指正。
——————————
2021-5-9 更新
大无语事件，论文查重会检查到这篇，“我抄我自己”
毕业后再恢复。
——————————
2021-6-5更新
恢复文章。原先想在这篇介绍完贪心算法后再另开文章介绍SieveStreaming算法和SieveStreaming++算法的，但发现自己还是不敢说已经完全理解透了上述算法，所以就先搁置起来。（最后的毕设分数也是不尽如人意，只能怪自己实验确实做的一般般）。

参考资料

---

1. A. Krause and D. Golovin. Submodular function maximization. In Tractability: Practical Approaches to Hard Problems. Cambridge University Press, 2013. ↩︎

2. Lawrence N, Seeger M, Herbrich R. Fast sparse Gaussian process methods: The informative vector machine[C]//Proceedings of the 16th annual conference on neural information processing systems. 2003 (CONF): 609-616. ↩︎

3. Seeger, M. Greedy forward selection in the informative vector machine. Technical report, Technical report, University of California at Berkeley, 2004. ↩︎

4. G. L. Nemhauser, L. A. Wolsey, and M. L. Fisher. An analysis of approximations for maximizing submodular set functions - I. Mathematical Programming, 1978. ↩︎

5. Badanidiyuru A, Mirzasoleiman B, Karbasi A, et al. Streaming submodular maximization: Massive data summarization on the fly[C]//Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining. 2014: 671-680. ↩︎

6. Buschjäger S, Honysz P J, Morik K. Very Fast Streaming Submodular Function Maximization[J]. arXiv preprint arXiv:2010.10059, 2020. ↩︎

7. 杨瑞琪, 徐大川, 杜东雷,等. 次模函数最大化的流算法综述[J]. 运筹学学报, 2020(2):73-86 ↩︎