符号记法
- $[n]$：表示整数集合$\{1,\dots,n\}$
部分记法与原文有出入
概述
$(n,\lambda)$-MSM表示计算下列等式，其中$(P_1,\dots,P_n)$为$n$个群元素，$\{k_i\}_{i\in [n]}$为$\lambda$比特的标量
$$ Q = \sum_{i = 1}^n k_i\cdot P_i $$
实际MSM算法的$\lambda$范围可能在254~768之间，$n$则是百万级别的数，如果群元素$\{P_i\}$是椭圆曲线群上的点，则MSM的计算效率会更慢
MSM的一个主要的应用就是zkSNARK，用于计算证明与验证等式上的各个点，甚至部分zkSNARK方案中有超过七成的计算开销都来自于计算MSM，因此MSM算法的快慢很大程度上决定了这些zkSNARK方案的效率
Pippenger于1976年提出了一个MSM算法，至今仍然非常好用，具体实现可以看之前的那篇Blog
这里放一张比较形象的图，简单回顾一下这个方法，主要分为三个步骤
第一步是将原始任务拆分为多个小的子任务，这里需要设置一个窗口大小$s$，将每个$\lambda$比特的标量$k_i$拆分为$\frac{\lambda}{s}$个部分，每个部分为$s$比特的标量，记为$m_{ij}$，表示$k_i$拆分后的第$j$个子标量
对于划分后的子标量$m_{ij}$，可以得到对于任意的标量$k_i$，有
$$ k_i=\sum_{j=1}^{\frac{\lambda}{s}}2^{(j-1)s}\cdot m_{ij} $$
然后基于这些子标量$m_{ij}$，可以将第$i$个群元素$P_i$的计算转变为下列等式，也就是上图中左侧的中间部分（$j\in [\frac{\lambda}{s} ]$）
$$ G_j = \sum^n_{i = 1}m_{ij}\cdot P_i $$
然后原始MSM等式$Q = \sum_{i = 1}^n k_i\cdot P_i$就可以转变为下面这样
$$ \begin{aligned} \begin{aligned}{Q}=\sum_{i=1}^n{Q}_i=\sum_{i=1}^nk_i\cdot {P}_i\end{aligned}& =\sum_{i=1}^n\sum_{j=1}^{\frac{\lambda}{s}}(2^{(j-1)s}\cdot m_{ij})\cdot {P}_i \\ &=\sum_{j=1}^{\frac\lambda s}2^{(j-1)s}\cdot \sum_{i=1}^nm_{ij}\cdot {P}_i \\ &=\sum_{j=1}^{\frac\lambda s}2^{(j-1)s}\cdot {G}_j \end{aligned} \tag 1 $$
基于上面这个等式，算法的第二步是分别计算拆分后每个子任务的标量乘法，也即计算每个$G_j(j\in [\frac{\lambda}{s} ] )$
这里有个小trick：构造若干个桶（桶的数量取决于窗口大小$s$），然后把相关联的$P_i$放到对应的桶中
上面这个图的例子的窗口大小$s = 4$，因此桶的数量为$2^s -1= 2^4 -1= 15$个，并分别编号为$B_1,\dots, B_{15}$
然后观察$P_1,\dots,P_5$所对应的子标量（也就是上图左侧中间的三个框），首先是第一个子标量（红色框），子标量$(1101)_2$分别对应于群元素$P_1$和$P_5$，因此把$P_1,P_5$放到$B_{13}$这个桶里，剩余的两个标量同理，将$P_2$放入$B_{14}$，$P_3,P_4$放入$B_{15}$，这个过程对应于上图右上角的红色部分，另外两部分的子标量同理，分别对应于绿色和蓝色部分
这里放入桶中的意思就是点加法，$B_{13}$存放的实际上是$P_1+P_5$的结果
下面蓝色部分同理，$B_{10}$存放的是$P_1+P_2+P_4+P_5$
接下来就是将每个部分的桶求和，第$j$部分的桶求和方式为
$$ {G}_j=\sum_{l=1}^{2^s-1}l\cdot {B}_l^{(j)} $$
以上图的红色部分为例，实际上就是计算$G_3 = 13\cdot B_{13} + 14\cdot B_{14}+ 15\cdot B_{15} = 13\cdot (P_1+P_5) +14\cdot P_2 + 15\cdot (P_3+P_3)$
最后一步就是将每个桶的求和结果乘上一个幂系数并求和，幂系数取决于窗口大小$s$，求和方式如下
$$ {Q}=\sum_{j=1}^{\frac{\lambda}{s}}2^{(j-1)s}\cdot {G}_j=2^s\left(\cdots\left(2^s\left(2^s\cdot {G}_{\frac{\lambda}{s}}+{G}_{\frac{\lambda}{s}-1}\right)+{G}_{\frac{\lambda}{s}-2}\right)\cdots\right)+{G}_1 \tag 2 $$
上图中$s=4$，划分为了三个部分，因此幂系数分别为$2^0,2^4,2^8$，因此求和结果为左下角的$Q = 2^8\cdot G_3+2^4\cdot G_2+G_1$
Elastic MSM Algorithm
这里观察一下上面的式1的第一行
$$ {Q}=\sum_{i=1}^n{Q}_i=\sum_{i=1}^nk_i\cdot {P}_i=\sum_{i=1}^n\sum_{j=1}^{\frac{\lambda}{s}}(2^{(j-1)s}\cdot m_{ij})\cdot {P}_i $$
这里将$\frac{\lambda}{s}$拆分成两个整数$w,k$，满足$\frac{\lambda}{s} = w\cdot k$，此时$Q_i$的计算过程可以转变为一个由2的幂构成的行向量与子标量构成的列向量相乘的结果，如下
$$ {Q}_i=\begin{pmatrix}1&2^s&2^{2s}&\cdots&2^{(wk-1)s}\end{pmatrix}\cdot\begin{pmatrix}m_{i1}{P}_i\\m_{i2}{P}_i\\m_{i3}{P}_i\\\vdots\\m_{i(wk)}{P}_i\end{pmatrix} $$
接下来定义$M_{i(l,t)} = m_{i((l-1)k+t)}$，其中$l\in [w],t\in [k]$，则$Q_i$的计算方式可以转变为下面这样
$$ {Q}_i=\sum_{j=1}^{wk}\left(2^{(j-1)s}m_{ij}\right)\cdot {P}_i=\sum_{l=1}^w\sum_{t=1}^k2^{((l-1)k+(t-1))s}\cdot M_{i(l,t)}\cdot {P}_i $$
基于上面这个等式，$Q_i$的计算方式可以变成下面这个矩阵乘法
$$ \begin{pmatrix}1&2^{ks}&\cdots&2^{(w-1)ks}\end{pmatrix}\cdot{P}_i\cdot\begin{pmatrix}M_{i(1,1)}&M_{i(1,2)}&\cdots&M_{i(1,k)}\\M_{i(2,1)}&M_{i(2,2)}&\cdots&M_{i(2,k)}\\\vdots&\vdots&&\vdots\\M_{i(w,1)}&M_{i(w,2)}&\cdots&M_{i(w,k)}\end{pmatrix}\cdot\begin{pmatrix}1\\2^s\\\vdots\\2^{(k-1)s}\end{pmatrix} \tag 3 $$
这里定义$M_{i(l,t)}$的意思是将$m_{ij}$重新排列成一个$w$行$k $列的矩阵，没有这么复杂
接下来再定义下面三个辅助变量
$$ \begin{aligned} {P}_{ij}&=2^{(j-1)ks}\cdot {P}_i\\ {G}_{il}&=\sum_{j=1}^wM_{i(j,l)}\cdot {P}_{ij}\\ N_{ij}&=\sum_{l=1}^k2^{(l-1)s}\cdot M_{i(j,l)} \end{aligned} $$
这里解释一下三个辅助变量的含义，这里有$i\in [n],j\in [w]$
- ${P}_{ij}$：$P_i$乘上一个2的幂，其实就是式3中$P_i$与行向量$\begin{pmatrix}1&2^{ks}&\cdots&2^{(w-1)ks}\end{pmatrix}$第$j$项的乘积
- ${G}_{il}$：式3的$M$矩阵中第$l$列中所有元素与$P_{ij}$相乘再求和
- $N_{ij}$：式3的$M$矩阵中第$j$列与列向量$\begin{pmatrix}1&2^s&\cdots &2^{(k-1)s}\end{pmatrix}^\intercal$的乘积
基于这三个辅助变量，式3中计算$Q_i$的方式可以改写为如下形式
$$ \begin{aligned} {Q}_i=\sum_{j=1}^wN_{ij} \cdot {P}_{ij}& =\sum_{j=1}^w\sum_{l=1}^k2^{(l-1)s}\cdot M_{i(j,l)}\cdot {P}_{ij} \\ &=\sum_{l=1}^k2^{(l-1)s}\cdot \sum_{j=1}^wM_{i(j,l)}\cdot {P}_{ij} \\ &=\sum_{l=1}^k2^{(l-1)s}\cdot {G}_{il} \end{aligned} $$
然后再把所有的$Q_i$求和就可以得到原始的结果，也即$Q = \sum^n_{i= 1} Q_i$
下面这张图解释了如何从原版的MSM转换为Elastic版本的MSM，所有的$P_{ij}$都可以在于处理阶段进行计算
The Elastic Pippenger Algorithm
前面提到了MSM是zkSNARK中最耗时的操作，为了提高算法的效率，并行执行算法是一个非常好的选择
正好Pippenger算法的并行执行比其他MSM算法性能更好，这一节介绍一下如何将Elastic MSM应用到两类并行的Pippenger算法和一种更先进的MSM算法中，相当于将Elastic MSM作为一个模块来替换其他算法中MSM的计算部分
Naive Approach 1
第一种并行方法最朴素，观察式1的计算过程
$$ Q =\sum_{j=1}^{\frac\lambda s}2^{(j-1)s}\cdot {G}_j $$
由于计算过程被分成了$\frac{\lambda}{s}$个部分，因此最简单的办法就是并行$\frac{\lambda}{s}$个任务计算对应的$G_j$，最后再将每个任务的计算结果乘上对应的系数$2^{(j-1)s}$再相加即可
这个方法简单是简单，但是加速效果至多只能达到$\frac{\lambda}{s}$倍，实际应用中$\lambda$的取值取决于所选的椭圆曲线，即便是合理选择分块大小$s $，加速效果也有限，而且这个方法也不能完全使用显卡的算力，因此不是一种好办法
Naive Approach 2
这种方法也比较朴素，假设$n$可以被$t$整除，直接将$Q = \sum^n_{i = 1}k_i\cdot P_j$拆分成$t $个子任务$Q_j,j\in [t]$，最后再合并每个子任务的计算结果$Q = \sum^t_{k = 1}Q_j$
这个方法不错，因为只要拆分的任务大小比较合理，就可以在显卡中并行足够多的MSM子任务，直到占满显卡的算力为止
Naive Approach 3
第三种方法相当于缝合了上面两种方法，既拆成$\frac{\lambda}{s}$个对应的$G_j$，又拆成$t$个子任务
$$ {Q}=\sum_{i=1}^nk_i{P}_i=\sum_{j=1}^{t/\frac\lambda s}{Q}_j=\sum_{j=1}^{t/\frac\lambda s}\sum_{l=1}^{\frac\lambda s}2^{(l-1)s}{G}_{j,l} $$
Benchmark
后续的效率分析可以看原文$[1]$