引言

近年来，随着深度学习技术的快速发展，大模型（Large Language Models, LLMs）在自然语言处理（NLP）领域取得了显著的进展。这些模型通过大规模数据训练，能够生成高质量的文本、完成复杂的语言理解任务，并在多个领域展现出强大的能力。然而，随着模型规模的不断扩大，如何高效地处理长序列数据、降低计算复杂度以及减少内存占用成为了研究的热点问题。

在大模型中，Attention 机制（注意力机制）是核心组件之一。自 Transformer 模型提出以来，Attention 机制凭借其强大的序列建模能力，逐渐取代了传统的循环神经网络（RNN）和卷积神经网络（CNN）。然而，标准的 Attention 机制（如 Multi-Head Attention, MHA）在处理长序列时面临着计算复杂度和内存占用的双重挑战，这使得其在处理长文本时效率低下。

为了解决这些问题，研究者们提出了多种改进的 Attention 机制，例如 FlashAttention、Multi-Query Attention (MQA)、Group Query Attention (GQA) 以及 Sparse Attention 等。这些方法通过优化计算流程、减少内存访问次数或引入稀疏性，显著提升了 Attention 机制的效率。此外，DeepSeek-V2 提出的 Multi-Head Latent Attention (MLA) 进一步通过低维映射和算子融合技术，降低了计算复杂度和内存占用，为大模型的规模化应用提供了新的思路。

本文将深入探讨这些 Attention 技术的原理、实现细节及其在大模型中的应用。我们将从标准的 MHA 出发，逐步介绍其变体（如 MQA 和 GQA），并详细分析 FlashAttention 如何通过分块计算和重计算技术优化性能。此外，我们还将探讨 Sparse Attention 和 DeepSeek-MLA 的创新设计，以及它们如何在大规模语言模型中实现高效的长序列处理。

通过对这些技术的全面分析，本文旨在为读者提供一个系统的视角，理解大模型中 Attention 机制的演进及其在实际应用中的优化策略。我们相信，随着这些技术的不断发展，大模型将在更多领域展现出其强大的潜力，并为人工智能的未来开辟新的可能性。

Attention 技术

MHA(Multi Head Attention)

Attention is all you need.

Single Head

在基础的 MHA，给定一个 Sequence 的表示\(x\in\mathbb{R}^{N\times d}\)，\(N\)是 Sequence 的长度，\(d\)是表征的维度。

构造三个 Tensor, \(W^Q,W^K,W^V\in\mathbb{R}^{d\times d}\), 那么对应的，

\[ Q=xW^Q, K=xW^K, V=xW^V \]

Attention 的结果\(O\)则变成了：

\(Q=g(QK^T)V\)，这边\(g\)代表按行 softmax 操作。

Multi Head

Single Head 往 Multi Head 的推广则只需要做如下的改变就可以了\(W^Q, W^K, W^V\in\mathbb{R}^{d\times\frac{d}{h}\times h},\)\(h\)是 Head 的个数。下面这段提供了 Multi Head 的 Demo 实现的两种方式：

1. 使用 Single Head 的方式，然后拼接起来
2. 直接完整的张量计算方式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49

import torch 

batch_size = 8
seq_len = 16
d = 128 
h = 4
sub_d = d // h

x = torch.randn(size=(batch_size, seq_len, d))

w_qs = [torch.randn(size=(d, sub_d)) for _ in range(h)]
w_ks = [torch.randn(size=(d, sub_d)) for _ in range(h)]
w_vs = [torch.randn(size=(d, sub_d)) for _ in range(h)]


print("========From Single Head START=========")
os = []
for head in range(h):
    q = torch.matmul(x, w_qs[head])
    k = torch.matmul(x, w_ks[head])
    v = torch.matmul(x, w_vs[head])
    att = torch.einsum('bnd,bmd->bnm', q, k) / (sub_d ** 0.5)
    att_sep.append(att)
    att = torch.softmax(att, dim=-1)
    att_sep_soft.append(att)
    o = torch.einsum('bnm,bmd->bnd', att, v)
    os.append(o)
output = torch.cat(os, dim=-1)
print("output shape: ", output.shape)
print("========From Single Head END=========")


print("========VEC START=========")
w_qs_vec = torch.stack(w_qs, dim=-1) # (d, sub_d, h)
w_ks_vec = torch.stack(w_ks, dim=-1)
w_vs_vec = torch.stack(w_vs, dim=-1)
q_vec = torch.einsum('bnd,dsh->bnsh', x, w_qs_vec) # (b, n, d, h)
k_vec = torch.einsum('bnd,dsh->bnsh', x, w_ks_vec)
v_vec = torch.einsum('bnd,dsh->bnsh', x, w_vs_vec)
att = torch.einsum('bnsh,bmsh->bnmh', q_vec, k_vec) / (sub_d ** 0.5)
att = torch.softmax(att, dim=-2)
o = torch.einsum('bnmh,bmsh->bnhs', att, v_vec)
output2 = torch.reshape(o, shape=(batch_size, seq_len, d))
print("output2 shape: ", output2.shape)
print("========VEC END=========")


err = torch.sum(torch.abs(output - output2))
print("err: ", err)

上面的具体计算流程可以归纳为如下

从这边我们可以看到传统的 MHA 针对\(N\)的复杂度\(\sim \mathcal{O}(N^2)\)

总体的 FLOPS 可以总结成如下的表格

Flash Attention

Flash Attention：从算法实现的角度，提出了一种IO高效，无损，快速的MHA计算方法。

对于高性能计算程序来说,往往需要从下面两个角度去考虑程序的性能瓶颈:

1. Compute Bottleneck:

   1. 处理器本身的计算能力成了短板。
   2. 在深度学习场景，例如大矩阵的乘法。

2. Memory and IO Bottleneck

   1. 程序大量和内存，显存或者磁盘交互。
   2. 在深度学习场景，各类 element-wise 操作，各类 Reduce 操作。

那么在 MHA 的计算场景，我们从上面的算法中可以看到，整个流程需要频繁写入 HBM，并且立马从 HBM 再次读取的操作，这样往往拖累整个程序的性能。我们简单看一下 GPU 的结构就可以理解为什么会这样。

从上面的这张图中，我们可以看到整个 GPU 的内存，分成了三个层次:

1. SRAM: 静态随机存储器，价格昂贵，容量小，但是速度快
2. HBM：即常说的显存
3. CPU DRAM：GPU 机器上的 CPU 的内存

我们可以看到 SRAM 的容量只有 20MB，因此往往不能存储下一个完整的超大 Tensor。在对大矩阵进行计算的时候，会将其分块，并分给不同 CUDA 核心去操作，然后通过 HBM 作为相互通信的中间 Buffer。

所以为了提高传统 MHA 的计算效率，需要遵循下面的原则：

尽可能减少和 HBM 的交互次数，对整个计算做融合(Fuse)

因此在 FlashAttention 的算法中，作者提供了两个模块来完成上述的这件事情：

1. Tiling: 使用分块和流式的方式，降低和 HBM 的交互频率
2. Recomputing: 通过存储 Normalization Factor 的方式，省去了存储超大的中间 Attention Matrix，一样可以保证整体梯度回传的正确性

备注：第二点其实在统计物理中非常常用，一个统计分布的配分函数(partition function),即 Normalization Factor,蕴含了该分布的非常多的信息。

Forward Pass

Tiling 的核心思想就是流式分块计算 softmax。

对于给定一个向量\(x\in\mathbb{R}^{N}\),那么起 softmax 可以按照如下的方式计算：

\[ g(x)_i=\frac{e^{x_i-m(x)}}{\sum_j e^{x_j-m(x)}}, m(x)=\max_i\{x_i\} \]

这边同时减去\(m(x)\)只是从数值稳定性的角度出发。

那么流式分块计算是什么意思，假设我们没有办法一下得到整个\(x\),而是分块得到\(x^{(i)}\),其中\(x=(x^{(1)},x^{(2)},\cdots,x^{(n)}),x^{(i)}\in\mathbb{R}^{d_i},\sum_i d_i = N\)，在这种情况下，我们需要找到一个算法，完美地计算出\(g(x)\)

为了方便讨论，我们考虑最简单的情况，\(x\)被分成了两个 block, \(x^{(1)},x^{(2)}\)。那么整个计算可以按照如下的步骤进行：

1. 当我们拿到第一个 block,\(x^{(1)}\)之后，我们先计算这个 block 本身自己的 softmax，同时保留如下的变量:

   1. \[ m(x^{(1)}),l^{(1)}=\sum_j e^{x^{(1)}_j-m(x^{(1)})} \]
   2. \(g(x)\)的前\(d_1\)个元素更新成\(\frac{e^{x^{(1)}_j}-m(x^{(1)})}{l^{(1)}}\)

2. 当我们拿到第二个 block 的时候，我们知道\(x^{(2)}\)对\(x^{(1)}\)那部分 softmax 值的影响，只有两个方面：

   1. 更新最大值\(m(x)=\max(m(x^{(1)}),m(x^{(2)}))\)，这边\(m(x^{(1)})\)之前所有的最大值，一直通过一个变量维护着，\(m(x^{(2)})\)使用新来的数据\(x^{(2)}\)计算
   2. 更新整体的 normalization factor，\(l(x)=e^{m(x^{(1)})-m(x)}l^{(1)}+e^{m(x^{(2)})-m(x)}l^{(2)}\),这边\(l^{(2)}\)可以使用新来的数据\(x^{(2)}\)临时计算
   3. 更新\(g(x)\),\(g(x)\)前\(d_1\)个元素，可以乘上\(\frac{l^{(1)}}{l(x)}e^{m(x^{(1)})-m(x)}\)这个因子去更新，\(g(x)\)后面新的\(d_2\)个元素，则可以直接赋值\(\frac{e^{x^{(2)}-m(x)}}{l(x)}\)

如果理解清楚上面这个简单的\(x\in\mathbb{R}^{N}\)的情况，那么不能扩展到 Attention 真正面临的场景。

1. 从\(\mathbb{R}^{N}\)求 softmax 扩展到任意\(B_r\)行,\(B_c\)列的矩阵\(\mathbb{R}^{B_r\times B_c}\)的按行求 softmax。

   1. \(l,m\)变成一个\(\mathbb{R}^{B_r}\)的向量，\(g\)变成一个\(\mathbb{R}^{B_r\times B_c}\)的矩阵
   2. \(\frac{1}{l(x)}\)的操作变成矩阵求逆操作

2. Atttention Matrix \(QK^T\)按照硬件参数拆分成个个小的 Block

3. 意识到\(V\)的乘法，只和 Block 本身的位置又关系，softmax 本身的流式更新并不影响。

那么在这种情况下面，我们就会知道在 Forward Pass 过程中，FlashAttention 整体的算法可以表述成下图:

Algorithm 1 \(\mathcal{O}(N^2d)\) FLOPs and requires \(\mathcal{O}(N)\) additional memory beyond inputs and output.

FLOPS 主要来自于\(Q_iK_j^T\)的矩阵乘法: \(T_c T_r B_rB_cd \sim \mathcal{O}(N^2d)\)

额外的\(\mathcal{O}(N)\)的显存开销来自于\(l,m\)向量的存储

标准MHA需要和HBM交互\(\mathcal{O}(Nd+N^2)\)次，FlashAttention需要和HBM交互\(\mathcal{O}(N^2d^2M^{-1})\)次

标准 MHA Load \(Q,K,V\)需要\(\mathcal{O}(Nd)\), Attention Matrix 的读写需要\(\mathcal{O}(N^2)\)

FlashAttention 主要 Load \(T_r T_c\)次\(M\)大小量级的 Tensor, 因此\(T_rT_cM \sim \mathcal{O}(N^2d^2M^{-1})\)

在一般情况下面\(M\gg d^2\),因此 FlashAttention 的交互次数会被大幅度压缩。

当然在 LLM 的训练过程中，因为有 Sequence Mask 的存在，往往很多 Attention Block 不需要计算，可以直接认为是一个\(-\infty\)，因此作者证明了，对于 Sparse 结构的 Attention，整体与 HBM 交互次数可以压缩到\(\mathcal{O}(N^2d^2M^{-1}s)\),\(s\)用来很衡量稀疏度。

Backward Pass

在分析完了 Forward Pass 之后，我们需要证明：

维护\(Q,K,V,l\)就足够完成梯度计算，不需要维护中间的Attention Matrix\(P:=g(QK^T)\)

假设梯度回传到\(O\)的时候为\(O^{\prime} \in \mathbb{R}^{N\times d}\), 因为\(O=PV\)，根据链式法则我们可以有:

\(V'_{ij} = \frac{\partial O_{rs}}{\partial V_{ij}} O'_{rs}=P_{ri}\delta_{sj}O'_{rs}=P^T_{ir}O'_{rj}=\sum_r\frac{e^{q_rk_i^T}}{l_r}O'_{rj}\)，即\(V'=P^TO'\).

因为关于\(V\)的梯度，可以仿照从 V 到 O 的过程，使用类似 ForwardPass 的方式进行计算从\(O'\)到\(V'\).

下面我们还需要确认一下\(Q',K'\)是不是也可以比较方便的计算出来。根据\(O=PV\),我们可以类似的得到\(P'=O'V^T\),即我们可以用上面的方式计算出来\(P'\)。

\[ Q'_{i,:}=P'_{is}\frac{\partial P_{is}}{\partial Q_{i:}}=P'_{is}\frac{\partial P_{is}}{\partial S_{it}}\frac{\partial S_{it}}{\partial Q_{i:}}=(\delta_{st}P_{is}-P_{is}P_{it})P'_{is}K_{t:}=\sum_tP_{it}(O'_{i:}V_{j:}-O'_{i:}O_{i:})K_{t:} \]

V2

FlashAttention-V2 的版本在 V1 的基础上做了一些小的优化:

1. 原来在更新O的过程中，每一步都是用老的\(l_{old}\)和新的\(l_{new}\)去做scale，将整个softmax计算结果保持到根据已有数据能算到最准确的结果。但是实际过程中，我们只关心最终结果的准确性，所以中间过程，完全不需要去对O做scale，只需要不断维护最新的\(l\),等到计算结束之后统计做scale
2. 在反向过程中不需要维护完整的\(m\)\(和\)\(l\)，只需要维护\(L=m+\log l\)，但是感觉这边本身的消耗就应该很小？

算法的核心改动就是红色圈出来的两部分。

V3

针对 Hopper GPU，例如 H100 的优化

MQA(Multi Query Attention)

MQA 算法，是对 MHA 一种近似替代，正如上面看到的，在 MHA 中,\(W^Q, W^K, W^V\in\mathbb{R}^{d\times\frac{d}{h}\times h},\)MQA 则是对于\(W^K,W^V\)做了共享的优化，所有的 head 共享相同的参数，变成\(W^K,W^V\in\mathbb{R}^{d\times\frac{d}{h}}\)，以此来降低整个算法的计算复杂度。

GQA(Group Query Attention)

MQA 相比于 MHA 跨度过大，在很多实验中，模型性能下降很明显，因此后续研究者，有提出来一个鉴于 MHA 和 MQA 中的方法，Key 和 Value 从不同 head 之间完全独立和完全共享，选择了一种分组共享的折中方式，以此平衡整个效率和性能。

如果我们\(h\)个 head 分成了\(h_g\)个 Group，我们可以看看整体 FLOPS 减少到多少。首先对比这个表格大模型里的 Attention 机制，我们可以发现，GQA 只是降低了第一步的 FLOPS，

\[ 6Nd^2 \to 2Nd^2+4Nd^2\frac{h_g}{h} \]

不过在 GQA 的文章中，作者提出来一个 Uptraining 的思路，做法很简单，不过开辟了从基于 MHA 训练得到 checkpoints 之后可以自由切换到 MQA 或者 GQA 去进行 finetuning，这样就可以极大可能的保证性能和效率。

做法就是很简单的对于基于 MHA 训练好的 chekpoints 的\(W^K,W^V\)做分组然后去平均。

Sparse Attention

SWA(Sliding Window Attention)

SWA(Sliding Window Attention),该方法是一种“倒退”的想法。我们知道在 CNN 类型的结构中，每个卷积操作都是 Local 信息的融合，需要不断加深网络的深度，才能慢慢捕捉 Global 的信息。

VisionTransformer 的提出，则是想借助 Attention 这种全局的机制，去直接捕获 Global 的信息。

SWA 则是觉得全局 Attention 的机制计算复杂度太高，通过使用 Mask 的方式，将每层的 Attention 限制在局部。

最左边则是正常的 Attention 对应的 Casual Mask，是一个下三角矩阵。中间则是 SWA 限制了 Attention 范围的 Mask 矩阵，最右边则是 SWA 通过深度不断加深之后，整体个模型视野图。

从上面这个图上，我们可以看到 SWA 带来的好处有下面两个点：

1. 从 FLOPS 的角度来说，MHA 计算中很多\(\mathcal{O}(N^2)\)的操作，都不会随着\(N\)再以平方的速度增长。
2. 从 KV Cache 的角度来说，当超过 Window 长度的 KV 就不再需要缓存，大幅度节省了 Cache。

LongFormer

和 SWA 这种借鉴 CNN 思想的 Attention 机制，后面又有一些类似的 Sparse Attention 机制被提出来，总体结构和下面的图类似。

Deepseek-MLA(MultiHead Latent Attention)

MLA 原理

Deepseek-V2 在技术报告中提出来他们的 MLA 技术，主要以此来降低整个 KV Cache 的消耗，同时降低 FLOPS。

从上面我们知道，在原始的 MHA 中，所有的\(Q,K,V\)都来自于上一层的输出表示\(\mathbb{R}^d\), 同时 KV Cache 也需要 Cache \(\mathcal{O}(2Nd)\)。在 MLA 当中，作者采用了类似 LORA 的一些做法，先将\(x\in\mathbb{R}^d\)映射到一个低维的空间，

\[ \begin{gather} c^{KV}=xW_{D}^{KV}, c^Q=xW_{D}^{Q}\\ W_{D}^{KV}\in\mathbb{R}^{d\times d_c^{KV}}, W_{D}^{Q}\in \mathbb{R}^{d\times d_{c}^Q}\\ d_c^{KV},d_c^{Q} \ll d \end{gather} \]

那么这个时候，所有的\(K,V\)都是基于这个共有的\(c\)再升维上去，即：

\[ \begin{gather} k=c^{KV}W_{U}^{K},v=c^{KV}W_{U}^{V},W_{U}^{K}, q=c^QW_{U}^Q\\ W_{U}^V \in \mathbb{R}^{d_c^{KV}\times d}, W_{U}^Q\in \mathbb{R}^{d_c^Q \times d} \end{gather} \]

那我们分别从 FLOPS 和 KV Cache 的角度来看这么做的好处：

1. FLOPS: 对于\(Q,K,V\)的计算的需要的 FLOPS 从 \(6Nd^2\)变成了\(6Ndd_c^{KV}+4Ndd_c^Q\)，而\(d_c^{KV},d_c^{Q} \ll d\)，因此 FLOPS 的得到了大幅度减少。
2. KVCache: 原始 MHA 的 KV Cache 每层需要 \(2Nd\)，MLA的KV Cache只需要\(Nd_c^{KV}\) 因此也大幅度减少了。

那么当然第一个疑问是不是如果我们只缓存\(c^{KV}\)，那么在推理的时候，我们仍然需要实时去计算\(c^{KV}W_U^K, c^{KV}W_U^V\)这两个矩阵乘法，那么 MLA 是不是一种以时间换取空间的策略呢？

当然，这边作者们发现在推理阶段，可以对模型算子做融合，从而避免上面的计算开销。我们以\(qk^T\)的计算为例，看看如何进行融合掉\(W_U^K\)。

\[ qk^T=c^QW_U^Q(W_U^K)^T(c^{KV})^T:=c^QW_U^{QK}(c^{KV})^T \]

因此，我们可以将\(W_U^K\)融入到\(W_U^Q\)当中去。类似的，我们也可以将\(W_U^V\)融合到\(W^O\)，即 Attention 的 FFN 的权重当中去，避免重复计算。

因此在推理阶段，如果考虑算子的融合 FLOPS 会进一步从 \(6Ndd_c^{KV}+4Ndd_c^Q\)压缩到\(2N(dd_c^{KV}+dd_c^Q+d_c^Qd_c^{KV})\)

MLA 的 RoPE

我们知道在 Transformer 架构中，都会加入 Position Encoding 来让模型感知到每个 Token 在时间序列中所处的位置。当然，对于 Position Encoding 的也有很多讨论，不过现在的主流 LLM 都几乎使用 RoPE，这边我们就简单介绍一下 RoPE 的原理。

在笔者看来，RoPE 的特点可以总结成如下的一句话：

RoPE 寻求一种绝对位置的 PE 表示方法，但是又可以保证在计算\(qk^T\)的时候能保持相对位置的一致性。

而\(qk^T\)作为一种内积的计算，满足上面的要求的表示，很显然就是旋转矩阵。因此，对于传统的 MHA，我们则可以乘上一个旋转矩阵,

\[ q_i^R=q_iR(i,\theta),k_j^R=k_jR(j,\theta) \]

这边\(i,j\)代表的是第\(i,j\)个 token。那么，

\[ q_i^R(k_j^R)^T=q_iR(i,\theta)R^T(j,\theta)k_j^T=q_iR(i-j,\theta)k_j^T \]

保证了相对位置的一致性。

那么如果我们直接将这种方式从 MHA 推广到 MLA，会遇到怎么样的问题呢？

\[ q_i^R(k_j^R)^T=c_i^QW_U^QR(i,\theta)R^T(j,\theta)(W_U^K)^T(c_j^{KV})^T \]

从上面的表达式，我们可以看到两个旋转矩阵\(R(i,\theta),R(j,\theta)\)插在了两个\(W\)中间，而这两个旋转矩阵又和绝对的位置\(i,j\)有关系，因此我们就没有办法去做算子融合，这就让 MLA 在推理阶段的性能优势大幅度降低。

那么在 Deepseek 的 V2 中，采取了一种折中的处理方式，值得注意的是，Deepseek 处理方式其实是在一定程度上，破坏了 RoPE 的设计理念的，即：

DeepSeek-MLA-RoPE 是不能保持相对位置的一致性的。

Deepseek 额外引入了\(q^R,k^R\)使用 RoPE 编码位置信息，即:

\[ q^R:=\mathrm{RoPE}\left(c^QW^{QR}\right),k^R:=\mathrm{RoPE}\left(xW^{KR}\right),W^{QR}\in\mathbb{R}^{d_c^Q\times d_R}, W^{KR}\in\mathbb{R}^{d\times d_R} \]

然后使用了暴力拼接的方法得到了\(q,k\):

\[ q\leftarrow [q, q^R], k \leftarrow [k, k^R] \]

参考文献

Ainslie, Joshua, James Lee-Thorp, Michiel de Jong, Yury Zemlyanskiy, Federico Lebrón, and Sumit Sanghai. “GQA: Training Generalized Multi-Query Transformer Models from Multi-Head Checkpoints.” arXiv, December 23, 2023. https://doi.org/10.48550/arXiv.2305.13245.

Beltagy, Iz, Matthew E. Peters, and Arman Cohan. “Longformer: The Long-Document Transformer.” arXiv, December 2, 2020. https://doi.org/10.48550/arXiv.2004.05150.

Dao, Tri, Daniel Y. Fu, Stefano Ermon, Atri Rudra, and Christopher Ré. “FlashAttention: Fast and Memory-Efficient Exact Attention with IO-Awareness.” arXiv, June 23, 2022. https://doi.org/10.48550/arXiv.2205.14135.

Lu, Enzhe, Zhejun Jiang, Jingyuan Liu, Yulun Du, Tao Jiang, Chao Hong, Shaowei Liu, et al. “MOBA: MIXTURE OF BLOCK ATTENTION FOR LONG-CONTEXT LLMS,” n.d.

Shah, Jay, Ganesh Bikshandi, Ying Zhang, Vijay Thakkar, Pradeep Ramani, and Tri Dao. “FlashAttention-3: Fast and Accurate Attention with Asynchrony and Low-Precision.” arXiv, July 12, 2024. https://doi.org/10.48550/arXiv.2407.08608.

Shazeer, Noam. “Fast Transformer Decoding: One Write-Head Is All You Need.” arXiv, November 5, 2019. http://arxiv.org/abs/1911.02150.

Vaswani, Ashish, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. “Attention Is All You Need.” arXiv, August 2, 2023. https://doi.org/10.48550/arXiv.1706.03762.

Yuan, Jingyang, Huazuo Gao, Damai Dai, Junyu Luo, Liang Zhao, Zhengyan Zhang, Zhenda Xie, et al. “Native Sparse Attention: Hardware-Aligned and Natively Trainable Sparse Attention.” arXiv, February 16, 2025. https://doi.org/10.48550/arXiv.2502.11089.

DeepSeek-AI, Aixin Liu, Bei Feng, Bin Wang, Bingxuan Wang, Bo Liu, Chenggang Zhao, et al. “DeepSeek-V2: A Strong, Economical, and Efficient Mixture-of-Experts Language Model.” arXiv, June 19, 2024. http://arxiv.org/abs/2405.04434.

DeepSeek-AI, Aixin Liu, Bei Feng, Bing Xue, Bingxuan Wang, Bochao Wu, Chengda Lu, et al. “DeepSeek-V3 Technical Report.” arXiv, December 27, 2024. https://doi.org/10.48550/arXiv.2412.19437.

Touvron, Hugo, Louis Martin, and Kevin Stone. “Llama 2: Open Foundation and Fine-Tuned Chat Models,” n.d.

Overfitting Testing

在LLM盛行的当下，研究人员出发点都变成了构建一个超级大的模型，去拟合整个人类社会的数据(pre-training)，并通过一些技术向人类演化出来的交互方式去做alignment(post-training)。大家似乎对于过拟合这一个原来机器学习中非常基础的议题不再关注，

Position Based Dynamics

🖌 PBD 这一类物理模拟方法，不同于基于经典力学那套方法，写出动力学方程，然后数值求解 PDE 的范式。而是直接只考虑外力作用，将内力全部归结于不同的约束条件考虑。使用投影的方法去解决外力作用更新之后的约束之间的冲突问题。

Basic info

PBD 相关的物理模拟方法

• Author: anschen

• Reading cycle

  • Start date 2024/05/23
  • End date 2024/06/10

Notations

三维空间的向量 对时间求全导数

\[ \mathbf{x} \]

\[ \dot{x} \]

Position Based Dynamics

在这个章节，我们将主要介绍原始的 PBD 算法的细节和原理，从而构建一个完整的 PBD 这套方法论的大框架。

Model

我们考虑如下建模场景，对于一个物体，有\(N\)个顶点， \(M\)个约束条件，每个顶点的质量，位置，速度分别为\(m_i,\mathbf{x}_i,\mathbf{v}_i\)。这边的约束条件可以按照如下的语言描述：

1. A cardinality \(n_j\)
2. A function \(C_j:\mathbb{R}^{3n_j} \to \mathbb{R}\)
3. A set of indices \(\{i_1, \cdots, i_{n_j}\}, i_k \in [1, \cdots, N]\)
4. A stiffness parameter \(k_j \in [0, 1]\)
5. A type of equality or inequality , i.e. \(C_j(\mathbf{x}_{i_1}, \cdots, \mathbf{x}_{i_{n_j}}) = 0;\)or \(C_j(\mathbf{x}_{i_1}, \cdots, \mathbf{x}_{i_{n_j}}) \ge 0\)

Algorithm

算法总结：

1. 第 5 步，只考虑外力的作用，相对简单。
2. 第 7-8 两步，使用临时的位置去生成新的碰撞约束。
3. 第 9-11 步，则是整个算法的核心，如何通过投影的方式来解决整个约束之间的冲突关系。
4. 第 12-16 步，则通过使用解决完冲突之后的位置，计算速度，以此保证整个过程的一致性。

Constraint Projection

在这个小节，我们将详细介绍如何通过迭代投影的方法去解决约束问题。

首先在经过算法第 5-6 步之后，外力对系统的影响，已经被归纳进入临时的位置\(\mathbf{p}_i\)，那么在约束投影这一步，我们需要我们的系统满足下面两条物理定律（下面的两个方程都忽略了固定时间间隔\(\Delta t\)）：

1. 动量守恒 \(\sum_i m_i \Delta \mathbf{p}_i = 0\)
2. 角动量守恒 \(\sum_i \mathbf{r}_i \times m_i \Delta\mathbf{p}_i = 0\)，其中\(\mathbf{r}_i\)为\(\mathbf{p}_i\)到任意旋转中心的矢量。

对于刚体系统，我们知道任意的约束\(C(\mathbf{p}), \mathbf{p}=(\mathbf{p}_1^T, \cdots, \mathbf{p}_n^T)^T\)都应该满足平移不变性和旋转不变性。因此，这边我们可以借鉴经典力学中的诺特定理[https://en.wikipedia.org/wiki/Noether%27s_theorem]的推导方式，得出如下重要的结论。

NOTE 如果所有顶点的质量\(m_i\)都相等，那么沿着\(\nabla_{\mathbf{p}}C(\mathbf{p})\)更新\(\Delta\mathbf{p}\)，更新后的系统自动满足动量守恒和角动量守恒。

证明分为两部分：动量守恒和角动量守恒

动量守恒

\[ C(\mathbf{p}+\Delta\mathbf{p}) = C(\mathbf{p}) + \nabla_{\mathbf{p}}C(\mathbf{p})\Delta\mathbf{p}+o(\Delta\mathbf{p}^2) \]

如果\(\Delta\mathbf{p}_i = \delta\mathbf{r}, \forall i \in [1, n] \)，即所有顶点平移一个相同的无限小位移\(\delta\mathbf{r}\)，那么平移不变性要求如下的方程成立:

\[ \sum_i \nabla_{\mathbf{p}_i}C(\mathbf{p})\delta\mathbf{r} = 0, \forall \delta\mathbf{r} \]

因此，如果\(m_i = m\)，我们可以得到动量守恒：\(0 = m\sum_i\nabla_{\mathbf{p}_i}C(p) = \sum_im_i \nabla_{\mathbf{p}_i}C(p)\)

角动量守恒

对于旋转不变性，我们可以仿照平移不变性，对\(C(\mathbf{p})\)施加一个类似的无穷小旋转变换，如上图所示。那么这个时候，无穷小旋转变化的对应的偏移矢量\(\delta\mathbf{r}=\delta\mathbf{\Theta} \times \mathbf{r}\)，因此

\[ \sum_i \nabla_{\mathbf{p}_i}C(\mathbf{p}) \delta\mathbf{\Theta} \times \mathbf{r}_i = \delta\mathbf{\Theta} \cdot\sum_i \mathbf{r}_i \times \nabla_{\mathbf{p}_i}C(\mathbf{p}) = 0, \forall \delta\mathbf{\Theta} \]

因此，如果\(m_i = m\)，我们可以得到动量守恒：\(0 = m \sum_i \mathbf{r}_i \times \nabla_{\mathbf{p}_i}C(\mathbf{p}) = \sum_i \mathbf{r}_i \times m_i\nabla_{\mathbf{p}_i}C(\mathbf{p})\)

更新尺度

在我们确定了\(\Delta\mathbf{p}\)的更新方向之后，接下来确定更新的步幅则可以完成整个更新工作，最起码对于质量是一样的情况是这样。考虑 equality 类型的约束，即\(C(\mathbf{p}+\Delta\mathbf{p})=0, \Delta\mathbf{p}=\lambda \nabla_{\mathbf{p}}C(p)\)，我们可以容易求解出来：

\[ \lambda = -\frac{C(\mathbf{p})}{\|\nabla_\mathbf{p}C(\mathbf{p})\|^2} \]

不同质量

在上面的讨论中，我们完整了讨论了，如果质量是相同情况，投影更新规则。但是如果质量各不相同，则上面所说的动量和角动量守恒则会被破坏。以动量守恒为例，平移不变性只能得到\(\sum_i\nabla_{\mathbf{p}_i}C(p)=0\)，那么这个时候如果我们仍然要求\(\sum_im_i \Delta\mathbf{p}_i = 0，\Delta\mathbf{p}_i = \lambda_i \nabla_{\mathbf{p}_i}C(\mathbf{p})\)，一个比较好的处理方式：根据各个顶点的质量\(m_i\)去调整更新的步幅\(\lambda_i\)。

\[ \lambda_i = -\frac{w_i C(\mathbf{p})}{\sum_i w_i \|\nabla_{\mathbf{p_i}}C(\mathbf{p})\|^2},w_i=\frac{1}{m_i} \]

对于这样的更新法则，可以同时保证三件事：

1. 动量守恒
2. 角动量守恒
3. \(C(\mathbf{p}+\Delta\mathbf{p})=0\)

Inequality 约束条件

上面的讨论都是针对了 Equality 类型的约束条件进行的，即\(C(\mathbf{p}+\Delta\mathbf{p})=0\)，对于 Inequality 类型的约束条件，我们则可以比较好的处理，因为\(\nabla_{\mathbf{p}_i}C(\mathbf{p})\)代表\(C(\mathbf{p})\)增加最快的方向，那么如果我们的要求是\(C(\mathbf{p}+\Delta\mathbf{p})\ge 0\)，我们则可以选择只对那些约束条件还小于 0 的约束进行投影更新。

Stiffness

刚度(Stiffness)是指物体抵抗形变的能力。在物理模拟中，刚度决定了物体在受到外力作用时的形变程度。

• 高刚度：如果一个物体或系统的刚度很高，这意味着它不容易形变。在模拟中，高刚度的物体在受到力的作用时会表现出较小的形变，这会使得物体看起来更硬或更紧绷。例如，一根硬木棍在受力时几乎不会弯曲。

• 低刚度：相反，低刚度的物体在受力时会表现出较大的形变，这使得物体看起来更软或更容易弯曲。例如，一根橡皮筋在拉伸时会明显延长。 因此，Stiffness 对于整个模拟的影响，可以直接作用到位移\(\Delta\mathbf{p}\)上面，即\(\Delta\mathbf{p} \leftarrow k\Delta\mathbf{p}\)，但是因为我们整个 Constraint Projection 是一个迭代过程，因此最终的 Stiffness 的影响，变成了非线性的。当然在 PBD 原文中，提出来用\((1-k_s)^{n_{iter}}:=1- k\)求解出来的\(k_s\)作为每个迭代步骤中 Stiffness 产生的影响。不过个人对于这个做法并不是特别的赞同：

  1. \(n_{iter}\)对于不同的模拟和物体，为了保证模拟的精度，是一个可变的量。
  2. \((1-k_s)^{n_{iter}}:=1- k\)的由来，是作者认为误差项为\((1-k)^{n_{iter}}\)，对于该结论持保守观点。

个人观点：原始的 PBD，并不能很好地处理 Stiffness 对于整个模拟的影响，尤其是“物理”意义上的正确处理这个过程。

Collision-碰撞

对于碰撞来说，一般需要考虑的碰撞情形，可以简单分成两大类：

1. 和静态物体的碰撞
2. 和动态物体的碰撞

静态碰撞

静态约束碰撞即检测被模拟对象和周围环境中静态障碍物之间的碰撞情况。这个过程可以简单抽象成如下的检测：

如果和静态障碍发生碰撞，则在模拟过程\(\mathbf{x}_i \to \mathbf{p}_i\)的更新穿过穿过某个法向量为\(\mathbf{n}_s\) 表面，因此该约束的数学表示为: \(C(\mathbf{p}_i) = (\mathbf{p}_i - \mathbf{p}_s) \cdot \mathbf{n}_s \ge 0\)，其中\(\mathbf{p}_s\)为射线\(\mathbf{x}_i \to \mathbf{p}_i\)和面的投影。

动态碰撞

所谓动态碰撞，即被模拟对象之间发生了碰撞。这个过程可以用如下的语言描述:

假设一个物体的顶点\(\mathbf{q} \to \mathbf{p}\)，与另外一个物体发生了碰撞，则另外一个物体存在一个三角面其顶点为\(\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3\)被穿过，其数学表示为: \(C(\mathbf{p},\mathbf{q},\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3) := \pm (\mathbf{p}-\mathbf{q})\cdot ((\mathbf{p}_1-\mathbf{p}_2)\times(\mathbf{p}_1-\mathbf{p_3})) \ge 0\)，这边\(\pm \)取决于最终法向朝向表面的那一侧。

Damping

Damping 是指抵抗运动的力，它消耗掉系统的能量。在物理模拟中，阻尼用于模拟那些减缓运动和振动的效果。

• 高阻尼：高阻尼意味着运动能量会快速消耗，导致物体的运动或振动很快减弱直至停止。这在模拟如防震器或是其他需要快速平息振动的场景时非常有用。
• 低阻尼：低阻尼意味着运动能量消耗得慢，物体的运动或振动会持续较长时间。这可以用来模拟如弹簧振动等持续较长时间的动态效果。

在动力学模拟中，适当的阻尼可以防止系统过度振动或“爆炸”，帮助保持模拟的稳定性和真实性。

比较简单来描述 Damping 的方式，可以用下面的算法:

其中\(\tilde{\mathbf{r}}_i \in \mathbb{R}^{3\times 3}, \tilde{\mathbf{r}}_i \mathbf{v} = \mathbf{r}_i\times\mathbf{v} \)。那么对于两种极限的高低阻尼情况:

1. \(k_{damping} = 0\) 整个过程保证完全没有阻尼，没有阻尼这部分的更新
2. \(k_{damping}=1\) 整个过程，所有顶点的速度，都为质心平移速度加上顶点的转动速度，整个模式表现为刚体模式。

从上面的算法中我们也可以看出来，阻尼只是在限制局部的运动，整体的运动并不受影响。

XPBD

通过对 PBD 的介绍吗，我们可以知道 PBD 的拥有简洁的算法逻辑，但是正如前面我们提到的，PBD 至少存在如下几个方面的问题：

1. PBD 总体上是一个 explicit 的算法，因此非常依赖\(\Delta t\)是一个非常小的值，因此需要较长的步数
2. PBD 对于 stiffness 的处理，缺乏较强的物理依据，同时还是依赖于不确定的参数 \(n_{iter}\)

XPBD 的提出，则是为了缓解上面两个问题带来的模拟效果退化。

XPBD 的总体思路，回归了原始的牛顿力学，即将约束重构成势能，去求解经典力学的方程。在 PBD 中，约束只是作为等式或者不等式的数学条件存在，在每次模拟过程中，通过数学更新去让新的位置满足对应的约束条件，当然这个过程中保证了一些的物理规律，例如动量守恒，角动量守恒等等。但是正如上面所说，还是存在很多没有那么物理的地方。

经典力学离散化

在迭代过程中，位置的信息用如下的序列表示\(\mathbf{x}^1, \cdots, \mathbf{x}^i, \mathbf{x}^{i+1},\cdots\)，那么牛顿第二定律可以很简单地写成如下形式:

\[ \mathbf{M}\frac{\mathbf{x}^{i+1}-2\mathbf{x}^i+\mathbf{x}^{i-1}}{\Delta t^2}=-\nabla U(\mathbf{x}^{i+1}) \]

值得注意的是，在方程的右侧我们使用了\(U(\mathbf{x}^{i+1})\)而非\(U(\mathbf{x}^i)\)，以保证整个过程是一个 implicit 的过程，来缓解对\(\Delta t\)要求很小的严苛条件。当然这边的势能\(U(\mathbf{x})\)是用来描述约束条件的，可以表示成如下形式：

\(U(\mathbf{x}) = \frac{1}{2}C^T(\mathbf{x}) \boldsymbol{\alpha}^{-1} C(\mathbf{x})\)，

其中是一个对角矩阵，用来描述 stffiness 强度的倒数，即所谓的 compliant matrix. 因此，

\[ \mathbf{M}(\mathbf{x}^{i+1}-2\mathbf{x}^i+\mathbf{x}^{i-1})=-\nabla C^T(\mathbf{x}^{i+1})\boldsymbol{\alpha}^{-1}C(\mathbf{x}^{i+1})\Delta t^2 := \nabla C^T(\mathbf{x}^{i+1}) \boldsymbol{\lambda}^{i+1} \]

在上面的方程中\(\boldsymbol{\lambda}:=-\tilde{\boldsymbol{\alpha}}^{-1}C(\mathbf{x})\)，\(\mathbf{x}^{i+1}, \boldsymbol{\lambda}^{i+1}\)为需要求解的变量。同时我们知道\(\tilde{\mathbf{x}}^i:=2\mathbf{x}^{i}-\mathbf{x}^{i-1}=\mathbf{x}^i+\mathbf{v}^i\Delta t\)。最终需要求解的方程组为如下：

\[ \mathbf{M}(\mathbf{x}-\tilde{\mathbf{x}})-\nabla C^T(\mathbf{x})\boldsymbol{\lambda} = 0 \\ C(\mathbf{x})+\tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda} = 0 \]

其中\(\tilde{\boldsymbol{\alpha}} = \frac{\boldsymbol{\alpha}}{\Delta t^2}\)，同时省略了上标。对于如何求解上面的方程有很多的方法， 这边采用最简单原始的方法，即牛顿迭代法[https://zh.wikipedia.org/wiki/%E7%89%9B%E9%A1%BF%E6%B3%95], 因此我们直接线性化上面的方程：

\[ \left[\begin{array}{cc}\mathbf{M} - (\nabla\nabla^T) C^T(\mathbf{x}^i) \boldsymbol{\lambda}^i & -\nabla C^T(\mathbf{x}^i) \\ \nabla C(\mathbf{x}^i) & \tilde{\boldsymbol{\alpha}}\end{array}\right]\left[\begin{array}{c}\Delta\mathbf{x}\\ \Delta\boldsymbol{\lambda}\end{array}\right]=-\left[\begin{array}{c}\mathbf{M}(\mathbf{x}^i-\tilde{\mathbf{x}})-\nabla C^T(\mathbf{x}^i)\boldsymbol{\lambda}^i \\ C(\mathbf{x}^i)+\tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda}^i \end{array}\right] \]

XPBD 在处理上面的方程，为了更加容易求解，引入了如下两个假设： - \[(\nabla\nabla^T) C^T(\mathbf{x}^i) \boldsymbol{\lambda}^i \approx 0\] - \[\mathbf{M}(\mathbf{x}^i-\tilde{\mathbf{x}})-\nabla C^T(\mathbf{x}^i)\boldsymbol{\lambda}^i \approx 0\]

对于被忽略的两项，我们可以知道误差都在\(o(\Delta t^2)\)这个量级上。同时整个方程被简化成：

\[ \left[\begin{array}{cc}M & -\nabla C^T(\mathbf{x}^i) \\ \nabla C(\mathbf{x}^i) & \tilde{\boldsymbol{\alpha}}\end{array}\right]\left[\begin{array}{c}\Delta\mathbf{x}\\ \Delta\boldsymbol{\lambda}\end{array}\right]=-\left[\begin{array}{c} 0 \\ C(\mathbf{x}^i)+\tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda}^i \end{array}\right] \]

使用 Shur Complement[https://en.wikipedia.org/wiki/Schur_complement] 我们可以轻松解出上面的方程:

\[ \Delta\mathbf{x}=\mathbf{M}^{-1}\nabla C^T(\mathbf{x}^i)\Delta\boldsymbol{\lambda}\\ \Delta\boldsymbol{\lambda}=(\nabla C(\mathbf{x}^i)\mathbf{M}^{-1}\nabla C^T(\mathbf{x}^i)+\tilde{\boldsymbol{\alpha}})^{-1}(C(\mathbf{x}^i+\tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda}^i)) \]

Connection with PBD

上面直接给出来了\(\Delta\mathbf{x},\Delta\boldsymbol{\lambda}\)的解析更新公式。但是实际使用时，直接对矩阵求逆不是一个很好的操作，更多的时候会采用一些其他的算法来求解线性方程组，例如经典的 Gauss-Seidel 算法。

然而，在原始论文里，作者直接得出了如下的结论:

\[ \Delta\boldsymbol{\lambda}_j=\frac{-C_j(\mathbf{x}^i)-\tilde{\boldsymbol{\alpha}}_j\boldsymbol{\lambda}_j^i}{\nabla C_j(\mathbf{x}^i)M^{-1}\nabla C_j^T(\mathbf{x}^i)+\tilde{\boldsymbol{\alpha}}_j}` \]

这边是十分不理解为什么可以讲不同约束之间的相互影响给独立处理。

当然这边如果处理完毕之后，我们会发现在\(\tilde{\boldsymbol{\alpha}}_j=0\)的情况下面，整个过程退化到 PBD 的场景了。

Damping

XPBD 中使用了 Rayleigh dissipation potential 来处理能量耗散的问题：

\[ D(\mathbf{x}):=\frac{1}{2}\dot{C}^T(\mathbf{x})\boldsymbol{\beta}\dot{C}(\mathbf{x})=\frac{1}{2}\mathbf{v}^T\nabla C^T(\mathbf{x})\boldsymbol{\beta}\nabla C(\mathbf{x})\mathbf{v} \]

\(\boldsymbol{\beta}\)是一个对角矩阵，用来衡量各个约束对于能量耗散的贡献程度。这种形式的耗散是能产生的力的作用:

\[ \mathbf{F}=-\nabla_{\mathbf{v}}D = -\nabla C^T(\mathbf{x})\boldsymbol{\beta}\nabla C(\mathbf{x})\mathbf{v} \]

如果这个时候，我们定义一个新的\(\boldsymbol{\mu}:=-\tilde{\boldsymbol{\beta}}\nabla C(\mathbf{x})\mathbf{v}\)，其中\(\tilde{\boldsymbol{\beta}}=\boldsymbol{\beta}\Delta t^2\)，我们就会发现整个过程和上面的推导没有本质的区别，只需要如下的变化：

\[ \mathbf{M}(\mathbf{x}-\tilde{\mathbf{x}})-\nabla C^T(\mathbf{x})(\boldsymbol{\lambda}+\boldsymbol{\mu}) = 0 \\ C(\mathbf{x})+\tilde{\boldsymbol{\alpha}}\tilde{\boldsymbol{\beta}}\nabla C(\mathbf{x})\mathbf{v}+\tilde{\boldsymbol{\alpha}}(\boldsymbol{\lambda}+\boldsymbol{\mu}) = 0 \]

如果我们将\(\boldsymbol{\lambda}+\boldsymbol{\mu} \to \lambda; \mathbf{v}=\frac{\mathbf{x}-\mathbf{x}_{n}}{\Delta t}\)这样看待，去线性化上面的方程同时使用相同的近似条件，我们会得到:

\[ \left[\begin{array}{cc}\mathbf{M}  & -\nabla C^T(\mathbf{x}^i) \\ \nabla C(\mathbf{x}^i)+\tilde{\boldsymbol{\alpha}}\tilde{\boldsymbol{\beta}}\nabla C(\mathbf{x})\frac{1}{\Delta t}  & \tilde{\boldsymbol{\alpha}}\end{array}\right]\left[\begin{array}{c}\Delta\mathbf{x}\\ \Delta\boldsymbol{\lambda}\end{array}\right]=-\left[\begin{array}{c}0\\ C(\mathbf{x}^i)+\tilde{\boldsymbol{\alpha}}\boldsymbol{\lambda}^i+\tilde{\boldsymbol{\alpha}}\tilde{\boldsymbol{\beta}}\nabla C(\mathbf{x}^i)\frac{\mathbf{x}^i-\mathbf{x}_n}{\Delta t} \end{array}\right] \]

当然这个时候，如果继续使用作者上面奇怪的处理方式，可以得到\(\Delta\boldsymbol{\lambda}\)每个分量上的更细法则:

\[ \Delta\boldsymbol{\lambda}_j=\frac{-C_j(\mathbf{x}^i)-\tilde{\boldsymbol{\alpha}}_j\boldsymbol{\lambda}_j^i-\boldsymbol{\gamma}_j\nabla C(\mathbf{x}^i)(\mathbf{x}^i-\mathbf{x}_n)}{(1+\boldsymbol{\gamma}_j)\nabla C_j(\mathbf{x}^i)M^{-1}\nabla C_j^T(\mathbf{x}^i)+\tilde{\boldsymbol{\alpha}}_j} \]

这边\(\boldsymbol{\gamma}_j=\frac{\tilde{\boldsymbol{\alpha}_j}\tilde{\boldsymbol{\beta}_j}}{\Delta t}\)

Constraints

这个章节，我们将介绍几种常见的约束条件，以及其常见的更新条件，即主要给出约束条件关于坐标的梯度表达式。

Length Constraint(Stretch Constraint)

长度的约束，用来保证相邻两点之间的长度和初始长度接近，因此可以描述成如下的方式:

\[ C(\mathbf{x}_i,\mathbf{x}_j):=\|\mathbf{x}_i-\mathbf{x}_j\|-l_{ij} = 0 \]

\[\nabla_{\mathbf{x}_i}C(\mathbf{x}_i,\mathbf{x}_j)=\frac{\mathbf{x}_i-\mathbf{x}_j}{\|\mathbf{x}_i-\mathbf{x}_j\|}\]

Bend Constraint

Bend Constaint 用来两个相邻三角形之间的夹角接近初始的夹角。

假设两个三角形公共的边为\(\mathbf{x}_1\to\mathbf{x}_2\)，其中第一个三角形另外一个顶点为\(\mathbf{x}_3\)，第二个三角形的另外一个顶点为\(\mathbf{x}_4\),那么可以定义我们的约束条件为下面的形式:

\[ C(\mathbf{x}_1,\mathbf{x}_2,\mathbf{x}_3,\mathbf{x}_4):=\arccos\left(\frac{\mathbf{p}_2\times\mathbf{p}_3}{\|\mathbf{p}_2\times\mathbf{p}_3\|}\cdot\frac{\mathbf{p}_2\times\mathbf{p}_4}{\|\mathbf{p}_2\times\mathbf{p}_4\|}\right)-\theta:=\arccos I - \theta \]

这边\(\mathbf{p}_2:=\mathbf{x}_2-\mathbf{x}_1,\mathbf{p}_3:=\mathbf{x}_3-\mathbf{x}_1,\mathbf{p}_4:=\mathbf{x}_4-\mathbf{x}_1\).

对于这个复杂表达式子，我们要想获得其对应的梯度，核心则需要计算下面的子表达式的梯度，其他的则使用链式法则就可以，(这边为了方便省去了粗体)。

\[ \nabla_x \frac{(x\times y)\cdot z}{\|x\times y\|}=\nabla_x \frac{(y\times z)\cdot x}{\|x\times y\|}=\left(x\cdot(y\times z)\right)\left(\nabla_x\frac{1}{\|x\times y\|}\right)+\frac{y\times z}{\|x\times y\|} \]

\(\nabla_x\|x\times y\|=\frac{\nabla_x\left((x\times y)\cdot(x\times y)\right)}{2\|x\times y\|}=\frac{y\times(x\times y)}{\|x\times y\|}:=y\times n\) 这边我们定义\(n:=\frac{x\times y}{\|x\times y\|}\)，那么

\[ \nabla_x \frac{(x\times y)\cdot z}{\|x\times y\|}=\frac{(n\cdot z)(n\times y)+y\times z}{\|x\times y\|} \]

那么我们定义:

\[ \mathbf{n}_1:=\frac{\mathbf{p}_2\times\mathbf{p}_3}{\|\mathbf{p}_2\times\mathbf{p}_3\|};\mathbf{n}_2:=\frac{\mathbf{p}_2\times\mathbf{p}_4}{\|\mathbf{p}_2\times\mathbf{p}_4\|} \]

1. \[\nabla C = -\frac{1}{\sqrt{1-I^2}}\nabla I\]
2. \[\nabla_{\mathbf{x}_3}I=-\frac{I\mathbf{n}_1\times\mathbf{p}_2+\mathbf{p}_2\times\mathbf{n}_2}{\|\mathbf{p}_2\times\mathbf{p}_3\|}\]
3. \[\nabla_{\mathbf{x}_4}I=-\frac{I\mathbf{n}_2\times\mathbf{p}_2+\mathbf{p}_2\times\mathbf{n}_1}{\|\mathbf{p}_2\times\mathbf{p}_4\|}\]
4. \[\nabla_{\mathbf{x}_2}I=\frac{I\mathbf{n}_1\times\mathbf{p}_3+\mathbf{p}_3\times\mathbf{n}_2}{\|\mathbf{p}_2\times\mathbf{p}_3\|}+\frac{I\mathbf{n}_2\times\mathbf{p}_4+\mathbf{p}_4\times\mathbf{n}_1}{\|\mathbf{p}_2\times\mathbf{p}_4\|}\]
5. \[\nabla_{\mathbf{x}_1}I=-\nabla_{\mathbf{x}_2}I-\nabla_{\mathbf{x}_3}I-\nabla_{\mathbf{x}_4}I\]

时间序列阅读笔记

引言

通用的时间序列建模方法

操作步骤：

• 画出主要feature的图，并查看
  • 是否有一个趋势(Trend)
  • 是否有周期的成分(Period)
  • 是否有突变(Sharp change)
  • 是否有异常值(Outlier)
• 移除趋势和周期项，获得一个接近静态的时间序列
• 建模拟合剩余的残差项
• 预测残差，并加上趋势和周期项，获得最终的预测值

静态模型

Def: Weakly Stationary A time series \({X_t}\) is called weakly stationary, if 1. \(\mu_X(t)\) is independent of time \(t\) 2. \(\gamma_X(t+h,t)\) is independent of time \(t\) for each \(h\), where \(\gamma_X(r, s) = Cov(X_r, X_s)\)

对于一个静态序列来说，\(\gamma_X(h):=\gamma_X(0,h)=\gamma_X(t,t+h)\), 此时，我们称\(\gamma_X(h)\) 为自协方差(auto covariance function).

Def: AutoCorrelation Function(自相关函数) \(\rho_X(h) = \frac{\gamma_X(h)}{\gamma_X(0)}\)

估计趋势和周期

假设一个时间序列由三个部分构成： \[ X_t = m_t + s_t + Y_t \] 其中\(m_t, s_t\)分别表示确定性的趋势项，周期项。\(Y_t\)剩余的一个带有随机性的静态序列。这节我们简单介绍一下一些简单的方法估计\(m_t, s_t\)

只有趋势项的时候

如果我们的模型只有趋势项，如下 \[ X_t = m_t + Y_t \] 那么常见的用来估计趋势项的方法：

• 移动平均平滑
  • 等权平均
  • 指数平均
• 过滤高频成分
  • 傅立叶变换
  • 其他的谱分析方法
• 参数化模型并拟合
  • 多项式拟合
  • 机器学习，深度学习方法

当然，我们也可以反其道而行，除了我们通过上述的方法，去获取一个预测趋势的模型，我们也可以试图消除趋势，直接获得静态序列\(Y_t\)。 消除趋势最常见的方法则是： - 差分

Def: Lag operator \[ BX_t = X_{t-1} \]

Def: Difference operator \[\nabla = 1 - B\]

因此如果\(m_t\)是一个\(k\)阶的多项式，那么\(\nabla^k X_t\) 则足以消除趋势。

同时有趋势项和周期项

上面介绍完了，在只有趋势项的是，我们如何去处理一个时间序列，这边我们考虑的模型更加完整一点，即趋势项和周期项同时存在。 假设周期为\(d\), 一般我们可以遵从下面的步骤：

• 使用当前时间点附近\(d\)个时间的数据做移动平均，大概估计出一个\(\hat{m}_t\)
• 对于周期内的任意一个时间点\(k, 1 \le k \le d\), 计算\(\{x_{k+jd}-\hat{m}_{k+jd}\}, 1 \le k+jd \le n\)的平均，计作\(w_t\), 为了保证周期震荡总体围绕0点震荡，我们可以使用\(\hat{s}_t = w_t - \frac{1}{d}\sum_{k=1}^d w_k\) 作为一个不错的周期项估计
• 现在我们获得了去周期的数据\(d_t = x_t -\hat{s}_t\), 问题则变成了，上面只有趋势项的问题了，我们就可以上一小段列举的方法解决这个问题了。

Stationary Process

Def: Nonnegative definite A real-valued function \(\kappa\) defined on integers is non-negative definite, if \[ \sum_{i,j}a_i \kappa(i-j) a_j \ge 0 \] for all positive integer \(i,j\) and vector \(a=(a_1, a_2, \cdots)^T\)

Portfolio Optimization

简介

投资组合的优化是量化投资中非常重要的一个环节。假设我们有\(N\)支资产，每个资产收集了\(T\)时间的数据，用 \(X \in \mathrm{R}^{N \times T}\) 表示整个风险收益样本。在量化投资中最核心的两个问题就是：

• 给定历史轨迹的条件下，预测未来时间的收益 \(\mathbb{E}[R_{t+1} |\mathcal{F}_{t}]\)。
• 给定历史轨迹的条件下，预测未来时间的各个资产之间的correlation,例如:\(\Sigma_{t+1} | \mathcal{F}_{t}\)。

如果能够很好得完成第一个任务，则就可以追求到人们梦寐以求的超额收益，于此我们也可以想象第一个问题的难度是如何地高，也超过了本篇Blog所想要讨论的范畴。

我们这边主要会将讨论的内容局限在第二个问题上，假设存在一个上帝已经帮助我们解决了第一个问题，即我们已经拥有了一个相对不错的对未来收益估计的算法，使用记号\(\hat{\mu}_{t+1}\) 代表我们对于 $ [R_{t+1} |_{t}]$ 的估计。那么在投资组合这一步，我们往往关注的是下面这个优化问题， \[ \max_{w} \hat{\mu}_{t+1}^T w, \\ f_i(\hat{\mu}_{t+1}, \hat{\Sigma}_{t+1}, w) \le 0, i = 1,2,\cdots, k \\ g_j(\hat{\mu}_{t+1}, \hat{\Sigma}_{t+1}, w) = 0, j = 1,2,\cdots, l \]

用通俗的语言描述，就是在风险可控并且符合法规(例如只允许做多)的条件下，去最大化预期收益。想要控制的风险和其他相关的约束都是用\(f, g\)两种形式描述。那么为了很好地解决这个问题，很明显我们有两个问题需要解决：

• 准确估计 \(\hat{\Sigma}_{t+1}\)， 即上面的第二个任务
• 开发一个高效的优化算法，解决上面这个优化问题。

Correlation估计

在这个部分我们将会介绍一下，关于Correlation估计的一些基本的算法。针对这个问题的建模假设，也是被分成了两部分。第一，假设Correlation是静态的；第二，假设Correlation是动态。第一种情况下面，我们拥有一些比较良好的数学保障，但是和现实背离得稍微远一点。第二种情况，则是和现实更近，但是解决该问题的难度也会高上很多。

静态Correlation估计

对于静态的Correlation的估计，首先我们会想到一个最基本的方法，就是用样本的协方差矩阵估计全体的协方差矩阵。所谓的样本协方差矩阵，可以写成下面的形式: \[ S = \frac{1}{T}(XX^T - \frac{1}{T}X\mathbf{1}\mathbf{1}^TX) = \frac{1}{T}X(\mathbf{I} - \frac{1}{T}\mathbf{1}\mathbf{1}^T)X^T \] 从这个表达式，我们就可以看出来一个较为严重的问题, \(S\)的rank不会高于\(\mathbf{I} - \frac{1}{T}\mathbf{1}\mathbf{1}^T\), 即不会高于\(T\)。那么当你所采集到的数据时间窗口\(T\)和\(N\)在同一个量级，甚至更低的时候，该估计矩阵很有可能是singluar的，这会极大的影响到后续优化问题的求解。该问题在现实投资中还是比较明显的，尤其是当你的投资频率较低时，该问题会越发显著。那么主流的解决上面这个问题的方法也分成两大主流的流派，当然这两个流派也不完全独立，也是可以混合使用。

• Shrinkage方法
• Multi-factor估计方法

Shrinkage

Shrinkage的想法非常的简单，就是既然\(S\)本身作为\(\Sigma\)的估计在上述 \(N,T\) 比较接近的时候有一些问题，那么我们可以再额外另外一个估计 \(F\), 使用这两个估计的组合，例如\((1-\delta)S + \delta F\), 作为\(\Sigma\)的估计呢，这边 \(\delta\) 是一个描述shrinkage程度的参数。在具体讨论这个\(F\)之前，我们可以思考一下引入这个\(F\), 我们需要解决什么，即期望\(F\)能够拥有哪些性质。下面是我自己个人的一些想法：

首先，\(S\)矩阵对于\(\Sigma\)的估计无偏的，但是如果样本量不够的情况下面，其会具有很高的 variance，即统计学习里，常说的 bias-variance trade off。那么很显然，我们是希望引入的\(F\)对于 biasness 的容忍度可以高一点，但是希望其有更低的 variance。除此之外，我们期望\(F\)的引入可以缓解 singularity 的问题。

虽然\(F\)确实构造的自由度很高，你可以选择任何的方法去构造一个合适的估计算子，但是就个人理解来说，上面这段话应当做一个基本的参考准则。我想这也是Ledoit,Wolf 在他们最初的两个工作里选择 \(F\)的原则 ^1,2

• 认为任何两个股票之间的correlation都是独立的，采取全局平均的相关程度作为\(F\)
• 采用单因子模型去估计协方差矩阵

最优的 \(\delta\)

这边我们会介绍一下如何获取最优的 \(\delta\)。根据上面的介绍，我们知道找到最优的 \(\delta^*\) 等价于求解下面这个问题 \[ \min_{\delta} \mathbb{E}\|(1-\delta) S + \delta F - \Sigma\|_2^2 \] 该问题的解为 \[ \delta^* = \frac{\sum_{ij}\mathbb{Var}[s_{ij}]- \mathbb{Cov}[f_{ij}, s_{ij}]}{\sum_{ij}\mathbb{Var}[s_{ij}-f_{ij}] + (\phi_{ij}-\sigma_{ij})^2} \] 其中 \(\Phi\) 为构造出来的估计算子 \(F\)的真实值。具体可以参考附录1

很显然这个式子有两个问题：

• 依赖\(F, S\)的真实统计，方差和协方差
• 依赖真实未知的 \(\Sigma, (\sigma_{ij})\)

对于第一个问题，我们考虑在\(T\)比较大的时候，真实统计和样本估计之间的误差是什么样子的。如果误差较小，那么我们则可以对于 \(\delta^*\) 构造一个渐进有效的表达式。 首先对于第一个问题来说，中心极限定理的直觉告诉我们，无论是\(\mathbb{Var}[s_{ij}],\mathbb{Cov}[f_{ij}, s_{ij}], \mathbb{Var}[s_{ij}-f_{ij}]\) 其中哪一个，他跟随\(T\)的scaling应该都是在 \(\frac{1}{T}\) 这量级上，也是就说 \(\delta^* \propto \frac{\alpha}{T}\), 那么一个最合理，也是最符合直觉的猜测，就是我们可以认为： \[ \delta^* = \frac{1}{T}\frac{\pi - \rho}{\gamma} + o(\frac{1}{T^2}) \] 这边, \[ \pi = \lim_{T \to \infty} \sum_{ij} \mathbb{Var}[\sqrt{T}s_{ij}],\\ \rho = \lim_{T \to \infty} \sum_{ij} \mathbb{Cov}[\sqrt{T}f_{ij}, \sqrt{T}s_{ij}],\\ \gamma = \sum_{ij}(\phi_{ij} - \sigma_{ij})^2 \] 严格的数学分析可以参考附录1。

这边我们对 \(\delta^*\) 做一些简要的说明，因为 \(\delta^* \sim \frac{1}{T}\), 那么我们可以知道 \(T \to \infty\) 的时候，\(\delta^* \to 0\), 整个估计收敛到 \(S\), 这也是和我们统计学上所学到的知识是吻合的。当我们拥有足够多的样本是，使用样本协方差矩阵作为真实协方差矩阵的估计是一个有效的估计。

Multi-factor

动态Correlation估计

优化算法

附录

最优的 \(\delta\)

\[ \|(1-\delta) S + \delta F - \Sigma\|_2^2 = \sum_{i,j=1}^N \left((1-\delta) s_{ij} + \delta f_{ij} - \sigma_{ij}\right)^2 \\ = \sum_{i,j=1}^N \delta^2 \left(s_{ij}^2 + f_{ij}^2 - 2s_{ij}f_{ij}\right) - 2\delta \left(s_{ij}^2 + f_{ij}\sigma_{ij}-s_{ij}\sigma_{ij}-s_{ij}f_{ij}\right) + C \] 因为上式为一个开口向上的二次函数，最后一项跟极值点没有太大关系。同时我们可以发现 \[ \delta^* = \frac{\sum_{ij}\mathbb{E}\left[s_{ij}^2 + f_{ij}\sigma_{ij} -s_{ij}\sigma_{ij}-s_{ij}f_{ij}\right]}{\sum_{ij}\mathbb{E}[s_{ij}^2 + f_{ij}^2 - 2s_{ij}f_{ij}]} \] 因为\(S\)是一个无偏的估计,因此 \(\mathbb{E}[S] = \Sigma\), 同时我们假设构造出来的估计 \(F\)在统计意义上也会收敛到 \(\Phi\)。那么我们可以知道， \[ \delta^* = \frac{\sum_{ij}\mathbb{Var}[s_{ij}]- \mathbb{Cov}[f_{ij}, s_{ij}]}{\sum_{ij}\mathbb{Var}[s_{ij}-f_{ij}] + (\phi_{ij}-\sigma_{ij})^2} \] 下面我们将以 \(\sum_{ij}\mathbb{Var}[\sqrt{T}s_{ij}]\) 为例证明其渐近收敛到 \(\pi\)，以此说明三个二次项的scaling均为\(\frac{1}{T}\)。 显然\(s_{ij}=\frac{1}{T}\sum_t x_{it}x_{jt} - \frac{1}{T}\sum_t x_{it}\frac{1}{T}\sum_{t}x_{jt}\)。对于第二项，中心极限定理告诉我们 \(\sqrt{T}(\frac{1}{T}\sum_{t}x_{it} - \mathbb{E}[x_i]) \to \mathcal{N}(0, \mathbb{Var}[x_i])\), 假设\(X\)中的所有的随机变量都是有限的4阶以下的moment，那么对于\(\sqrt{T}s_{ij}\)中的第二项在\(T\)很大的时候，则不会有实质的贡献，因此我们只需要把重点放在第一项。对于第一项，我们还知道\(\sigma_{ij}=\mathbb{E}[x_i x_j]-\mathbb{E}[x_i]\mathbb{E}[x_j]\)是一个确定性的数字，对比两项的差距，我们发现只要分析\(\frac{1}{T}\sum_{t}x_{it}x_{ij} - \mathbb{E}[x_i]\mathbb{E}[x_j]\)的性质。同样，中心极限定理告诉我们\(\sqrt{T}(\frac{1}{T}\sum_{t}x_{it}x_{ij} - \mathbb{E}[x_i]\mathbb{E}[x_j] - \sigma_{ij}) \to \mathcal{N}(0, \rho_{ij})\)

参考文献

Successful Algorithm Trading 读后感

概览

算法交易和主观交易的对比

优点:

• 提高了投研效率
• 没有主观输入，受情绪影响较小
• 评价指标更加多维
• 更高的交易频率，人工操作不来

缺点:

• 需要一定的资金量
• 需要较高的科学素养和编程水平

散户和机构的对比

散户的优势:

• 小资金量更加灵活
• 对冲基金之间人员流动，会导致策略趋同，造成交易拥挤，散户的交易系统和这些关联较小
• 散户的交易对市场价格的冲击可以忽略

散户的劣势:

• 对杠杆的使用有更多的限制
• 散户的交易在券商的排队往往优先级较低，这会造成实盘和回测有较大差距
• 信息数据滞后

散户对于风险管理相比于基金来说，会更加自由，完全取决于个人的风险偏好。这有好有坏，往往这也会让散户忽略风险管理，把绝大部分时间放在构建交易策略上。

散户不需要关注交易曲线，不需要和投资者打交道，不需要定期公布事项满足监管，不需要和同行比，不需要和被动投资比。 散户只需要关注的就是绝对收益。

回测

为什么需要回测

回测可以帮助做到下面几点：

• 筛选：可以快速帮我筛选策略
• 建模：可以快速验证对市场的想法
• 优化：可以精细打磨自己的策略

回测过程中常见的偏差(Bias)

我们通过回测慢慢迭代自己的策略的时候，往往也引入了很多的偏差(Bias)，常见的有：

• 优化偏差(Optimisation Bias)
  • 过拟合数据，包括模型，策略超参数等等
• 未来信息(Look ahead Bias)
• 幸存偏差(Survivorship Bias)
  • 每年有大量的股票退市，但是往往数据集里面只会包含存活至今的股票数据
• 主观偏差(Cognitive Bias)
  • 面对回测曲线里面最大回撤在25%，可能较容易接受。但是实盘面对如此的回测往往难以忍受，因此很多策略终止在低点。

交易损耗

交易过程中带来的损耗包括但不限于：

• 手续费
  • 印花税
  • 佣金
  • ...
• 滑点(Slippage)
  • 交易信号对应的价格和实际价格之间的差距
  • 主要因素包括：波动性，延迟，交易频率
• 市场冲击

交易风格

在确定自己的交易风格之前，我们需要从下面几个方面出发： - 你的性格是什么样子的 - 你所允许的交易时间是怎么样的

时间序列分析

这一节关于时间序列分析的相关算法，作者讲得非常简洁，不值得记录。

预测

这边作者只是简单减少了一些常见的机器学习算法，例如逻辑回归，线性回归，支持向量机，随机森林等等，以及机器学习当中常见的对于结果判别的方法。

业绩归因 && 风险管理

从各个level的业绩归因，可以更好地回溯整个交易系统的问题。一般可以从如下几个方面看待整个业绩的表现:

• 策略
  • 策略在回测中的模拟
  • 策略在模拟盘中的表现(模拟盘交易误差几乎可以忽略)
• 交易执行
  • 衡量策略提供的交易信号和现实交易执行产生的误差
• 投资组合管理
  • 是否存在更加合理的投资组合管理，可以降低交易的损耗

那么何为上述所说的"业绩"呢，一般我们需要从下面几个角度去考虑:

• 收益
  • 总收益: 计算的时候需要考虑杠杆，做空等等因素
  • 复合年化收益
• 风险
  • 回撤(DrawDown)
  • 收益的波动率
• 风险收益
  • Sharpe Ratio: \(\sqrt{N}\frac{\mathbb{E}[R_{\alpha}- R_{\beta}]}{\sqrt{\mathbb{Var}[R_{\alpha} - R_{\beta}]}}\), 这边\(N\)为一年内交易日数目，如果你的收益是按照天计算的话。\(R_{\alpha} - R_{\beta}\) 代表超额收益。其中Benchmark可以选择为对应的ETF.
  • Sortino Ratio: \(\sqrt{N}\frac{\mathbb{E}[R_{\alpha}- R_{\beta}]}{\sqrt{\mathbb{Var}[R_{\alpha} - R_{\beta}]}_d}\). 这边的下标\(d\)代表我们只考虑，当超额收益为负的时候波动率，因为亏损的波动，更加容易交易情绪的波动。
  • Calmar Ratio: \(\sqrt{N}\frac{\mathbb{E}[R_{\alpha}- R_{\beta}]}{\max \mathbf{drawdown}}\).
• 交易指标
  • PnL
  • Average Period PnL
  • Maximum Period Profit
  • Maximum Period Loss
  • Average Period Profit
  • Average Period Loss
  • Wining Period
  • Losing Period
  • Percentage Win/Loss Periods

Prepare Work Env for Linux

Prepare your work environment for linux system

线性回归漫谈

线性回归是一项假设因果线性的数据拟合方法，如今也是被大家广泛使用。但是其中对应underlying的条件假设都被大家忽略了，这篇文章尝试详细整理线性回归的方方面面，方便之后回顾和参考。大部分的资料都可以在参考文献中找到。

Time Series Analysis

This is unfinised time series analysis note, which is wished to be done around 2023 March

Diffusion map method

本篇介绍一种半监督的算法，以及其在理论物理里面的应用

非线性因果分析技术

在基础的统计学中，我们学过用协方差矩阵来表征两个高维随机变量\(X=(X_1,X_2,\cdots,X_n)^T, Y=(Y_1,Y_2,\cdots, Y_m)^T\)之间的相关性。但是，我们知道相关性并不意味着独立性，更加不能准确描述因果性。这篇Blog，我们将持续更新一些因果性分析工具。