程序设计实习2023复习

程序设计实习2023复习

这个东西也许和 OI 有关?反正也要复习,干脆写篇博客。

introduction

传统算法设计一般考虑有效算法(多项式时间),以及追求“精确解”,尤其是“经典”算法或者数据结构。这节课一般考虑那些不那么传统的、更加偏向现代的算法设计。

一种就是近似算法,比如探索 NP-hard 问题在多项式时间内可以近似到多么精确,或者是本来存在多项式时间算法的问题,能否使用近似解改良复杂度。

对于最小化问题,一般使用近似比来衡量一个解的质量,如果对于某个输入 \(I\),你得到的解为 \(\mathrm{ALG}(I)\) 而最优解为 \(\mathrm{OPT}(I)\),那么我们可以通过 \(\max_{I}\frac{\mathrm{ALG}(I)}{\mathrm{OPT}(I)}\) 来衡量你给出的解的质量,这个值的取值范围显然是 \([1,+\infty)\),越小说明质量越好。

称该算法是 \(\alpha-\) 近似的,当对任意 \(I\),都有 \(\mathrm{ALG}(I)\le \alpha \cdot \mathrm{OPT}(I)\)

例如 TSP 问题(旅行商问题,就是经过所有点并且回到出发地的最短长度),找到最小生成树上的最优解就是一个 \(2-\) 近似算法;找平面点集直径随便选择一个点然后找到离这个点最远的也是一个 \(2-\) 近似算法。

另一种意义上的近似是算法本身具有随机性,对于确定的输入得到随机的输出,比如有一定概率出错,但是不出错时有很好的性能保证。

近似 + 随机是现代算法设计的大趋势。

random

随机数

要在程序中生成随机数,可以考虑 C++ 中的 random 库,预设的几个标准随机数生成器如基于线性同余的 minstd_rand0minstd_randknuth_b,基于 Mersenne Twister 的 mt19937 ,还有基于 Lagged Fibonacci 的 ranlux24ranlux48,一般使用 minstd_rand 就行了。

然后还有从很多不同的特殊分布中采样,例如均匀独立采样 uniform_int_distributionuniform_real_distribution 分别是在一个区间里面均匀采整数或者实数,其他例如两点分布的 bernoulli_distribution\(p\) 的概率取 \(1\) 其余取 \(0\),还有二项分布 binomial_distribution ,另外一个较常用到的是正态分布(高斯分布) normal_distribution

例如下面程序可以生成 \(10\) 个从标准正态分布 \(\mathcal{N}(0,1)\) 中的采样。

#include<random>
using namespace std;
int main(){
    std::minstd_rand gen(114514);
    std::normal_distribution<double> dis(0.0,1.0);
    for(int i=1;i<=10;++i)
        printf("%.5lf\n",dis(gen));
    return 0;
}

而更加一般的离散分布采样就是 discrete_distribution ,例如下图:

#include<random>
using namespace std;
int main(){
    std::minstd_rand gen(114514);
    std::discrete_distribution<> dis({1,4,6,4,1});
    for(int i=1;i<=10;++i)
        printf("%d\n",dis(gen));
    return 0;
}

这里传入的参数是一个长度为 \(n\) 的数组 \(a\),下标从 \(0\)\(n-1\),其中 \(i\) 数字生成的概率为 \(\frac{a_i}{\sum_{j=0}^{n-1}a_j}\)

随机算法

考虑一个随机算法,如果有 \(\delta\) 的出错概率,也就是说正确概率为 \(1-\delta\),如果算法只需要运行一次 \(\delta\) 取足够小的常数就行了,如果要多次运行且每次成功,可能就需要让 \(\delta\) 与运行次数相关。如果运行 \(q\) 次且每次失败概率为 \(\delta\),根据 Union Bound \(\mathrm{Pr}[\bigcup_i X_i]\le \sum_i \mathrm{Pr}[X_i]\),那么存在一次失败概率至多 \(q\delta\)

算法具体分析的时候也许需要用到各种理论推导得出需要运行的次数,但是实际上只要试试够用就行了。

一些典型随机算法

测试矩阵相乘结果是否正确

给定三个 \(n\times n\) 的矩阵 \(A,B,C\),问是否有 \(AB=C\)

考虑生成一个长度为 \(n\) 的随机向量 \(x\),每个值从集合 \(\{0,1\}\) 中等概率选取,如果 \(AB=C\) 那么必须有 \(ABx=Cx\),如果 \(AB\ne C\),那么考虑 \(AB\)\(C\) 中至少有一列不完全相同,不妨假设是第 \(t\) 列,先确定 \(x\) 其他位置的值,然后考虑 \(t\) 位置的值,容易发现 \(t=0\)\(t=1\) 中至少有一个会导致 \(ABx\ne Cx\),所以可以看出此时 \(ABx\ne Cx\) 的概率不低于 \(0.5\)

重复随机 \(T\) 次,那么容易发现出错概率不超过 \(0.5^T\),更一般地,如果有一次出错概率 \(p\le 0.5\),那么 \(T\) 次都出错的概率就是 \(p^T\),要有 \(\delta\) 的出错概率或者说 \(1-\delta\) 的成功率就只需要有 \(T=O(\log (1/\delta))\) 即可。

\(\delta\) 取常数的话这个算法便可以做到 \(O(n^2)\) 判断。

Median Trick

现在有⼀个⿊盒能以 \(p>0.5\) 的概率正确回答某个 Yes/No 问题的答案,能否将其强化成任意的 \(1-\delta\)

同样的,重复随机 \(T\) 次,取出现次数较多的作为答案,要选多大才行?

Chernoff Bound:设独立随机变量 \(X_1,\dots,X_n\in [0,1]\),令 \(X=\sum_{i=1}^nX_i,\mu=\mathrm{E}[X]\),则对 \(t\in [0,1]\)\(\mathrm{Pr}[|X-\mu|\ge t\mu]\le 2\exp(-t^2\mu/3)\)

随机 \(T\) 次回答正确次数期望时 \(pT\),所以也就要使得 \(\mathrm{Pr}[|X-pT|\ge (p-0.5)T]\le \delta\),视 \(p\) 为常数容易发现 \(T\) 取到 \(O(\log (1/\delta))\) 即可。

Hoeffding inequality:设独立随机变量 \(X_1,\dots,X_n\in [a,b]\),令 \(X=\sum_{i=1}^nX_i,\mu=\mathrm{E}[X]\),则对任意 \(t\)\(\mathrm{Pr}[X-\mu\ge t]\le 2\exp\left(-\frac{2t^2}{n(b-a)^2}\right)\)

recursion(递归)

这里是介绍经典算法,没啥好说的()

hash

哈希算法

\(h:[n]\to[m]\) 是随机哈希,指对于每个 \(i\in [n]\),有 \(h(i)\)\([m]\) 上的均匀随机分布。一般要求是对于 \(i\ne j\),有 \(\mathrm{Pr}[h(i)=h(j)]\le 1/m\),其中每个 \(j\in [m]\) 称为 bucket。

考虑分布式环境的负载均衡问题,有 \(n\) 个请求要丢给 \(m\) 台服务器处理,如何使得每台服务器处理的请求个数尽量平均?可以使用上面讲到的随机哈希。

那么如果需要新增服务器呢?重构哈希函数太麻烦。

可以使用一致性哈希算法:构造一个足够大的哈希函数 \(h:\mathrm{int} \to \mathrm{int}\),然后把请求和服务器都哈希一下到一个环上,然后每个请求丢给这个环上对应后面第一个服务器做,只需要一个支持增删和查询后继操作的有序表即可。可以证明大概率负载均衡。

Count-min Sketch

考虑这样一个问题:输入 \(N\)\([n]\) 上的整数,输入结束后给定 \(i\),要求近似回答 \(i\) 出现的次数。

Count-min Sketch 是一个数据结构,支持对于任意的误差参数 \(\epsilon\in (0,1)\),满足插入 \(N\)\([n]\) 上的元素后给定 \(x\in[n]\) 查询 \(x\) 共出现多少次,设 \(\tilde{C}\) 表示估计量,\(C\) 表示真实量,可以以大概率(\(1-1/\mathrm{poly}(n)\) 的概率)满足 \(C\le \tilde{C}\le C+\epsilon\cdot N\)

空间使用 \(\frac{\mathrm{poly}\log n}{\epsilon}\),且单次查询和插入单个元素的时间为 \(\frac{\mathrm{poly}\log n}{\epsilon}\)

大体思路就是设置一个随即哈希函数 \(h:[n]\to [m]\),然后对每一个 \(j\in [m]\) 统计 \(C[j]\) 为满足 \(h(x)=j\) 的数量。容易发现如果不存在哈希冲突最后 \(C[h(x)]\) 就是 \(x\) 出现次数,如果出现冲突那么 \(C[h(x)]\) 只会更大。

\(x\) 真实出现次数为 \(C_x\),那么有:

\[\begin{aligned} \mathrm{E}[C[h(x)]]&=\mathrm{E}\left[C_x+\sum_{y\ne x,h(y)=h(x)}C_y\right]\\ &=C_x+\sum_{y\ne x,h(y)=h(x)}C_y\\ &=C_x+\sum_{y\ne x}\mathrm{Pr}[h(x)=h(y)]C_y\\ &\le C_x+N/m \end{aligned} \]

Markov inequality:若 \(X\ge 0\),有 \(\mathrm{Pr}[X\ge a]\le \mathrm{E}[X]/a\)

\(C[h(x)]-C_x\) 用 Markov inequality,有 \(\mathrm{Pr}[C[h(x)]-C_x> \epsilon\cdot N]\le \frac{N/m}{\epsilon\cdot N}\le \frac{1}{\epsilon m}\),这里令 \(m\)\(\frac{2}{\epsilon}\) 即可得到 \(0.5\) 的概率满足条件。

由于 \(C[h(x)]\) 总是高估,可以多次实验取最小值,每次至少 \(0.5\) 概率满足条件,所以独立运行 \(O(\log n)\) 次即可以 \(1-1/\mathrm{poly}(n)\) 的概率满足条件。

所以整体思路就是设置 \(T=O(\log n)\) 个独立随机哈希 \(h^{(i)}:[n]\to [m]\),其中 \(m=2/\epsilon\),并初始化 \(T\)\(m\) 元 counter,然后插入删除元素就在每个 counter 里面都做就行了,查询就返回所有 counter 对应位置的最小值。

并且可以发现由于我们仅仅只是统计了元素数量,所以 count-min sketch 是可以接受相加相减的。

Count-min Sketch 应用

近似 k-HH

考虑 k-Heavy hitter(k-HH) 问题:给定长度为 \(N\) 的在 \([n]\) 上的数组,求出现次数前 \(k\) 多的元素。这个问题的应用例如用来查找那些异常流量或者访问较多的页面。

近似 k-HH 问题:你需要找到所有出现至少 \(N/k\) 次的元素,并且你还可以返回其他元素,返回的其他元素出现不少于 \(N/k-\epsilon N\) 次。

基于 count-min 的方法做数据流上的近似 k-HH,直接使用 count-min sketch 来存每个元素出现次数就行了,每次新加一个元素就查一下出现次数看是否满足条件即可,你找到的所有元素数量大概率是不超过 \(k\) 的,所以直接用表记录这些元素即可,动态增加 \(N\) 也没有任何问题。

近似 Inner Product

考虑近似 Inner Product 问题:给定两个在 \([n]\) 上的 frequency vector \(a\)\(b\),估计 \(\langle a,b\rangle\)。frequency vector 就是一个 \(n\) 为向量,统计每个 \(i\in [n]\) 出现多少次。

直接维护对应 \(a\)\(b\) 的 count-min sketch,然后对于每个 \(1\le i\le T\) 直接求对应 bucket 的 inner product,得到 \(\hat{C}^{(i)}:=\sum_{j=1}^mC_a^{(i)}[j]\cdot C_b^{(i)}[j]\),直接令 \(\hat{C}:=\min_{1\le i\le T}\hat{C}^{(i)}\) 作为最后估计。最后有保证大概率 \(\langle a,b\rangle\le \hat{C}\le \langle a,b\rangle+\epsilon\cdot \lVert a\rVert _1 \lVert b \rVert _1\)

实现哈希函数(Universal Hashing)

要构造 \(h:[n]\to [m]\) 使得 \(\forall x\ne y\in [n],\mathrm{Pr}[h(x)=h(y)]\le \frac{1}{m}\)

不妨 \(m\) 为素数,然后令 \(k=\lceil\log_m n\rceil\),取 \(k\)\([0,m)\) 的均匀独立随机数 \(a_i\),给定 \(x\in [n]\),把 \(x\) 表成 \(m-\) 进制数,定义 \(h(x):=\sum_{i}a_ix_i \bmod m\),容易证明满足条件。

当然这样不能保证 \(h(x)\)\(h(y)\) 完全独立,但是做 count-min 足够了。

Bloom Filter

维护一个集合 \(S\),全集为 \([n]\),支持插入操作和询问 \(i\in [n]\) 是否在 \(S\) 内。

Bloom filter 以 \(O(n)\) 空间(乘以极小的常数),和 \(O(1)\) 时间,大概率 \(1-\delta\) 成功回答询问。

简单来说就是维护一个长度 \(m\) 的 0-1 数组 \(A\),用 \(T\) 个随机哈希函数 \(h^{(i)}:[n]\to [m]\),插入 \(x\) 时把 \(A[h^{(i)}(x)],i=1,2,\dots,T\) 都置为 \(1\),然后认为 \(x\) 在集合内当且仅当对于任意的 \(i=1,2,\dots,T\) 都有 \(A[h^{(i)}(x)]=1\)

如果不插入 \(x\),对于每个 \(i\),考察 \(A[h^{(i)}(x)]=0\) 的概率。当插入某个 \(y\ne x\) 时,\(T\) 个哈希函数都有可能将其置为 \(1\),所以仍然为 \(0\) 的概率就是 \((1-1/m)^{T}\),最多插入 \(n\) 个,所以最终 \(A[h^{(i)}(x)]=0\) 的概率不超过 \((1-1/m)^{Tn}\approx \exp(-nt/m)\)。于是错误回答 Yes 的概率至多为 \((1-\exp(-nt/m))^T\)\(n,m\) 固定求导可以得到 \(T=\frac{m}{n}\ln 2\) 时概率最小,令失败概率为 \(\delta\) 得到 \(m=\frac{n\ln (1/\delta)}{(\ln 2)^2},T=\log_2(1/\delta)\)

如果要支持删除就不能仅仅维护 0-1 变量而要维护出现次数。

Karp-Rabin Hash

就是普通的字符串哈希。

分析略。

distance

集合相似度度量

两个集合 \(A,B\sube [n]\),至少有一个非空集,定义两个集合的 Jaccard Similarity 为 \(J(A,B)=\frac{|A\cap B|}{|A\cup B|}\)

如果数据有很多集合,要快速查询两个集合间的相似度,有什么好做法?

MinHash

考虑构造一个随机哈希函数 \(h:[n]\to [0,1]\),注意这里映射到实数,然后对于两个集合 \(A,B\),显然有 \(\mathrm{Pr}[\min_{x\in A}h(x),\min_{y\in B}h(y)]=J(A,B)\)

所以完整做法就是采用 \(T\) 个独立随机哈希 \(h^{(1)},\dots,h^{(T)}\),然后对于每个 \(A_i\)\(h^{(t)}\),计算 \(h^{(t)}_{\min}(A_i)=\min_{x\in A_i}h^{(t)}(x)\)。查询 \(J(A_i,A_j)=\frac{1}{T}\sum_{t=1}^T[h_{\min}^{(t)}(A_i)=h_{\min}^{(t)}(A_j)]\)

分析略,对于确定误差范围,有 \(T=O(\log n)\) 时失败概率 \(1/\mathrm{poly}(n)\)

向量相似度度量

考虑两个向量 \(x,y\in \mathbb{R}^d\),定义 cosine similarity 为 \(\sigma(x,y):=\frac{\langle x,y\rangle}{\lVert x\rVert _2\lVert y\rVert _2}\)

SimHash

类似 MinHash,尝试找到一个 \(h\)\(\sigma(x,y)\)\(\mathrm{Pr}[h(x)=h(y)]\) 相关。

生成一个 \(d\) 维随机高斯向量 \(w\)(每一维独立从标准正态分布中取然后归一化,等价于在 \(d\) 为单位球面取一个均匀随机点),定义 \(h(x):=\mathrm{sgn}(\langle w,x\rangle)\)\(\mathrm{sgn}\) 是符号函数,这里不妨令 \(\mathrm{sgn}(x)=[x\ge 0]\)),那么有 \(\forall x\ne y\in \mathbb{R}^d,\mathrm{Pr}[h(x)\ne h(y)]=\frac{\theta(x,y)}{\pi}\),其中 \(\theta(x,y)\) 为两个向量的夹角。

同样的,取 \(T\) 个独立的 SimHash 取平均值,也就是认为 \(\frac{\theta(x,y)}{\pi}=\frac{1}{T}\sum_{t=1}^T[h^{(t)}(x)\ne h^{(t)}(y)]\),令 \(f(x):=(h^{(1)}(x),\dots,h^{(T)}(x))\),可以看作是将 \(\mathbb{R}^d\) 上的点对应到 \(T\) 维的 Hamming 距离,也就是一个二进制串,两者距离定义为不同的位置数量,即可以认为有 \(\mathrm{dist}_H(f(x),f(y))\approx \theta(x,y)/\pi\)

Hamming空间近似最近邻查询

梦回 NOI2021 day2t1()

\(n\)\(d\) 维 Hamming 空间 \(\mathbb{H}^d\) 上的点,误差参数 \(\epsilon > 0\),预处理时间 \(O(dn+n^{1+1/(1+\epsilon)})\),给定任何 \(q\in \mathbb{H}^d\) 要求在 \(O(n^{1/(1+\epsilon)})\) 时间内找到 \((1+\epsilon)-\) 近似最近邻查询。

预处理过程,设 \(N=O(n^{1/(1+\epsilon)})\),找到 \(N\)\(d\) 个维度的随机排列 \(\sigma _1,\dots,\sigma _N\),然后对每个排列 \(\sigma _i\),将所有坐标按照 \(\sigma _i\) 排列后的 \(n\) 个数据点按字典序排序。

查询过程,假如查询串是 \(q\),对每个排列 \(\sigma _i\),将 \(q\) 坐标按照 \(\sigma _i\) 排列后找到它前后字典序相邻最近的串,然后看这两个距离更新答案。

分析略。

Locality Snesitive Hashing (LSH)

称一个哈希函数 \(h\) 是一种 LSH,若 \(\mathrm{Pr}[h(x)=h(y)]\) 可以反映 \(x\)\(y\) 的一种相似度。容易发现前面所说的 MinHash 和 SimHash 都是 LSH 的一种。

一般 LSH 都可以映射到 Hamming 空间做。具体做法就是考虑一个映射到 \(\{0,1\}\) 的随机哈希 \(b:X\to \{0,1\}\),然后将 \(b\)\(h\) 复合。我们知道如果 \(x=y\)\(\mathrm{Pr}[b(x)=b(y)]=1\) 否则 \(\mathrm{Pr}[b(x)=b(y)]=0.5\),所以有 \(\mathrm{Pr}[b(h(x))=b(h(y))]=\frac{1}{2}(1+\mathrm{Pr}[h(x)=h(y)])\)。于是用 \(T\)\(b\)\(h\) 的复合即可映射到 Hamming 空间了。

euclidean(欧式距离)

欧氏距离,也即 \(\ell_2-\) 范数,最常用的一个距离:

\[\forall x,y\in \mathbb{R}^d,\lVert x-y\rVert _2=\sqrt{\sum_{i=1}^d(x_i-y_i)^2} \]

格点离散化

最简单最直接的近似思想,把整个空间划分成 \(\overbrace{\ell\times\ell\times\cdots\times\ell}^{d}\) 的小正方形,然后把所有数据点 round 到它所在格子的正中心,容易发现每个点最多移动 \(\sqrt{d}\ell\),并且在一个 \(\overbrace{R\times R\times\cdots\times R}^{d}\) 的区域内最多有 \(\lceil R/\ell\rceil^d\) 个正方形。

另外欧式空间还有一个比较重要的性质就是一个半径 \(r\)\(d\) 维球中最多可以存在 \((O(r/t))^d\) 个两两距离 \(\ge t\) 的点,这个用体积的思想很好说明。

近似欧式点集直径:先找到一个直径的 \(2-\) 近似值 \(T\),然后做 \(\ell:=\epsilon\cdot T/\sqrt{d}\) 的网格,在新的点集的暴力求直径,总的复杂度为 \(O(nd)+\epsilon^{-O(d)}\),于是可以在线性时间内求出 \((1+\epsilon)-\) 近似直径。

类似算法可以用来求近似最小包围球。

四分树

四分树是二维情况下的特殊称法,这里以二维为例。

构建方法很简单,就从最大的 \(\Delta\times\Delta\) 正方形开始,每次正中间划开两刀分成四个部分,一直递归下去构造。构建的时空复杂度可以为 \(O(nd\log \Delta)\) 或者 \(O(n2^d\log \Delta)\)

四分树可以用来做 \(\epsilon-\) 近似最近邻查询,问题的准确描述是给定点集 \(P\sube [\Delta]^d\),对于任意的 \(q\in [\Delta]^d\) 你要回答出 \(\hat{q}\in P\) 使得 \(\lVert \hat{q}-q\rVert \le (1+\epsilon)\cdot\min_{q\prime\in P}\lVert q\prime-q\rVert\)

算法很简单,就是在四分树上每个节点维护一个代表点,每次到一个节点就看一下代表点到 \(q\) 的最短距离,尝试更新答案,接着往下搜索,如果再往下的节点无法让答案更优就剪掉,可以证明查询时间复杂度是 \(O(\epsilon^{-d}\log\Delta)\)

点对离散化:WSPD

输入点集 \(P\sube[\Delta]^d\),设 \(|P|=n\),考虑所有点对距离有 \(O(n^2)\) 对,如果两个点集 \(S\)\(T\),它们间的最短距离要比内部距离大太多,那么我们就可以把 \(S\times T\) 上的距离等同于同一个值。

对于一个点集 \(P\),一个 \(\epsilon-\) WSPD 是一系列点集对 \(\{(A_i,B_i)\}_i\),满足 \(\max(\mathrm{diam}(A_i),\mathrm{diam}(B_i))\le \epsilon\cdot\mathrm{dist}(A_i,B_i)\),并且 \(\forall p\ne q\in P,\exists i,\mathrm{s.t.} p\in A_i,q\in B_i\)

基于四分树构造 WSPD 十分简单,就是直接递归构造两个点集 \(S\)\(T\) 之间的 WSPD,其中 \(S\)\(T\) 对应的都是四分树上的节点,如果 \((S,T)\) 可以直接放入 WSPD 就直接放,否则就把直径较大的点集拆分成它的子节点,然后递归构造即可。可以证明采取压缩四分树得到的 WSPD 共有 \(O(n\epsilon^{-d}\log \Delta)\) 个。

使用 WSPD 可以精确求最近点对,构造 \(\epsilon<1\) 的 WSPD,然后可能成为最近点对的 \((p,q)\) 显然只有可能是那些单点对着单点的点对。

另外 WSPD 好像也可以用来做之前第五届图灵杯的 T3,到时候再补吧。

WSPD 也可以近似直径,构造 \(\epsilon/2-\) WSPD,然后对于每个点集对,任意取左边一个点和右边一个点看看距离更新答案,就可以得到直径的 \((1-\epsilon)-\) 近似。

WSPD 还可以用来求 \((1+\epsilon)-\) Spanner。\(t-\) Spanner 是数据集 \(P\sube [\Delta]^d\) 的一个欧式子图,边集是 \(P\times P\) 的子集并且边权是两端点的欧氏距离,并且满足 \(\forall x,y\in P\)\(\lVert x-y\rVert\le \mathrm{dist}_H(x,y)\le t\cdot \lVert x-y\rVert\)。做法也很简单,就是求 \(\epsilon/4-\) WSPD,然后每个点集对任意找代表点连边,可以证明满足条件。

在构造得到的 \((1+\epsilon)-\) spanner 上求 MST 即可得到 \(P\) 上的一个 \((1+\epsilon)-\) 近似的 MST。

Tree Embedding

Tree embedding 是将数据集 \(P\)\(\mathbb{R}^d\) 映射到一棵树 \(T(V,E)\) 上,其中 \(P\sube V\),目的是用树上的距离来近似欧式距离,首先要求 \(\forall x,y\in P,\mathrm{dist}_T(x,y)\ge \lVert x-y\rVert\),性能衡量为 Distortion(\(\max_{x\ne y\in P}\frac{\mathrm{dist}_T(x,y)}{\lVert x-y\rVert}\)),越小越好。

随机平移构造一个四分树,然后只有父亲孩子连边,其中第 \(i\) 层边权重为 \(\sqrt{d}\cdot 2^i\),结论是 \(\forall x,y\in P,\mathrm{dist}_T(x,y)\ge \lVert x-y\rVert,\mathbb{E}[\mathrm{dist}_T(x,y)]\le O(dh)\cdot \lVert x-y\rVert\),其中 \(h\) 是四分树的高度。

虽然这个近似比大,但是每对点期望距离都能保持,并且构造的是一棵树,而且也没有 \(2^d\)

dim reduction - JL

简要介绍

在做各类问题的时候,可能会遇到维数爆炸(Curse of Dimensionality)的问题,就是维数大一点可能就会导致算法运行效率等急速增加,所以我们希望将维数降下来。

Johnson-Lindenstrauss Transform 是一种针对距离的有效降维方法,可以维持 \(n\)\(\mathbb{R}^d\) 上的数据点集 \(P\) 中任意两点的距离。

具体来说,给定误差 \(\epsilon\),JL 给定了一个映射 \(f:\mathbb{R}^d\to \mathbb{R}^m\),这里 \(m=O(\epsilon^{-2}\log n)\),使得 \(\forall x,y\in P,\lVert f(x)-f(y) \rVert _2\in (1\plusmn \epsilon)\cdot \lVert x-y \rVert _2\),并且可以在 \(O(dm)\) 的时间计算 \(f(x)\)

可以发现非常神奇的一点就是降到的维度和 \(d\) 无关,而只和 \(n\)\(\epsilon\) 有关,这说明对于近似算法来说,“高维”就是 \(O(\log n)\) 维。

构造方法

把 JL 限制放大,令 \(\forall x,y,\lVert f(x)-f(y)\rVert ^2\in (1\plusmn \epsilon)\cdot\lVert x-y\rVert ^2\)

对于一对点 \(x,y\in \mathbb{R}^d\),考虑 \(v=x-y\),设计一个线性随机哈希(随机映射)\(g:\mathbb{R}^d\to \mathbb{R}\),使得 \(\mathbb{E}[g(v)^2]=\lVert v\rVert ^2\),然后多次实验,对于 \(O(n^2)\) 对距离应用该结论,保证所有距离都成立。

具体而言,构造 \(m\) 个独立的向量 \(w_1,\dots,w_m\),其中 \(w_i\) 的每一维都是从标准正态分布 \(\mathcal{N}(0,1)\) 中采样,然后令 \(f(v):=\frac{1}{\sqrt{m}}(\langle w_1,v\rangle,\dots,\langle w_m,v\rangle)\)(注意这里是直接拼接)。

固定 \(v\),首先考虑对于某个 \(w_i\),有 \(\mathbb{E}[\langle v,w_i\rangle]=\mathbb{E}[\sum_{j=1}^dv_jw_{i,j}]=0\),而 \(\mathbb{E}[\langle v,w_i\rangle ^2]=\mathbb{E}[\sum_{j=1}^dv_jw_{i,j}\sum_{k=1}^dv_kw_{i,k}]=\sum_{j=1}^d\sum_{k=1}^dv_jv_k\mathbb{E}[w_{i,j}w_{i,k}]\),而我们知道 \(\mathbb{E}[w_{i,j}w_{i,k}]=[j=k]\),所以有 \(\mathbb{E}[\langle v,w_i\rangle ^2]=\lVert v\rVert ^2\)

于是有 \(\mathbb{E}[f(v)^2]=\frac{1}{m}\sum_{i=1}^m\mathbb{E}[\langle v,w_i\rangle ^2]=\lVert v\rVert ^2\),起到了多次实验取平均的效果。

JL 的核心结论是,\(\forall v\in \mathbb{R}^d,\mathrm{Pr}[\lVert f(v)\rVert ^2\in (1\plusmn \epsilon)\cdot \lVert v\rVert ^2]\ge 1-\exp(-O(\epsilon ^2m))\)

然后对 \(n^2\) 对点对距离使用 Union Bound,一对失败概率 \(\exp(-O(\epsilon ^2m))\)\(n^2\) 失败概率不超过 \(\exp(-O(\epsilon ^2m))\cdot O(n^2)\),令该值小于 \(\delta\) 为常数,则有 \(m=O(\epsilon ^{-2}\log n)\)

当然实际中没人管 \(m\) 具体要多大,一般是试试哪个好用能用就行。

更简单的 JL 是将随机向量 \(w_i\) 替换成每维是独立 \(\{-1,1\}\) 均匀随机数的向量,因为这个随机也满足期望为 \(0\),方差为 \(1\)

应用:Linear regression

最小二乘法,给定矩阵 \(X\in \mathbb{R}^{n\times d}\) 和向量 \(y:=(y_1,\dots,y_n)^T\),则要求向量 \(w\in \mathbb{R}^d\) 使得 \(\lVert Xw-y\rVert ^2\) 尽量小,相当于是要解方程 \(X^TXw=X^Ty\)。当 \(X^TX\) 可逆时有唯一解 \(w=(X^TX)^{-1}X^Ty\)

方便起见不妨假设 \(X^TX\) 可逆,那么直接做的复杂度为 \(O(nd^2+d^3)\)。优化复杂度考虑降维,这里的 \(Xw\)\(y\) 都是 \(n\) 维向量,我们优化的目标是这两个向量的距离,想到构造一个 JL 矩阵 \(A\in \mathbb{R}^{m\times n}\),使得 \(\lVert AXw-Ay\rVert ^2\in (1\plusmn \epsilon)\lVert Xw-y\rVert ^2\)。这样可以让 \(O(nd^2)\) 复杂度降为 \(O(md^2)\)

但是考虑到我们相当于是要对所有向量使得它们有距离保证,不能再用 Union Bound 来分析 \(m\) 取多大,不过我们有子空间版本的 JL,设 \(\mathcal{U}\)\(\mathbb{R}^n\) 的一个 \(d\) 维子空间,则存在 \(m=O\left(\frac{d\log(\epsilon^{-1})+\log (\delta^{-1})}{\epsilon^2} \right)\) 的 JL 矩阵 \(A\in \mathbb{R}^{m\times n}\) 满足 \(\forall v\in \mathcal{U},\mathrm{Pr}[\lVert Av\rVert ^2\in (1\plusmn \epsilon)\cdot \lVert v\rVert ^2]\ge 1-\delta\)。如果 \(\epsilon,\delta\) 都认为是常数,那么就可以认为 \(m=O(d)\),这样带入求 \(X^TX\) 的复杂度就可以降为 \(O(d^3)\),但是如果得到的 JL 十分稠密还是要求 \(AX\)\(O(nd^2)\) 的,没有改进。

(其实这里不太严谨,因为我们要最小化的是 \(Xw\)\(y\) 的距离,而所有可能的 \(Xw-y\) 不能说是一个 \(d\) 维子空间,而应该是一个 \(d\) 维子空间整体往一个方向平移了部分,不过不影响结论成立)

使用一个随机 \(S\in \mathbb{R}^{m\times n}\) 的矩阵,这个矩阵每列只有一个均匀随机位置非 \(0\),且该元素均匀随机取 \(1\)\(-1\),然后用 \(\sqrt{\frac{n}{m}}S\) 替代原来的 JL 矩阵,可以证明 \(m=O(\epsilon ^{-2}d\mathrm{poly}\log (\epsilon ^{-1}d))\) 时有类似上面子空间版本的 JL 矩阵保证。

这个矩阵是十分稀疏的,以至于我们可以以 \(O(\mathrm{nnz}(X))\) 的优秀复杂度计算 \(AX\),其中 \(\mathrm{nnz}\) 表示 \(X\) 的非 \(0\) 元素数量(number of non-zero)。

dim Reduction - PCA

简要介绍

上面提到的 JL 降维方法容易发现是和数据无关的,这是好事但是也说明降维过程并没有用到数据本身的隐藏的依赖关系,这里要讲的 PCA 降维是通过数据间的隐含关系实现的。

PCA 的大体目标是给定 \(n\)\(d\) 维数据点 \(x_1,\dots,x_n\in \mathbb{R}^d\),你要找到 \(m\)\(d\) 维“基” \(v_1,\dots,v_m\in \mathbb{R}^d\),使得每个数据点可以近似用这些基表示,或者说分布在这个基张成的子空间附近。

更加具体地说,目的是找到 \(m\) 个单位正交(orthonormal)的向量 \(v_1,v_2,\dots,v_m\in \mathbb{R}^d\),最大化 \(\frac{1}{n}\sum_{i=1}^n\sum_{j=1}^m\langle x_i,v_j\rangle ^2\)。最优的这 \(m\) 个向量叫做 \(x_1,x_2,\dots,x_n\) 的前 \(m\) 个主成分(principal components,其实 PCA 全称就是 principal component analysis 主成分分析)。

PCA 求解

先考虑 \(m=1\) 的情况,即找到单位向量 \(v\),并最大化 \(\frac{1}{n}\sum_{i=1}^n\langle x_i,v\rangle ^2\),把 \(x_i\) 拍成一列得到一个矩阵 \(X\in \mathbb{R}^{n\times d}\),则就是要最大化 \(\lVert Xv\rVert ^2=(Xv)^TXv=v^TX^TXv\),令 \(A=X^TX\),则目标等价于最大化 \(v^TAv\)

容易发现 \(A\) 是实对称矩阵,并且 \(A\) 还是半正定矩阵,所以 \(A\) 可以正交对角化成 \(D=\mathrm{diag}(\lambda _1,\lambda _2,\dots,\lambda _d)\),且 \(\lambda _1\ge \lambda _2\ge \cdots\ge \lambda _d\ge 0\)。存在正交矩阵 \(Q\),有 \(A=QDQ^T\),于是最大化 \(v^TAv=v^TQDQ^Tv=(Q^Tv)^TD(Q^Tv)\),而 \(Q^Tv=\epsilon _1\) 显然最优,所以 \(m=1\) 是就是取 \(Q\) 的第一列即可。

事实上,前 \(m\) 个主成分就是 \(Q\) 的前 \(m\) 列。

\(A=X^TX\) 的所有特征向量,精确求解就是 \(O(nd^2)\) 的。

近似求解

近似求解可以做到线性时间 \(O(\mathrm{nnz}(X))\)。考虑到 \(A\) 本质上就是把单位球拉伸成椭球,所以任取一个向量 \(u\),多次应用 \(A\) 之后那么 \(u\) 的方向就会十分贴近最大的主成分的方向,所以算法就是随便选一个方向向量 \(u_0\),然后令 \(u_i=X^TXu_{i-1}\) 等到 \(u_i/\lVert u_i\rVert\approx u_{i-1}/\lVert u_{i-1}\rVert\) 即可。

结论是若 \(v_1\) 是最大的主成分,那么对于任意的 \(t\ge 1\),随机 \(u_0\) 每一维独立取 \(\mathcal{N}(0,1)\) 中的一个采样然后归一化,有常数概率使得 \(|\langle u_t,v_1\rangle|\ge 1-O(d)\left(\frac{\lambda _k}{\lambda _1}\right)^t\)(这里 \(k\) 是满足 \(\lambda _k<\lambda _1\) 的最小的 \(k\))。

如果要找第二大的主成分就先找最大的,然后数据点 \(x_i\gets x_i-\langle x_i,v_1\rangle v_1\),将该分量消除,然后再找最大的即可。前 \(m\) 大同理。

也可以把 \(u_0\) 取成每一维均匀 \(\plusmn 1\) 然后归一化,同样有类似保证。

linear programming(线性规划)

压缩感知简介

稀疏数据,非 \(0\) 元素较少的数据在现实中很常见,前面说的利用稀疏性都是拿到原始数据然后考虑降维,不过我们可以考虑不把所有数据读进来,而仅读需要的部分,也就是这里提到的压缩感知(compressive sensing)算法。

具体来说就是设计 \(m\) 个线性测量(linear measurements)\(a_1,\dots,a_m\in \mathbb{R}^n\),有一个未知信号 \(z\in \mathbb{R}^n\),给定线性测量上的值 \(b\),其中 \(b_i=\langle a_1,z\rangle\),要从 \(b\) 恢复 \(z\)。或者说构造矩阵 \(A\),从 \(b=Az\) 恢复出 \(z\)

事实上,求解 \(z\) 本质上相当于求 \(Ax=b\) 的解,如果 \(m=n\)\(A\) 可逆就可以精确解出 \(z\)。一般来说我们都是考虑 \(m<n\) 的情况,解不唯一,还需要利用 \(z\) 的稀疏性来找到 \(z\)

压缩感知定理是设 \(A\in \mathbb{R}^{m\times n}\),其中 \(m=O(k\log (n/k))\)\(A\) 的每个元素从 \(\mathcal{N}(0,1)\) 中独立采样,则 \(A\) 有高概率同时对所有 \(n\)\(k-\) 稀疏的 \(z\)\(b=Az\) 中精确恢复出 \(z\)\(k-\) 稀疏的定义是至多有 \(k\) 个非零位置,并且一般假定 \(k\ll n\)

\(z\) 中最大的 \(k\) 个元素绝对值和“远大于”其它元素绝对值和可以近似恢复出 \(z\)

具体如何恢复呢?当 \(m<n\)\(Ax=b\)\(x\) 是不唯一的,想到利用稀疏性我们便希望找到一个 \(x\) 满足非零坐标尽可能少,非零坐标数量也叫做 \(x\)\(\ell _0-\) 范数,所以这个优化问题也叫做 \(\ell _0-\) 最小化问题。

但是 \(\ell _0-\) 最小化问题是 NP-hard 的,考虑用 \(\ell _1-\) 范数替换,所以我们最终的问题就是给定 \(b\),找到 \(\lVert x\rVert _1\) 最小的 \(x\) 满足 \(Ax=b\)

线性规划

给定 \(n\) 个变量 \(x_1,\dots,x_n\),和若干个线性约束 \(\sum_{j=1}^na_{i,j}x_j\le b_i\),最小化目标函数 \(\sum_{j=1}^nc_jx_j\)

单纯性法不会,寄了()

streaming(数据流算法)

不会,反正不考()

MPC(大数据并行计算模型)

不会,反正不考()

final(期末总结安排)

直接开摆!!!

posted @ 2023-06-18 12:12  xiaolilsq  阅读(119)  评论(0编辑  收藏  举报