|
15 | 15 |
|
16 | 16 | # U统计量计算新算法
|
17 | 17 |
|
18 |
| -## U统计量的定义 |
| 18 | +## U统计量和V统计量 |
19 | 19 | 在统计理论中,$U$统计量一直是研究者的“心头好”。它最大的优势在于:能够轻松构造**无偏估计**,名字里的 $U$ 其实就是 “unbiased” 的缩写。
|
20 | 20 |
|
21 | 21 | 不过,虽然原理简单,$U$统计量的计算却常常令人头疼。先来看$m$阶$U$统计量的一般形式:
|
|
24 | 24 | \mathbb{U}_{m}(h) = \frac{1}{n(n-1)(n-2)\cdots(n-m-1)}\sum_{1 \le i_1 \neq i_2 \neq \cdots \neq i_m \le n} h(X_{i_1}, X_{i_2}, \cdots, X_{i_m}),\qquad\qquad (1)
|
25 | 25 | $$
|
26 | 26 |
|
27 |
| -其中,$X_1, X_2, \cdots, X_n \in \Omega \subseteq \mathbb{R}^{p}$ 是 $n$ 个独立同分布的样本;$h : \Omega^m \to \mathbb{R}$ 称为核函数;指标条件 $i_1 \neq i_2 \neq \cdots \neq i_m$ 表示取的$m$个样本互不相同因此互相独立。 |
| 27 | +其中,$X_1, X_2, \cdots, X_n \in \Omega \subseteq \mathbb{R}^{p}$ 是 $n$ 个独立同分布的样本;$h : \Omega^m \to \mathbb{R}$ 称为核函数;指标条件 $i_1 \neq i_2 \neq \cdots \neq i_m$ 表示取的$m$个样本**互不相同**因此互相独立。 |
28 | 28 |
|
29 | 29 | $U$统计量之所以重要,在于它天然满足无偏性:
|
30 | 30 |
|
|
39 | 39 | \mathbb{V}_{m}(h) = \frac{1}{n^m}\sum_{1 \le i_1, i_2,\cdots, i_m \le n} h(X_{i_1}, X_{i_2}, \cdots, X_{i_m}).\qquad\qquad (2)
|
40 | 40 | $$
|
41 | 41 |
|
42 |
| -注意$U$和$V$的区别就在于计算中指标是否存在限制,$V$统计量的计算中对指标$(i_1, i_2,\cdots, i_m)$没有任何限制,因此总共有$n^m$个值。 |
| 42 | +注意$U$和$V$的区别就在于计算中指标是否存在限制,$V$统计量的计算中对指标$(i_1, i_2,\cdots, i_m)$**没有任何限制**,因此总共有$n^m$个值。 |
43 | 43 |
|
44 |
| -很容易可以看出来,$V$统计量可以表示成一组$U$统计量的线性组合(就是加加减减乘乘),反过来$U$统计量也可以表示成一组$V$统计量的线性组合(例子在稍后展示)。所以只要有一方可以高效计算,另一方就可以通过线性组合得到,我们算法就是这样子做! |
| 44 | +很容易可以看出来,$V$统计量可以表示成一组$U$统计量的线性组合(就是加加减减乘乘),反过来$U$统计量也可以表示成一组$V$统计量的线性组合(例子在稍后展示)。所以只要有一方可以高效计算,另一方就可以通过线性组合得到,我们的算法正是这个朴素的想法。 |
45 | 45 |
|
46 | 46 | ---
|
47 | 47 |
|
|
71 | 71 | ---
|
72 | 72 |
|
73 | 73 | ## 计算难题
|
74 |
| -回到公式 (1)。如果直接用 for 循环枚举$m$阶组合来计算,那么复杂度是 $O(n^m)$(假设$h$的计算复杂度与$n$无关)。 作为对比,矩阵乘法或矩阵求逆的复杂度是 $O(n^3)$,因此,一旦阶数 $m > 3$,$U$统计量的计算往往会变得难以承受。 |
| 74 | +回到公式 (1)和(2)。如果直接用 for 循环枚举$m$阶组合来计算,那么复杂度是 $O(n^m)$(假设$h$的计算复杂度与$n$无关)。 作为对比,矩阵乘法或矩阵求逆的复杂度是 $O(n^3)$,因此,一旦阶数 $m > 3$,$U$统计量和$V$统计量的计算往往会变得难以承受。 |
75 | 75 |
|
76 | 76 | 而偏偏,我们关注的统计量就是高阶$U$统计量。James M. Robins 等人在 2008 年提出的 **高阶影响函数(Higher Order Influence Function, HOIF)**,在各种假设下均能达到最优估计速率。HOIF 是一个高阶 $U$统计量,最优阶数大约是 $m \sim \sqrt{\log(n)}$。这意味着,HOIF 的实际应用必须直面高阶 $U$ 统计量的计算挑战。
|
77 | 77 |
|
|
89 | 89 | ---
|
90 | 90 | ## 新算法
|
91 | 91 |
|
92 |
| -接下来我们介绍新算法的实现。核心工具是 Python 库中提供的强大函数 [numpy.einsum](https://numpy.org/doc/stable/reference/generated/numpy.einsum.html) 和 [torch.einsum](https://pytorch.org/docs/stable/generated/torch.einsum.html)。 |
93 |
| -这两个函数的底层实现都经过高度优化,尤其是 `torch.einsum`,可以方便地利用 **CPU 与 GPU 并行**,在工程上极具优势。 |
| 92 | +接下来我们介绍新算法的实现。核心的想法已经在前文提到,把$U$统计量拆成$V$统计量的组合,再高效计算$V$统计量。 |
| 93 | + |
| 94 | +计算$V$统计量的核心工具是 Python 库中提供的强大函数 [numpy.Einsum](https://numpy.org/doc/stable/reference/generated/numpy.Einsum.html) 和 [torch.Einsum](https://pytorch.org/docs/stable/generated/torch.Einsum.html)。 |
| 95 | +这两个函数的底层实现都经过高度优化,尤其是 `torch.Einsum`,可以方便地利用 **CPU 与 GPU 并行**,在工程上极具优势。 |
94 | 96 |
|
95 | 97 | ---
|
96 | 98 |
|
97 | 99 | ### Einsum
|
98 | 100 |
|
99 | 101 | 那么,Einsum 究竟做了什么?
|
100 |
| -简单来说,Einsum 是一种针对张量的高效操作方式:输入若干张量,指定索引规则,对部分或全部指标执行无约束求和,最终输出一个新的张量或常数。 |
| 102 | +简单来说,Einsum 是一种针对张量的高效操作方式:输入若干张量,指定索引规则,对部分或全部指标执行**无约束**求和,最终输出一个新的张量或常数。 |
101 | 103 |
|
102 | 104 | 下表给出几个常见例子:
|
103 | 105 |
|
104 | 106 | | Einsum 调用格式 | 对应结果 |
|
105 | 107 | |:---------------------------------|:---------------------------------------------|
|
106 |
| -| `np.einsum('ij,jk->ik', A, B)` | $D_{ik} = \sum_j A_{ij} B_{jk}$ | |
107 |
| -| `np.einsum('ijk->i', X)` | $D_i = \sum_{j,k} X_{ijk}$ | |
108 |
| -| `np.einsum('ij,jk,kl-> ', A,B,C)`| $D = \sum_{i,j,k,l} A_{ij} B_{jk} C_{kl}$ | |
109 |
| -| `np.einsum('ijk,jjk,kkl-> ', A,B,C)`| $D = \sum_{i,j,k,l} A_{ijk} B_{jjk} C_{kkl}$| |
| 108 | +| `Einsum('ij,jk->ik', A, B)` | $D_{ik} = \sum_j A_{ij} B_{jk}$ | |
| 109 | +| `Einsum('ijk->i', X)` | $D_i = \sum_{j,k} X_{ijk}$ | |
| 110 | +| `Einsum('ij,jk,kl-> ', A,B,C)`| $D = \sum_{i,j,k,l} A_{ij} B_{jk} C_{kl}$ | |
| 111 | +| `Einsum('ijk,jjk,kkl-> ', A,B,C)`| $D = \sum_{i,j,k,l} A_{ijk} B_{jjk} C_{kkl}$| |
110 | 112 |
|
111 | 113 | 在表达式 `'ij,j->j'` 中,`->` 左边指定输入张量的指标分布,右边表示最终保留的指标。如果某个指标只出现在左边而没出现在右边,它就会被求和“消去”。
|
112 | 114 | 这个记号起源于 **Einstein 求和约定**(1916 年提出),矩阵与高阶张量的绝大多数运算都能通过 Einsum 表达,而在处理三阶及以上张量时,它往往是最简洁、最高效的工具。
|
|
116 | 118 | 读到这里,你可能已经注意到 Einsum 和 $V$ 统计量的天然联系。事实上,所有的 $V$ 统计量都可以直接用 Einsum 表示:
|
117 | 119 | 只要核函数 $h$ 能够写成乘法分解形式,$V$ 统计量的计算就可以转化为一个张量求和问题;如果 $h$ 不能分解,那就相当于一个整体张量,复杂度无法进一步降低。
|
118 | 120 |
|
119 |
| -一个典型的例子是高阶影响函数(HOIF)中的核函数: |
| 121 | +回到高阶影响函数(HOIF)中关心的核函数: |
120 | 122 |
|
121 | 123 | $$
|
122 | 124 | \begin{aligned}
|
|
164 | 166 | $$
|
165 | 167 | 它退化成了一个一阶 $V$ 统计量。
|
166 | 168 |
|
167 |
| -接下来看三阶的情况。为了简化符号,我们暂时省略掉指标范围 $1\sim n$ 以及归一化因子(如 $n(n-1)(n-2),n^3$ 等): |
| 169 | +接下来看三阶的情况。为了简化符号,我们暂时省略掉指标范围 $1\sim n$ 以及$\mathbb{U}$和$\mathbb{V}$的归一化因子(如 $\frac{1}{n(n-1)(n-2)},\frac{1}{n^3}$ 等): |
168 | 170 |
|
169 | 171 | $$
|
170 | 172 | \begin{aligned}
|
171 | 173 | \mathbb{U}_{3}[h] & = \sum_{i_1 \neq i_2 \neq i_3} h(X_{i_1},X_{i_2},X_{i_3}) \\[6pt]
|
172 | 174 | & = \Big(\sum_{i_1,i_2,i_3} - \sum_{(i_1=i_2)\ne i_3 } - \sum_{(i_1=i_3)\ne i_2 } - \sum_{(i_2=i_3)\ne i_1 } - \sum_{i_1=i_2=i_3}\Big) h(X_{i_1},X_{i_2},X_{i_3}) \\[6pt]
|
173 |
| - & = \mathbb{V}_{3}[h] - \mathbb{U}_{2,1=2}[h] - \mathbb{U}_{2,1=3}[h] - \mathbb{U}_{2,2=3}[h] - \mathbb{U}_{1,1=2=3}[h]. \qquad (5) |
| 175 | + & = \mathbb{V}_{3}[h] - \mathbb{U}_{2,1=2}[h] - \mathbb{U}_{2,1=3}[h] - \mathbb{U}_{2,2=3}[h] - \mathbb{U}_{1,1=2=3}[h]. \qquad (4) |
174 | 176 | \end{aligned}
|
175 | 177 | $$
|
176 | 178 |
|
177 | 179 | 在这里,$(i_1=i_2)\ne i_3$ 表示 $i_1=i_2=i$ 且 $i\neq i_3$。也就是说只剩下两个指标 $(i,i_3)$,对应一个二阶 $U$ 统计量,记作 $\mathbb{U}_{2,1=2}$(忽略归一化因子 $\tfrac{1}{n(n-1)}$)。其他几个符号同理。注意一阶 $U$ 统计量和一阶 $V$ 统计量是一样的。
|
178 | 180 |
|
179 |
| -因此,$\mathbb{U}_3$ 被拆成了 $\mathbb{V}_3$ 和更低阶的 $U$ 统计量。而低阶 $U$ 我们已经能写成 $V$ 的组合,所以顺着这个递推关系,就能得到一个一般的“$U$ 拆 $V$”算法。继续推下去: |
| 181 | +因此,$\mathbb{U}_3$ 被拆成了 $\mathbb{V}_3$ 和更低阶的 $U$ 统计量。而$2$阶 $U$和$1$阶 $U$ 我们已经能写成 $V$ 的组合,所以顺着这个递推关系,就能得到一个一般的“$U$ 拆 $V$”算法。继续推下去: |
180 | 182 |
|
181 | 183 | $$
|
182 | 184 | \begin{aligned}
|
|
192 | 194 |
|
193 | 195 | #### 更优雅的数学描述
|
194 | 196 |
|
195 |
| -但是,我们再仔细看看,能不能有个更优雅的数学刻画呢? 仔细观察公式 $(5)$,它其实表明:**一个 $V$ 统计量可以拆成若干 $U$ 统计量的和**,并且系数都是 $1$: |
| 197 | +但是,我们再仔细看看,能不能有个更优雅的数学刻画呢? 仔细观察公式 $(4)$,它其实表明:**一个 $V$ 统计量可以拆成若干 $U$ 统计量的和**,并且系数都是 $1$: |
196 | 198 |
|
197 | 199 | $$
|
198 | 200 | \mathbb{V}_{3}[h]= \mathbb{U}_{3}[h] + \mathbb{U}_{2,1=2}[h] + \mathbb{U}_{2,1=3}[h] + \mathbb{U}_{2,2=3}[h] + \mathbb{U}_{1,1=2=3}[h].
|
|
224 | 226 |
|
225 | 227 | ---
|
226 | 228 |
|
227 |
| -#### 一般公式 |
| 229 | +#### 更优雅的条件描述 |
228 | 230 |
|
229 | 231 | 上面的定义可以完美推广出去,有了这个数学的形式化,就能写出一个漂亮的恒等式。对任意 $m$ 阶:
|
230 | 232 |
|
|
241 | 243 | \pi_2 = \{ \{ 1\},\{ 2 \} \} > \{ \{ 1, 2\} \},\\
|
242 | 244 | \pi_3 = \{ \{ 1\},\{ 2 \}, \{3\} \} > \{ \{ 1, 2\}, \{3\} \} > \{ \{1,2,3\} \}.
|
243 | 245 | $$
|
244 |
| -所以精细的意思是,本来某个组里还可以再分,他却没有分,如果你把这一组又继续细分了,你就更精细,你就更大。但是并不是所有的分组方式都可以比较,比如: |
| 246 | +所以精细的意思是,本来某个组里还可以再分,他却没有分,如果你把这一组又继续细分了,你就更精细,你就更大。在$\{ \{1, 2\}, \{3\} \}$里继续切分$\{1,2\}$得到$\{ \{ 1\},\{ 2 \}, \{3\} \}$,你就得到了更精细的一个划分。但是并不是所有的分组方式都可以比较,比如: |
245 | 247 | $$
|
246 | 248 | \{ \{ 1, 2\}, \{3\} \},\{ \{ 1, 3\}, \{2\} \},\{ \{ 2, 3\}, \{1\} \} 均无法互相比较
|
247 | 249 | $$
|
|
254 | 256 | 并且很容易可以证明不只是$\pi_m$, 对任意的划分$\pi$, 都有这个关系:
|
255 | 257 |
|
256 | 258 | $$
|
257 |
| -\mathbb{V}[\pi](h) = \sum_{\rho \le \pi} \mathbb{U}[\rho](h), \forall \pi \in \Pi_m. \qquad (6) |
| 259 | +\mathbb{V}[\pi](h) = \sum_{\rho \le \pi} \mathbb{U}[\rho](h), \forall \pi \in \Pi_m. \qquad (5) |
258 | 260 | $$
|
259 | 261 |
|
260 | 262 | ---
|
|
282 | 284 | f(\pi) = \sum_{\rho \le \pi} g(\rho)\,\mu(\rho,\pi), \quad \forall \pi \in \Pi_m.
|
283 | 285 | $$
|
284 | 286 |
|
285 |
| -套用到我们的情形(令 $g(\pi) = \mathbb{V}[\pi](h), f(\pi) = \mathbb{U}[\pi](h)$),得到: |
| 287 | +套用到我们的情形(令 $g(\pi) = \mathbb{V}[\pi](h), f(\pi) = \mathbb{U}[\pi](h)$),我们已经有了公式$(5)$,由此得到: |
286 | 288 |
|
287 | 289 | $$
|
288 |
| -\mathbb{U}[\pi](h) = \sum_{\rho \le \pi} \mathbb{V}[\rho](h)\,\mu(\rho,\pi), \quad \forall \pi \in \Pi_m. \qquad (7) |
| 290 | +\mathbb{U}[\pi](h) = \sum_{\rho \le \pi} \mathbb{V}[\rho](h)\,\mu(\rho,\pi), \quad \forall \pi \in \Pi_m. \qquad (6) |
289 | 291 | $$
|
290 | 292 |
|
291 | 293 | 我们最终关心的是 $\pi=\pi_m$ 的情况:
|
|
300 | 302 | \mu(\pi,\pi_m) = \mu_{\pi} = (-1)^{m-K}\,(n_1-1)!\,(n_2-1)!\cdots(n_K-1)!.
|
301 | 303 | $$
|
302 | 304 |
|
303 |
| -比如回到$m=3$的例子,注意$n! = n(n-1)(n-2)\cdots 1,0! =1,(-1)^{0}=1$. |
| 305 | +比如回到$m=3$的例子,注意$n! = n(n-1)(n-2)\cdots 1,0! =1,(-1)^{0}=1$. |
304 | 306 | $$
|
305 | 307 | \begin{aligned}
|
306 |
| - \pi = \{\{1\},\{2\},\{3\}\}, \quad & \mu_{\pi} = (-1)^{3-3} (1 -1)! (1 -1)! (1 -1)! = + 1\\ |
307 |
| - \pi = \{\{1,2\},\{3\}\}, \quad &\mu_{\pi} = (-1)^{3-2} (2 -1)! (1 -1)! = - 1\\ |
308 |
| - \pi = \{\{1,3\},\{2\}\},\quad &\mu_{\pi} = (-1)^{3-2} (2 -1)! (1 -1)! = - 1\\ |
309 |
| - \pi = \{\{2,3\},\{1\}\}, \quad &\mu_{\pi} = (-1)^{3-2} (2 -1)! (1 -1)! = - 1\\ |
310 |
| - \pi = \{\{1,2,3\}\}, \quad & \mu_{\pi} = (-1)^{3-1} (3 -1)! = + 2\\ |
| 308 | + \pi = \{\{1\},\{2\},\{3\}\}, \quad & \mu_{\pi} = (-1)^{3-3} (1 -1)! (1 -1)! (1 -1)! = + 1,\\ |
| 309 | + \pi = \{\{1,2\},\{3\}\}, \quad &\mu_{\pi} = (-1)^{3-2} (2 -1)! (1 -1)! = - 1,\\ |
| 310 | + \pi = \{\{1,3\},\{2\}\},\quad &\mu_{\pi} = (-1)^{3-2} (2 -1)! (1 -1)! = - 1,\\ |
| 311 | + \pi = \{\{2,3\},\{1\}\}, \quad &\mu_{\pi} = (-1)^{3-2} (2 -1)! (1 -1)! = - 1,\\ |
| 312 | + \pi = \{\{1,2,3\}\}, \quad & \mu_{\pi} = (-1)^{3-1} (3 -1)! = + 2. |
311 | 313 | \end{aligned}
|
312 | 314 | $$
|
| 315 | + |
| 316 | +到此,我们的新算法在工程上已经清晰了,对于$m$阶$U$统计量,计算$\Pi_{m}$上的所有划分$\pi$和系数$\mu_{\pi}$,然后把$\mathbb{V}[\pi](h)$表示成Einsum带入 Python 的Einsum 函数,然后再汇合到一起就算出了$\mathbb{U}[h]$的值。 |
| 317 | + |
| 318 | +### 复杂度分析 |
| 319 | + |
| 320 | +接下里我们分析上述新算法的复杂度,这会将我们导入比莫比乌斯反演更美丽宏大的数学世界,简单图中的顶点消除(vertex elimination),树宽(treewidth), 商图(quotien graph)将完美刻画我们的计算过程与最终的计算复杂度。 |
0 commit comments