Aller au boulot!

Welcome to my notes space, I really appreciate if you find something helpful here.

Host on excel.kuizhiqing.cn.

chronicles

chronicles

Feb 13

Summarizey the story of invention of Blue LED.

The invention of the blue LED was a significant breakthrough in lighting technology, for which Isamu Akasaki, Hiroshi Amano, and Shuji Nakamura were awarded the Nobel Prize in Physics in 2014. Shuji Nakamura, a Japanese-born American electronic engineer, played a key role in this innovation. He worked at the Nichia Corporation, where he developed a bright blue LED using indium gallium nitride (InGaN) and a low-temperature heating method, significantly increasing its efficiency. This blue LED was nearly 100 times brighter than existing ones at the time, enabling its use in various applications. By combining the bright blue LED with fluorescent material, Nakamura enabled the creation of white light. Additionally, this development led to the production of power-efficient screens used in televisions, computer monitors, and smartphones. The blue LED also paved the way for the invention of the blue laser, which significantly increased storage capacity and led to the creation of new industry standards. Nakamura's invention of the blue LED was a game-changer in the lighting industry, with far-reaching implications for energy-efficient lighting and display technologies. Despite the immense impact of his invention, it is worth noting that Nakamura initially received only 180 for his discovery.

• https://www.nobelprize.org/prizes/physics/2014/press-release/
• https://www.invent.org/blog/inventors/shuji-nakamura-blue-led-lights
• https://en.wikipedia.org/wiki/Shuji_Nakamura

How pigment works ?

Pigments work by selectively absorbing certain wavelengths of light and reflecting or transmitting others, which determines their color. This process is known as subtractive color, and it's what happens with paints and dyes. The absorbed colors are the ones you don't see, and you see only the colors that come bouncing back to your eye. Pigments are insoluble and chemically unreactive substances that are used to add color or change the visual appearance of materials. They can be organic or inorganic, with inorganic pigments being brighter and longer-lasting than organic ones. The color of pigments arises because they absorb only certain wavelengths of visible light, and this absorption is the result of charge transfer between a metal for inorganic pigments, or conjugated double bonds for organic pigments. In summary, pigments work by selectively absorbing and reflecting light to produce the colors we see.

Is there blue pigment in the nature ? And how the blue pigment is different from others ?

Natural blue pigments are rare in nature. The color blue found in foods, plants, and animals often lacks a chemical compound that makes them blue, making natural blue pigment uncommon. For instance, the majority of natural blue colors in food are derived from the purple compound "anthocyanin," which can also produce rich red and black shades in addition to blue and purple. Blue colors in animals are often achieved through the manipulation of light, such as the layered nanostructures on the wing scales of blue-winged butterflies or the light-scattering structures in blue birds. The only known animal to produce a true blue pigment is the Obrina Olivewing butterfly, which uses a bile pigment called biliverdin to achieve its blue color. Therefore, true blue pigments are rare in nature, and the blue color in living organisms is often achieved through structural or optical means rather than through a chemical pigment.

Feb 10

Tweet

Elad Hazan: most exciting paper ever from our @GoogleAI lab at @Princeton : @naman33k @danielsuo @XinyiChen2 https://arxiv.org/abs/2312.06837 Convolutional filters predetermined by the theory, no learning needed!

Yann Lecun: Hi Elad, basically, you are backpedaling to the old days of statistical pattern recognition when the entire literature was devoted to designing fixed feature extractors by exploiting the properties of the signal.
Pretty much every signal classification system (including speech recognition) used to use spectral features of some kind (computed by banks of wavelet-like filters), followed by some sort of non-linearity (e.g. square complex modulus) and temporal pooling.
9 years ago, Stéphane Mallat made a bet with me: he claimed that he could use features from his scattering transform followed by a simple trainable classifier and approach the performance of convnets on ImageNet. The scattering transform architecture looks very much like a convnet (and is inspired by it), except that the kernels are fixed wavelets, whereas in convnets, the kernels are all learned. The prize was a meal at a 3-star restaurant. He admitted defeat after 5 years. He said, "I learned something while trying to understand why it didn't work."

Elad Hazan: Hi Yann, is this is an invitation for another bet? I'll take it, worst case I have to pay for dinner w. you ;-) These particular filters have a nice property: they can learn a symmetric LDS without dependence on the condition number (effective memory). The hope is they can help, but not replace, the learned component. My bet would be that this improves performance, in tasks requiring long memory. But there is no claim that it replaces learning completely.
You would agree that in certain cases, pre-fixed structures can help? That is the case, for example, in adaptive gradient methods for optimization, where a predefined rule is helpful.

• https://twitter.com/ylecun/status/1756249264500380129
• https://arxiv.org/abs/2312.06837

Feb 9

MiniCPM

• 2B/3B
• Deployment on Mobile Phone, ~10 tokens/s
• 面壁智能
• https://github.com/OpenBMB/MiniCPM

compare fine-tuned version v.s. sft version

projects

copilot

Goals

Run local copilot with llm.

References

• https://developer.apple.com/metal/pytorch/
• https://github.com/vllm-project/vllm

library

https://arxiv.org/pdf/2407.01219

Vector Database

• https://github.com/milvus-io/milvus
• https://github.com/qdrant/qdrant

papers

Pollux

Pollux: Co-adaptive Cluster Scheduling for Goodput-Optimized Deep Learning

Aurick Qiao et al. Petuum, Inc, CMU, UCB, MBZUAI.

Award: Jay Lepreau Best Paper

集群调度的目标是为任务分配资源使得

• 训练时间小
• 集群资源利用率高
• 同时保证公平

一般的调度场景是提交任务时声明所需要的资源， 但是在深度学习任务中，资源分配之后如何高效地利用这些资源取决于如 batch size 和 learning rate 之类的配置， 这些配置一般是不属于调度器管理的。

由此，推出

Pollux: Co-adaptive Cluster Scheduler for DL

Pollux, 一个深度学习模型训练的互适应集群调度器。

它自动且动态地，

• 从整个集群性能和公平性的角度进行资源分配
• 动态地为每个训练任务调整 batch size 和 learning rate

什么是互适用呢？

• 从整个集群的角度做调度决策，同时在任务级别来调整配置参数
• 单个任务会根据分配到的资源进行动态地调整

结果：减少用户对任务的手动配置，缩短平均训练时常达 37-50%.

背景

分布式-深度学习模型训练-数据并行

• mini-batch 组装 [batch size]
• 计算梯度 [计算密集]
• 梯度聚合 [通信密集]
• 更新参数 [learning rate]

batch size 对系统吞吐量的影响

总体而言，大 batch size 能够获得更高的吞吐，同时，这个关系

• 次/非线性 sub-linearly
• 一定量后，系统吞吐量不再随 gpu 数量提高

问题在与多少 gpu 合适呢？这取决于

• 模型结构
• 集群硬件
• batch size 等等

batch size 对实际效率（statistical efficiency）的影响

batch size 并不是越大越好，因为

• 需要同时调整 lr 确保模型质量（quality）[应该包括精度等指标]，而且很困难
• 增大 batch size 会降低训练的实际效率
• 大 batch size 训练的模型往往泛化能力较差

对于同样的数据吞吐量，大 batch size 所需要的训练步数较少（二者乘积一定）， 但从模型实际效率上看，大 batch size 需要更多步数才能取得同样的效果。

整体训练性能

经过以上分析可以看出，系统的吞吐和实际的训练效率相对于 batch size 存在经典的 trade-off 曲线关系， 所以，我们需要也可以找出他们的交叉点以获得优化的整体训练性能。

然而，这个优化的 batch size 在不同的训练阶段是不一样的（即随时间变化），最多能差 10 倍多 [McCandlish et al. 2018]， 即需要随训练进度动态地调整 batch size 以获得最优。

集群调度

• 所需的 GPU 的取决于 batch size，反过来，合适的 batch size 取决于分配到的 GPU 资源
• 集群的 GPU 分配是一个全局的决策，需要考虑公平和竞争
• batch size 的选取还取决于具体的任务，如拓展性和实际效率等

用户在提交任务时决定上述配置是困难的，所以由集群调度器来调整这些配置参数可以优化深度学习模型的训练。

Pollux 集群调度器

• 计算模型训练性能
• 计算最优的资源需求、batch size、learning rate
• 重新调整资源分配，相应地调整任务的 batch size 和 learning rate

Key Idea: Goodput, not Throughput

Pollux 的优化目标是新定义的深度学习训练性能指标 goodput

$$GOODPU{T}_t(a,m,s)=THROUGHTPUT(a,m,s)\times EFFICIENC{Y}_t(M)$$

其中，

• $a$ 资源分配向量，如 a_n = #GPU, 表示节点 n 上分配的 gpu 数
• $m$ 每个 gpu 的 batch size
• $s$ 梯度聚合步数
• $M$ 总 batch size，$M=|a|\times m\times s$

$(a,m,s)$ 会在训练过程中由 Pollux 自动计算, $THROUGHTPUT(a,m,s)$ 表示的是系统吞吐量（examples/second）， $EFFICIENC{Y}_t(M)$ 表示实际训练效率（progress/example）.

系统吞吐

$$T_{iter}(a,m,s)=s\times T_{grad}(a,m)+(T_{g}rad(a,m{)}^{\gamma }+T_{sync}(a{)}^{\gamma }{)}^{1/\gamma }$$

其中，

• $T_{iter}$ 每步训练的时间
• $T_{grad}$ 计算梯度的时间
• $T_{sync}$ 网络通信的时间
• $s$ 梯度聚合的步数
• $\gamma$ 计算和通信重合度

Pollux 自动，

• 确定合适的 gpu 数量和 batch size
• 使用梯度聚合提高 batch size 达到 gpu 显存上限
• 把任务尽可能放置（pack）在尽量少的节点上以减少网络负载

实际训练效率

每一个任务的实际训练效率可以表示为 $$EFFICIENC{Y}_t(M)=\frac{{\varphi }_t+M_0}{{\varphi }_t+M}$$

其中，

• $M_0$ 表示用户定义的 baseline batch size
• $M$ 表示 batch size
• ${\varphi }_t$ 梯度噪声 [McCandlish et al. 2018]

用户可以选择较小的初始 batch size $M_0$，Pollux 会选择不同的 bs 去平衡系统吞吐和实际训练效率。

关于梯度噪声 Gradient noise scale

• 较大的梯度噪声 -> 使用较大的 mini-batch 能够获得较高的实际效率
• 接近收敛的低信噪比 -> 更好的实际训练效率

Pollux 能够在不进行提前训练的情况下使用 $(a,m,s)$ 计算出任务的 GOODPUT.

任务优化

在特定分配 gpu 为 a 的前提下，计算最优

$$m^{\ast },s^{\ast }={\mathrm{argmax}}_{m,s}GOODPU{T}_t(a,m,s)$$

改变 batch size 的同时，learning rate 也需要同步改变。 Pollux 为用户提供更新策略

• Linear scaling
• Square-root scaling
• AdaScale (Johnson et al. 2020)

集群优化

优化目标

$$FITNES{S}_p(A)={\left(\frac{1}{J}\sum _{j=1}^{J}SPEEDU{P}_j(A_j{)}^p\right)}^{1/p}$$

其中，

$$SPEEDU{P}_j(A_j)=\frac{ma{x}_{m,s}GOODPU{T}_j(A_j,m_j,s)}{ma{x}_{m,s}GOODPU{T}_j(a_j,m_j,s)}$$

p 是可变参数，用于控制任务间的公平性。

找到分配矩阵 $A$, $A_{j}n$ 表示节点 n 上分配给任务 j 的 gpu 数量。

• 对 A 的寻找使用 metaheuristic algorithm
• 调度要避免频繁的重新分配
• 避免分布式任务共享节点

Pollux 效果评估

Pollux 带来的主要收益是在共享集群上自动配置任务。

重点评估目标：即使任务已经给定理想的静态配置，Pollux 仍然能够相比于传统集群调度器有所提升。 包括以下方面，

• 真实的 Microsoft 深度学习分布式训练集群 (Jeon et al. 2019).
• 不同场景训练任务混合：图像分类、目标检测、语音识别、问答、推荐
• 手动配置 gpu 数量、batch size、learning rate、梯度聚合参数 (不使用Pollux，设定强baseline)

实验数据表明 Pollux 能比专家配置的任务缩短 37-50% 的平均训练时间。

总结

• Pollux 同时从集群和任务的角度对任务参数进行优化
• Pollux 引入 goodput 概念，一种结合系统吞吐和实际效率的衡量标准
• Pollux 实测缩短 37-50% 的平均训练时间

Scaling distributed training with adaptive summation

Saeed Maleki et al. Microsoft Research

Key point

$$\vec{g}=\left(1-\frac{{\vec{g}}_1\cdot {\vec{g}}_2}{2|{\vec{g}}_1{|}^2}\right){\vec{g}}_1+\left(1-\frac{{\vec{g}}_1\cdot {\vec{g}}_2}{2|{\vec{g}}_2{|}^2}\right){\vec{g}}_2$$

Reference

• AdaSum with Horovod
• Arxiv

Adaptation 调研

分布式训练的参数调整基本有以下共识：

• learning rate 需要根据 batch size 调整，以保证训练速度和收敛速度
• learning rate 需要在不同的训练阶段进行调整以保证收敛速度 这个方向早期已经有比较多的研究，包括
• Adam 2015，在 variance 较大时使用较小的 lr
• LARS 2017，在每一层上根据 wieghts 和 gradient 适配 lr
• LARM 2019，LARS + Adam

Summary

最新的一些参数适配相关的研究进展：

An Empirical Model of Large-Batch Training Sam McCandlish et al. OpenAI 总结：文章提出了 gradient noise scale (GNS) 指标，

• GNS 较小，保持 batch size
• GNS 较大，增大 batch size GNS 定义如下：

$$B_{noise}=\frac{tr(H\Sigma )}{G^{T}HG}$$

其中：

• $G$ true gradient
• $H$ true Hessian at parameter values
• $\Sigma$ covariance matrix

AdaScale SGD: A Scale-Invariant Algorithm for Distributed Training

Tyler B. Johnson et al. APPLE

Key: Gradient variance

总结：adascale 根据 gradient variance 来调整 learning rate，提供稳定算法保证在不同的 batch size 下都能找到合适的 lr 保证快速收敛。

Scaling distributed training with adaptive summation

Saeed Maleki et al. Microsoft Research

总结：adasum 利用 gradients 自身的数学性质，提出了一种新的 combine 方法，使得 merge 结果受 outlier 影响较小，更加“合理”，从而加快收敛。 算法公式

$$\vec{g}=\left(1-\frac{{\vec{g}}_1\cdot {\vec{g}}_2}{2|{\vec{g}}_1{|}^2}\right){\vec{g}}_1+\left(1-\frac{{\vec{g}}_1\cdot {\vec{g}}_2}{2|{\vec{g}}_2{|}^2}\right){\vec{g}}_2$$

示意图

Large-Scale Distributed Second-Order Optimization Using Kronecker-Factored Approximate Curvature for Deep Convolutional Neural Networks

Kazuki Osawa et al. Tokyo Institute of Technology, NVIDIA

Key: Gradient second-order metrics

总结：通过计算 Fisher 矩阵二阶信息来更新梯度，每一步的计算速度会慢于sgd，但收敛所需步数减少，特别是在large scale 的场景，整体效率接近。

Resource Elasticity in Distributed Deep Learning

Andrew Or et al. Princeton University, Google AI

autoscaling engine/system

总结：文章对分布式场景深度学习场景下的弹性系统进行了比较全面的梳理，着墨较多在弹性触发条件，但是对系统实现没有什么深入描述，这块创新性存疑。文中提到了依赖 tensorflow 和 horovod，姑且认为，默认实现方案为 horovod 吧。 文中涉及到了类似慢节点检测的机制（STRAGGLER DETECTION）。

KungFu: Making Training in Distributed Machine Learning Adaptive

Luo Mai et al. Imperial College London

总结：kungfu 提供了一个统一的框架/库用来在分布式训练场景下进行不同 adaptation 操作，包括 api、monitor、control 等部分。

• 提供 adaption policies 来定义不同 adaption
• 内置 monitoring，使得依赖各种 metric 做 adapt 决策变得容易
• 分布式参数 adpating 机制：提供弹性、异步 collective

Pollux: Co-adaptive Cluster Scheduling for Goodput-Optimized Deep Learning

Aurick Qiao et al. Petuum, Inc, CMU, UCB, MBZUAI.

总结：提出 goodput 指标用于计算优化配置，包括资源和参数 lr、bs。兼容其他如 adascale 策略。

Efficient Strong Scaling Through Burst Parallel Training

Seo Jin Park et al.

DeepPool key ideas: 引入 foreground/background jobs

• burst parallelism
• GPU multiplixing ： gpu 共享

batch-optimal scaling : find (throughtput , sample efficienncy) for best time to accuracy

• Optimizing time-to-accuracy requires small per-GPU batches at large scale
• Strong scaling and small per-GPU batches are more effective with fast networks
• None of the approaches achieve perfect linear scaling

实现特点：

• 控制粒度到 layer 级别
• 针对静态图，生成 parallel training plan
• stage 级别分配 gpu
• backgroud job 为单机低优任务，现实场景中可能比较难有这样类型任务

总结：依靠并行调度和 gpu 共享技术，引入foreground/background 区分任务优先级，计算最优资源需求和参数配置以获得高资源利用率和任务完成时效。 有点偏向于自动并行和 gpu 共享的混合产物。

strong scaling : hold global batch size constant, decrease per-GPU batch size weak scaling: increase the global batch size correspondingly, per-GPU batch size kept constant

Gradient Descent

An overview of gradient descent optimization algorithms

Sebastian Ruder, Insight Centre for Data Analytics, NUI Galway Aylien Ltd., Dublin, 2017

Gradient Descent

Multi-variable function $f:{\mathbb{R}}^n\mapsto \mathbb{R}$, defined differentiable in a neighborhood of a point $x\in {\mathbb{R}}^n$, for $\lambda \in {\mathbb{R}}_{+}$ small enough,

$$x_{k+1}=x_k-\lambda \nabla f(x_k)$$

leads to $f(x_{k+1})\le f(x_k)$.

If $f$ convex and $\nabla f$ Lipschitz, $f_k$ converge to a local mimimum.

Optimization : Momentum

Let $\gamma <1,v_0=0$,

$$\{\begin{array}{ll}{x}_{k+1}& =x_k-v_k\\ v_k& =\gamma v_{k-1}+\lambda \nabla f(x_k)\end{array}$$

The momentum term increases for dimensions whose gradients point in the same directions and reduces updates for dimensions whose gradients change directions. As a result, we gain faster convergence and reduced oscillation.

Optimization : Nesterov

Version with correction,

$$\{\begin{array}{ll}{x}_{k+1}& =x_k-v_k\\ v_k& =\gamma v_{k-1}+\lambda \nabla f(x_k-\gamma v_{k-1})\end{array}$$

This anticipatory update prevents us from going too fast and results in increased responsiveness, which has significantly increased the performance of RNNs on a number of tasks.

Adagrad

It adapts the learning rate to the parameters, performing larger updates for infrequent and smaller updates for frequent parameters. For this reason, it is well-suited for dealing with sparse data.

$$x_{k+1}=x_k-\frac{\lambda }{\sqrt{G_k+ϵ}}\nabla f(x_k)$$

Application : learned to recognize cats in Youtube videos; GloVe word embeddings.

Adadelta

Adadelta is an extension of Adagrad that seeks to reduce its aggressive, monotonically decreasing learning rate. Instead of accumulating all past squared gradients, Adadelta restricts the window of accumulated past gradients to some fixed size w.

RMSprop

RMSprop and Adadelta have both been developed independently around the same time stemming from the need to resolve Adagrad’s radically diminishing learning rates. RMSprop in fact is identical to the first update vector of Adadelta.

$$v(x_k,t):=\gamma v(x_k,t-1)+(1-\gamma )(\nabla f(x_k){)}^2$$

$$x_{k+1}=x_k-\frac{\eta }{\sqrt{v(x_k,t)}}\nabla f(x_k)$$

Adam

Adaptive Moment Estimation (Adam) is another method that computes adaptive learning rates for each parameter. In addition to storing an exponentially decaying average of past squared gradients $v_k$ like Adadelta and RMSprop, Adam also keeps an exponentially decaying average of past gradients $u_k$, similar to momentum.

$$\begin{array}{rl}{u}_k& ={\beta }_1{u}_{k-1}+(1-{\beta }_1)\nabla f(x_k)\\ v_k& ={\beta }_2{v}_{k-1}+(1-{\beta }_2){\nabla }^{2}f(x_k)\end{array}$$

$$\begin{array}{rl}{\hat{u}}_k& =\frac{u_k}{1-{\beta }_1}\\ {\hat{v}}_k& =\frac{v_k}{1-{\beta }_2}\end{array}$$

$$x_{k+1}=x_k-\frac{\eta }{\sqrt{{\hat{v}}_k}+ϵ}\nabla f(x_k)$$

PyTorch Implementation

SGD, Momentum

The implementation of SGD with Momentum/Nesterov subtly differs from Sutskever et. al. and implementations in some other frameworks.

Considering the specific case of Momentum, the update can be written as

$$\begin{array}{rl}{v}_{t+1}& =\mu \ast v_t+g_{t+1},\\ p_{t+1}& =p_t-\text{lr}\ast v_{t+1},\end{array}$$

where $p$, $g$, $v$ and $\mu$ denote the parameters, gradient, velocity, and momentum respectively.

This is in contrast to Sutskever et. al. and other frameworks which employ an update of the form

$$\begin{array}{rl}{v}_{t+1}& =\mu \ast v_t+\text{lr}\ast g_{t+1},\\ p_{t+1}& =p_t-v_{t+1}.\end{array}$$

The Nesterov version is analogously modified.

# torch/optim/sgd.py

def _single_tensor_sgd(params: List[Tensor],
                       d_p_list: List[Tensor],
                       momentum_buffer_list: List[Optional[Tensor]],
                       *,
                       weight_decay: float,
                       momentum: float,
                       lr: float,
                       dampening: float,
                       nesterov: bool,
                       maximize: bool,
                       has_sparse_grad: bool):

    for i, param in enumerate(params):

        d_p = d_p_list[i]
        if weight_decay != 0:
            d_p = d_p.add(param, alpha=weight_decay)

        if momentum != 0:
            buf = momentum_buffer_list[i]

            if buf is None:
                buf = torch.clone(d_p).detach()
                momentum_buffer_list[i] = buf
            else:
                buf.mul_(momentum).add_(d_p, alpha=1 - dampening)

            if nesterov:
                d_p = d_p.add(buf, alpha=momentum)
            else:
                d_p = buf

        alpha = lr if maximize else -lr
        param.add_(d_p, alpha=alpha)

Adagrad

Algorithm

$$\begin{array}{rl}& \\ & \text{input}:\gamma \text{ (lr)},{\theta }_0\text{ (params)},f(\theta )\text{ (objective)},\lambda \text{ (weight decay)},\\ & \tau \text{ (initial accumulator value)},\eta \text{ (lr decay)}\\ & \text{initialize}:statesu{m}_0\leftarrow 0\\ & \\ & \text{for}t=1\text{to}\dots \text{do}\\ & g_t\leftarrow {\nabla }_{\theta }{f}_t({\theta }_{t-1})\\ & \gamma \leftarrow \gamma /(1+(t-1)\eta )\\ & \text{if}\lambda \ne 0\\ & g_t\leftarrow g_t+\lambda {\theta }_{t-1}\\ & statesu{m}_t\leftarrow statesu{m}_{t-1}+g_t^2\\ & {\theta }_t\leftarrow {\theta }_{t-1}-\gamma \frac{g_t}{\sqrt{statesu{m}_t}+ϵ}\\ & \\ & \mathbf{r}\mathbf{e}\mathbf{t}\mathbf{u}\mathbf{r}\mathbf{n}{\theta }_{\mathbf{t}}\\ & \end{array}$$

# torch/optim/adagrad.py

def _single_tensor_adagrad(params: List[Tensor],
                           grads: List[Tensor],
                           state_sums: List[Tensor],
                           state_steps: List[Tensor],
                           *,
                           lr: float,
                           weight_decay: float,
                           lr_decay: float,
                           eps: float,
                           has_sparse_grad: bool):

    for (param, grad, state_sum, step_t) in zip(params, grads, state_sums, state_steps):
        # update step
        step_t += 1
        step = step_t.item()

        if weight_decay != 0:
            if grad.is_sparse:
                raise RuntimeError("weight_decay option is not compatible with sparse gradients")
            grad = grad.add(param, alpha=weight_decay)

        clr = lr / (1 + (step - 1) * lr_decay)

        if grad.is_sparse:
            grad = grad.coalesce()  # the update is non-linear so indices must be unique
            grad_indices = grad._indices()
            grad_values = grad._values()
            size = grad.size()

            state_sum.add_(_make_sparse(grad, grad_indices, grad_values.pow(2)))
            std = state_sum.sparse_mask(grad)
            std_values = std._values().sqrt_().add_(eps)
            param.add_(_make_sparse(grad, grad_indices, grad_values / std_values), alpha=-clr)
        else:
            is_complex = torch.is_complex(param)
            if is_complex:
                grad = torch.view_as_real(grad)
                state_sum = torch.view_as_real(state_sum)
                param = torch.view_as_real(param)
            state_sum.addcmul_(grad, grad, value=1)
            std = state_sum.sqrt().add_(eps)
            param.addcdiv_(grad, std, value=-clr)
            if is_complex:
                param = torch.view_as_complex(param)
                state_sum = torch.view_as_complex(state_sum)

Adam

Algorithm

$$\begin{array}{rl}& \\ & \text{input}:\gamma \text{ (lr)},{\beta }_1,{\beta }_2\text{ (betas)},{\theta }_0\text{ (params)},f(\theta )\text{ (objective)}\\ & \lambda \text{ (weight decay)},\text{amsgrad},\text{maximize}\\ & \text{initialize}:m_0\leftarrow 0\text{ ( first moment)},v_0\leftarrow 0\text{ (second moment)},{\widehat{v_0}}^{max}\leftarrow 0\\ & \\ & \text{for}t=1\text{to}\dots \text{do}\\ & \text{if}\text{maximize}:\\ & g_t\leftarrow -{\nabla }_{\theta }{f}_t({\theta }_{t-1})\\ & \text{else}\\ & g_t\leftarrow {\nabla }_{\theta }{f}_t({\theta }_{t-1})\\ & \text{if}\lambda \ne 0\\ & g_t\leftarrow g_t+\lambda {\theta }_{t-1}\\ & m_t\leftarrow {\beta }_1{m}_{t-1}+(1-{\beta }_1)g_t\\ & v_t\leftarrow {\beta }_2{v}_{t-1}+(1-{\beta }_2)g_t^2\\ & \widehat{m_t}\leftarrow m_t/(1-{\beta }_1^t)\\ & \widehat{v_t}\leftarrow v_t/(1-{\beta }_2^t)\\ & \text{if}amsgrad\\ & {\widehat{v_t}}^{max}\leftarrow \mathrm{m}\mathrm{a}\mathrm{x}({\widehat{v_t}}^{max},\widehat{v_t})\\ & {\theta }_t\leftarrow {\theta }_{t-1}-\gamma \widehat{m_t}/(\sqrt{{\widehat{v_t}}^{max}}+ϵ)\\ & \text{else}\\ & {\theta }_t\leftarrow {\theta }_{t-1}-\gamma \widehat{m_t}/(\sqrt{\widehat{v_t}}+ϵ)\\ & \\ & \mathbf{r}\mathbf{e}\mathbf{t}\mathbf{u}\mathbf{r}\mathbf{n}{\theta }_{\mathbf{t}}\\ & \end{array}$$

# torch/optim/adam.py

def _single_tensor_adam(params: List[Tensor],
                        grads: List[Tensor],
                        exp_avgs: List[Tensor],
                        exp_avg_sqs: List[Tensor],
                        max_exp_avg_sqs: List[Tensor],
                        state_steps: List[Tensor],
                        *,
                        amsgrad: bool,
                        beta1: float,
                        beta2: float,
                        lr: float,
                        weight_decay: float,
                        eps: float,
                        maximize: bool):

    for i, param in enumerate(params):

        grad = grads[i] if not maximize else -grads[i]
        exp_avg = exp_avgs[i]
        exp_avg_sq = exp_avg_sqs[i]
        step_t = state_steps[i]
        # update step
        step_t += 1
        step = step_t.item()

        bias_correction1 = 1 - beta1 ** step
        bias_correction2 = 1 - beta2 ** step

        if weight_decay != 0:
            grad = grad.add(param, alpha=weight_decay)

        # Decay the first and second moment running average coefficient
        exp_avg.mul_(beta1).add_(grad, alpha=1 - beta1)
        exp_avg_sq.mul_(beta2).addcmul_(grad, grad.conj(), value=1 - beta2)
        if amsgrad:
            # Maintains the maximum of all 2nd moment running avg. till now
            torch.maximum(max_exp_avg_sqs[i], exp_avg_sq, out=max_exp_avg_sqs[i])
            # Use the max. for normalizing running avg. of gradient
            denom = (max_exp_avg_sqs[i].sqrt() / math.sqrt(bias_correction2)).add_(eps)
        else:
            denom = (exp_avg_sq.sqrt() / math.sqrt(bias_correction2)).add_(eps)



        step_size = lr / bias_correction1
        # param = param - step_size * (exp_avg / denom)
        # element-wise division
        param.addcdiv_(exp_avg, denom, value=-step_size)

AdamW is Adam with correct Weight Decay, when weight decay is 0, there is no difference between Adam and AdamW.

Auto Parallel

A Survey on Auto-Parallelism of Neural Networks Training

Peng Liang et. al. National University of Defense Technology

Abstract.

DL --> large model --> distributed training --> heterogeneous cluster --> auto-parallelism --> large scale DL model

• basic parallelism schemes, communication cost and memory consumption
• current works, strategies, common methods
• promising trends

Introduction

Parallelism strategy

• intra-operator parallelism: data parallelism, tensor parallelism (intra-layer model parallelism)
• inter-operator parallelism: inter-layer model parallelism, pipeline parallelism

Hybrid parallelism

• data + model + pipeline
• Megatron-LM, DeepSpeed(3D parallelism)

Manual --> Auto

All practicable works: a few combinations of parallelism schemes, weak scalability

• e.g. cost model
• automatic parallelism search space can be further expanded
• heterogeneous devices, communication pace/topology

Challenges

• detailed analysis of different parallelism schemes
• trade-offs between different parallelism schemes
• load-balance across heterogeneous devices
• optimization of network communication
• trade-off between runtime and strategy performance in finding strategy

Parallelism schemes

Data parallelism

• Vanilla DP
• ZeRO-Powered DP
• Communication of DP, Centralized/Decentralized architecture

ZeRO-DP

• three stages: 1 partition optimizer states, 2 partition gradients and optimizer states, 3 partition parameters
• stage 1 and 2: reduce-scatter accumulated gradients, stage 3: all-gather updated parameters
• solve redundancy problem with 50% more communication volume (all-gather)

Model parallelism

• Intra-layer MP, tensor parallelism, partition weight tensor
• Inter-layer MP

Pipeline parallelism

The partition pattern of PP is the same as that of MP

• inter-layer MP = PP
• PP = well scheduled pipelined MP
• overlap computation, solve low-utility of MP

PipeDream (Microsoft), GPipe (Google)

Strategy Searching Methods for Auto-Parallelism

• NP-hard problem
• classic-algorithm-based v.s. machine-learning-based

Classic-algorithm based methods

• recursive algorithm
• dynamic programming algorithm
• integer linear programming algorithm
• breath-first-search (BFS) algorithm

Machine-learning based methods

• Monte-Carlo Markov Chain (MCMC)
• Monte-Carlo Tree Search (MCTS)
• reinforcement learning

Conclusion

Accelerating strategy searching

• grouping
• profiling-base cost model
• using heuristics

Optimizing parallelism strategies

• topology-aware computation
• topology-aware communication

Supporting more parallelism schemes

Scheduling

OSDI 2021

Pollux: Co-adaptive Cluster Scheduling for Goodput-Optimized Deep Learning

Aurick Qiao et al. Petuum, Inc, CMU, UCB. Cite 14

Award: Jay Lepreau Best Paper

Oort: Efficient Federated Learning via Guided Participant Selection

Fan Lai et al. University of Michigan. Cite 17

PET: Optimizing Tensor Programs with Partially Equivalent Transformations and Automated Corrections

Haojie Wang et al. Tsinghua, CMU, FB. Cite 4

Privacy Budget Scheduling

Tao Luo et al. Columbia University, Microsoft Research. Cite 2

OSDI 2020

Providing SLOs for Resource-Harvesting VMs in Cloud Platforms

Pradeep Ambati et al. Microsoft Azure, Microsoft Research.

The CacheLib Caching Engine: Design and Experiences at Scale

Benjamin Berg et al. CMU, FB, MS

Twine: A Unified Cluster Management System for Shared Infrastructure

Chunqiang Tang et al. FB

FIRM: An Intelligent Fine-Grained Resource Management Framework for SLO-Oriented Microservices

Haoran Qiu et al. UIUC

Building Scalable and Flexible Cluster Managers Using Declarative Programming

Lalith Suresh et al. VMware, IST, UIUC ...

Protean: VM Allocation Service at Scale

Ori Hadary et al. MS

1-Bit Stochastic Gradient Descent and its Application to Data-Parallel Distributed Training of Speech DNNs

Frank Seide et. al. MRA, Tsinghua, MR. INTERSPEECH 2014

1-BitSGDwithErrorFeedback

总结：1 bit 量化的思路是把每个 32 位的值按照正负分别用 1 或 0 量化，然后通信，更新时 1 则更新 +1，0 则更新 -1，同时本地量化差异保留补偿。这样的效果应该相当于本地梯度累积，每次通信根据正负同步 1 个单位的量。当然，每次更新带上系数，最终得到收敛效果。

1-bit Adam: Communication Efficient Large-Scale Training with Adam’s Convergence Speed

Hanlin Tang et. al. Microsoft 2021

总结：adam 的梯度更新是非线性的，沿用 1bit sgd 的梯度累积或者说梯度补偿策略无法保证收敛性，文章提出了针对 adam 的保证收敛性的 1 bit 梯度量化方法并给出了理论证明。方法的核心在于在 warm-up 阶段计算 momentum 方差，用于后面进行误差补偿的计算。

Deep Gradient Compression: Reducing the Communication Bandwidth for Distributed Training

Yujun Lin et. al. Tsinghua University, ICLR 2018

DGC

Optimizing Network Performance for Distributed DNN Training on GPU Clusters: ImageNet/AlexNet Training in 1.5 Minutes

Peng Sun et. al. SenseTime 2019

GRACE: A Compressed Communication Framework for Distributed Machine Learning

Hang Xu et. al. 2021

GRACE

SIDCo An Efficient Statistical-based Gradient Compression Technique for Distributed Training Systems

Ahmed M. Abdelmoniem et. al. CEMSE, KAUST. MLSys 2021

SIDCo

PowerSGD: Practical Low-Rank Gradient Compression for Distributed Optimization

Thijs Vogels et. al. EPFL. NeurIPS 2019

PowerSGD

Don't Use Large Mini-Batches, Use Local SGD

Tao Lin et. al. EPFL. ICLR 2020

LocalSGD-Code

总结：为减轻同步模式中慢节点的影响，可以减少通信，这会带来精度损失。使用 local SGD 的方法可以现在节点内进行 SGD 更新，多步之后再同步各个节点上的参数。 post local SGD 的方法将训练过程分成两阶段：先使用同步 SGD，再增大同步间隔提高训练吞吐。

Adaptive Communication Strategies to Achieve the Best Error-Runtime Trade-off in Local-Update SGD

Jianyu Wang et. al. Carnegie Mellon University. SysML 2019

总结：动态调整参数同步间隔来平衡训练吞吐和精度。

Overlap Local-SGD: An Algorithmic Approach to Hide Communication Delays in Distributed SGD

Jianyu Wang et. al. ICASSP 2020

Overlap_Local_SGD

Deep Gradient Compression

Deep Gradient Compression: Reducing the Communication Bandwidth for Distributed Training

Yujun Lin et. al. Tsinghua University, ICLR 2018

Gradient exchange require hight bandwidth,

• high latency
• low throughput
• poor connnection

Since 99.9% of the gradient exchange in SGD are redundant, propose deep gradient compression (DGC) to reduce communication bandwidth.

To preserve accuracy,

• momentum correction
• local gradient clipping
• momentum factor masking, alleviate staleness
• warm-up training

improving local gradient accumulation and overcomming the staleness effect

DGC

• pushes the gradient compression ratio to up to 600×
• not need to change the model structure
• no loss of accuracy

How

• only gradients larger than a threshold are transmitted
• accumulate the rest of the gradients locally, local gradient accumulation is equivalent to increasing the batch size over time

DGC naively perform fine-grained (i.e., element-wise) top-k to select gradients, and thus the communication will suffer from increased allgather data volume as #nodes increases.

CSC modified the process with coarse-grained sparsification: gradients are partioned into chunks, allreduce the gradient chunks selected based on allreduced L1-norm of each chunk, which gets rid of the allgather and solves the problem.

deep-gradient-compression github

配置

# configs/dgc/__init__.py

训练流程

# train.py

from dgc.horovod.optimizer import DistributedOptimizer

# from dgc.compression import DGCCompressor
compression = configs.train.compression()
# cpr_parameters 即 dgc 处理的范围
compression.initialize(cpr_parameters.items())

# from dgc.optim import DGCSGD
optimizer = configs.train.optimizer(model.parameters())

# Horovod: wrap optimizer with DistributedOptimizer.
optimizer = DistributedOptimizer(
    optimizer, named_parameters=model.named_parameters(),
    compression=compression,
    backward_passes_per_step=configs.train.num_batches_per_step,
    op=hvd.Average
)

# 训练基本循环 zero_grad -> loss.backward -> optimizer.step
# 特别注意这里多次 backward 才走一次 step 更新
model.train()
for step, (inputs, targets) in enumerate(...):
    optimizer.zero_grad()

    # 这里用了内置循环累积梯度，比直接使用大 batch 剩显存
    # 注意这个 for 循环，对于 optimizer 里面理解 synchronize 过程非常重要
    for b in range(0, step_size, batch_size):
        _inputs = inputs[b:b+batch_size]
        _targets = targets[b:b+batch_size]
        _outputs = model(_inputs)
        _loss = criterion(_outputs, _targets)
        _loss.mul_(_r_num_batches_per_step)
        _loss.backward()
        loss += _loss.item()
    optimizer.step()

Optimizer

# dgc/horovod/optimizer.py

class _DistributedOptimizer(torch.optim.Optimizer):

    def __init__(self, ...):
        # 初始化最后注册 通信 hook
        self._register_hooks()
        
    def _register_hooks(self):
        for param_group in self.param_groups:
            for p in param_group['params']:
                if p.requires_grad:
                    # 注册函数只执行一次，这里 zero grad 不是每次调用 hook
                    p.grad = p.data.new(p.size()).zero_()
                    self._requires_update.add(p)
                    # 创建幽灵 tensor 来累积梯度，节点间同步，直至更新；
                    # p_tmp 和 p 使用同样的 storage，不占用额外显存
                    p_tmp = p.expand_as(p)
                    grad_acc = p_tmp.grad_fn.next_functions[0][0]
                    # 注册 _make_hook 这个关键 hook
                    grad_acc.register_hook(self._make_hook(p))
                    self._grad_accs.append(grad_acc)

    def _make_hook(self, p):
        # 这个 hook 有一个计数器，_allreduce_delay, 根据对象 p 不一样可以取不一样的值
        # 计数器不为零时跳过，这样可以让 grad 在本地累积，因为这个 hook 是做通信的
        # 效果为这个 hook 在多次调用才会被执行一次
        def hook(*ignore):
            handle, ctx = None, None
            self._allreduce_delay[p] -= 1
            if self._allreduce_delay[p] == 0:
                handle, ctx = self._allreduce_grad_async(p)
            self._handles[p] = (handle, ctx)
        return hook

    # 然后主要流程 step
    def step(self, closure=None):
        self.synchronize()
        return super(self.__class__, self).step(closure)

    # step 调用 synchronize, 可以有跳过逻辑
    def synchronize(self):
        # 处理 hook 注册不成功，或者说 hook 没有被调用
        # hook 被调用后会添加 self._handles
        missing_p = self._requires_update - set(self._handles.keys())
        for p in missing_p:
            handle, ctx = self._allreduce_grad_async(p)
            self._handles[p] = (handle, ctx)

        # handle 为 None 的 hook 跳过又不跳过了？
        # 需要注意 synchronize 函数每个 step 被调用，但不是每次 backward 都会被调用
        # 在之前的 train 中有每个 step 会多次 backward，所以 grad 的 hook 会被多次调用，次数匹配
        # 所以代码执行到这里 handle 应该是一次调用_allreduce_grad_async 如果不是就补上
        for p, (handle, ctx) in self._handles.items():
            if handle is None:
                handle, ctx = self._allreduce_grad_async(p)
                self._handles[p] = (handle, ctx)

        # for 循环处理异步通信的结果
        for p, (handle, ctx) in self._handles.items():
            output = self._synchronize_(handle)
            # 重置本地累积次数
            self._allreduce_delay[p] = self.backward_passes_per_step
            # 解压更新梯度
            p.grad.set_(self._compression.decompress(output, ctx))

        # 执行完毕，清理
        self._handles.clear()

    # 异步通信的 op，核心逻辑在 compression 中
    def _allreduce_grad_async(self, p):
        name = self._parameter_names.get(p)
        tensor_compressed, ctx = self._compression.compress(p.grad, name)

        handle = self._communicate_(tensor_compressed, name=name, op=self.op)
        return handle, ctx

• hook 函数是一次注册，多次调用，所以 self._handles 会不断被填充，每次 synchronize 后可以被 clear

# dgc/compression.py

class DGCCompressor:
    def __init__(self, ...):
        self.attributes = {}

    def initialize(self, named_parameters):
        # 工作范围
        for name, param in named_parameters:
            self.attributes[name] = (numel, shape, num_selects, num_samples, top_k_samples, sample_stride)

    def _sparsify(self, tensor, name):
        # 选出稀疏的 tensor 去通信
        # 原实现中比较复杂
        # 先随机选取部分梯度值的 TOPK 来计算阈值
        # 然后通过该阈值对原 tensor 做稀疏化
        importance = tensor.abs()
        mask = torch.ge(importance, threshold)
        indices = mask.nonzero().view(-1)
        num_indices = indices.numel()
        # 这里实现上有个 for 循环确保选出的 topk 满足要求
        indices = indices[:num_selects]
        values = tensor[indices]
        return values, indices, numel, shape, num_selects

    def compress(self, tensor, name):
        if self.compress_ratio < 1.0 and name in self.attributes:
            # compress
            tensor_compensated = self.memory.compensate(tensor, name, accumulate=True)
            values, indices, numel, shape, num_selects = self._sparsify(tensor_compensated, name)
            self.memory.update(name, (indices, ))
            return tensor, ctx
        else:
            return tensor, ctx

    def decompress(self, tensor, ctx):
        name, numel, shape, vdtype, idtype, grad = ctx
        if self.compress_ratio < 1.0 and name in self.attributes:
            # 这里的 tensor 是个 tuple
            # decompress
            values, indices = tensor
            # 把同步回来的稀疏 tensor 对应位置更新
            # accumulate=True 处理 indices 中有重复的情况
            grad.zero_().index_put_([indices], values, accumulate=True)
            if self.op == Average:
                grad.mul_(1. / self.world_size)
            return grad.view(shape)
        else:
            return self.memory.compensate(tensor, name, accumulate=False)

    # optimizer _communicate_
    def communicate(self, tensor_compressed, name, op):
        # 两个分支
        if self.compress_ratio < 1.0 and name in self.attributes:
            # dgc 分支，tensor_compressed 是 tuple，各个节点选的 topk index 不相同
            # 所以使用 allgather 交换，然后各自解压、更新
            return [allgather_async_(t, name=f'{name}.t{e}')
                    for e, t in enumerate(tensor_compressed)]
        else:
            # 普通分支，直接 allreduce 完整 tensor
            return allreduce_async_(tensor_compressed, name=name, op=op)

    # optimizer _synchronize_
    def synchronize(self, handle):
        # from horovod.torch.mpi_ops import synchronize as synchronize_
        if isinstance(handle, (tuple, list)):
            return [synchronize_(h) for h in handle]
        else:
            return synchronize_(handle)

为了保证精度文章中介绍了下面几种补偿策略

• momentum correction
• local gradient clipping
• momentum factor masking, alleviate staleness
• warm-up training

前三种策略在 Memory 实现

# dgc/memory.py

class DGCSGDMemory(Memory):

    def compensate(self, grad, name, accumulate=True):
        if self.gradient_clipping is not None:
            grad = self.gradient_clipping(grad)
        mmt = self.momentums[name]
        if accumulate:
            # Momentum Correction
            vec = self.velocities[name]
            if self.nesterov:
                mmt.add_(grad).mul_(self.momentum)
                vec.add_(mmt).add_(grad)
            else:
                mmt.mul_(self.momentum).add_(grad)
                vec.add_(mmt)
            return vec
        else:
            if self.nesterov:
                mmt.add_(grad).mul_(self.momentum)
                return mmt.add(grad)
            else:
                mmt.mul_(self.momentum).add_(grad)
                return mmt.clone()  # TODO: save this clone

    def update(self, name, ctx):
        indices = ctx[0]
        if self.momentum_masking:
            self.momentums[name].view(-1).index_fill_(0, indices, 0)
        self.velocities[name].view(-1).index_fill_(0, indices, 0)

Optimizing Network Performance for Distributed DNN Training on GPU Clusters: ImageNet/AlexNet Training in 1.5 Minutes

Peng Sun et. al. SenseTime 2019

Communication backend: GradientFlow

• ring-based allreduce
• mixed-precision training
• computation/communication overlap
• lazy allreduce: fusing multiple communication operations
• coarse-grained sparse communication: only transmitting important gradient chunks

and also,

• momentum SGD correction
• warm-up dense training

DGC naively perform fine-grained (i.e., element-wise) top-k to select gradients, and thus the communication will suffer from increased allgather data volume as #nodes increases.

CSC modified the process with coarse-grained sparsification: gradients are partioned into chunks, allreduce the gradient chunks selected based on allreduced L1-norm of each chunk, which gets rid of the allgather and solves the problem.

Flash Attention

Key Idea

• memory hierachy: GPU HBM (40G, 1.5TB/s) -> GPU SRAM (20MB, 19TB/s)
• tiling: split the input into blocks and make several passes over blocks, fit GPU SRAM
• recompute: recompute attention on-chip in the backward pass instead of retrieve from HBM

Conclusion: increase FLOPs, decrease Wall-clock time

Code Details

Modules

flash_attn

// csrc/flash_attn/fmha_api.cpp

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
    m.doc() = "Fused Multi-head Self-attention";
    m.def("fwd", &mha_fwd, "Forward pass");
    m.def("bwd", &mha_bwd, "Backward pass");
    m.def("fwd_block", &mha_fwd_block, "Forward pass (blocksparse)");
    m.def("bwd_block", &mha_bwd_block, "Backward pass (blocksparse)");
}

fused_dense_lib

// csrc/fused_dense_lib/fused_dense.cpp

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
  m.def("linear_bias_forward", &linear_bias_forward, "linear bias forward");
  m.def("linear_bias_backward", &linear_bias_backward, "linear bias backward");
  m.def("linear_bias_wgrad", &linear_bias_wgrad, "linear bias wgrad");
  m.def("linear_bias_residual_backward", &linear_bias_residual_backward, "linear bias residual backward");
  m.def("linear_gelu_forward", &linear_gelu_forward, "linear gelu forward");
  m.def("linear_gelu_linear_backward", &linear_gelu_linear_backward, "linear gelu linear backward");
  m.def("linear_residual_gelu_linear_backward", &linear_residual_gelu_linear_backward, "linear residual gelu linear backward");
}

fused_softmax

// csrc/fused_softmax/fused_softmax.cpp

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
  m.def("scaled_masked_softmax_forward", &multihead_attn::fused_softmax::scaled_masked_softmax::fwd, "Self Multihead Attention scaled, time masked softmax -- Forward.");
  m.def("scaled_masked_softmax_backward", &multihead_attn::fused_softmax::scaled_masked_softmax::bwd, "Self Multihead Attention scaled, time masked softmax -- Backward.");
  m.def("scaled_masked_softmax_get_batch_per_block", &multihead_attn::fused_softmax::scaled_masked_softmax::get_batch_per_block, "Return Batch per block size.");
  m.def("scaled_upper_triang_masked_softmax_forward", &multihead_attn::fused_softmax::scaled_upper_triang_masked_softmax::fwd, "Self Multihead Attention scaled, time masked softmax -- Forward.");
  m.def("scaled_upper_triang_masked_softmax_backward", &multihead_attn::fused_softmax::scaled_upper_triang_masked_softmax::bwd, "Self Multihead Attention scaled, time masked softmax -- Backward.");
}

layer_norm

// csrc/fused_softmax/fused_softmax.cpp

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
  m.doc() = "CUDA DropoutAddLayerNorm";
  m.def("dropout_add_ln_fwd", &dropout_add_ln_fwd, "Run Dropout + Add + LayerNorm forward kernel");
  m.def("dropout_add_ln_bwd", &dropout_add_ln_bwd, "Run Dropout + Add + LayerNorm backward kernel");
}

rotary

// csrc/rotary/rotary.cpp

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
  m.def("apply_rotary", &apply_rotary, "Apply rotary embedding");
}

xentropy

// csrc/xentropy/interface.cpp

PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
    m.def("forward", &softmax_xentropy_forward, "Softmax cross entropy loss with label smoothing forward (CUDA)");
    m.def("backward", &softmax_xentropy_backward, "Softmax cross entropy loss with label smoothing backward (CUDA)");
}

Flash Attention

Python API call chains overview

FlashSelfAttention
    flash_attn_func
        flash_attn_unpadded_qkvpacked_func
            FlashAttnQKVPackedFunc
                _flash_attn_forward  = flash_attn_cuda.fwd
                _flash_attn_backward = flash_attn_cuda.bwd
    
FlashCrossAttention
    flash_attn_unpadded_kvpacked_func
        FlashAttnKVPackedFunc
            _flash_attn_forward
            _flash_attn_backward
    
# ----
    flash_attn_unpadded_qkvpacked_split_func
        FlashAttnQKVPackedSplitFunc
            _flash_attn_forward
            _flash_attn_backward

FlashBlocksparseAttention
    flash_blocksparse_attn_func 
        FlashBlocksparseAttnFun 
            _flash_blocksparse_attn_forward  = flash_attn_cuda.fwd_block
            _flash_blocksparse_attn_backward = flash_attn_cuda.bwd_block

FlashSelfAttention

# flash_attn/flash_attn_interface.py

def flash_attn_unpadded_qkvpacked_func(qkv, cu_seqlens, max_seqlen, dropout_p, softmax_scale=None, causal=False, return_attn_probs=False):
    return FlashAttnQKVPackedFunc.apply(qkv, cu_seqlens, max_seqlen, dropout_p, softmax_scale, causal, return_attn_probs)

class FlashAttnQKVPackedFunc(torch.autograd.Function):
    @staticmethod
    def forward(ctx, qkv, cu_seqlens, max_seqlen, dropout_p, softmax_scale, causal, return_softmax):
        rng_state = torch.cuda.get_rng_state() if dropout_p > 0 else None
        out, softmax_lse, S_dmask = _flash_attn_forward(...)
        return (out, softmax_lse, S_dmask)

    @staticmethod
    def backward(ctx, dout, *args):
        _flash_attn_backward(...)
        return dqkv, None, None, None, None, None, None

def _flash_attn_forward(q, k, v, out, cu_seqlens_q, cu_seqlens_k, max_seqlen_q, max_seqlen_k,
                        dropout_p, softmax_scale, causal, return_softmax, num_splits=0,
                        generator=None):
    softmax_lse, *rest = flash_attn_cuda.fwd(...)
    return out, softmax_lse, S_dmask


def _flash_attn_backward(dout, q, k, v, out, softmax_lse, dq, dk, dv, cu_seqlens_q, cu_seqlens_k,
                         max_seqlen_q, max_seqlen_k, dropout_p, softmax_scale, causal, num_splits=0,
                         generator=None):
    _, _, _, softmax_d = flash_attn_cuda.bwd(...)
    return dq, dk, dv, softmax_d

mha_fwd + mha_bwd

// csrc/flash_attn/fmha_api.cpp

std::vector<at::Tensor>
mha_fwd(const at::Tensor &q,         // total_q x num_heads x head_size, total_q := \sum_{i=0}^{b} s_i
        const at::Tensor &k,         // total_k x num_heads x head_size, total_k := \sum_{i=0}^{b} s_i
        const at::Tensor &v,         // total_k x num_heads x head_size, total_k := \sum_{i=0}^{b} s_i
        at::Tensor &out,             // total_q x num_heads x head_size, total_k := \sum_{i=0}^{b} s_i
        const at::Tensor &cu_seqlens_q,  // b+1
        const at::Tensor &cu_seqlens_k,  // b+1
        const int max_seqlen_q_,
        const int max_seqlen_k_,
        const float p_dropout,
        const float softmax_scale,
        const bool zero_tensors,
        const bool is_causal,
        const bool return_softmax,
        const int num_splits,
        c10::optional<at::Generator> gen_) {

    Launch_params<FMHA_fprop_params> launch_params(dprops, stream, is_dropout, return_softmax);
    run_fmha_fwd(launch_params);
    return result;
}

void run_fmha_fwd(Launch_params<FMHA_fprop_params> &launch_params) {
    if (launch_params.params.d <= 32) {
        run_fmha_fwd_hdim32(launch_params);
    } else if (launch_params.params.d <= 64) {
        run_fmha_fwd_hdim64(launch_params);
    } else if (launch_params.params.d <= 128) {
        run_fmha_fwd_hdim128(launch_params);
    }
}

std::vector<at::Tensor>
mha_bwd(const at::Tensor &dout,  // total_q x num_heads, x head_size
        const at::Tensor &q,   // total_q x num_heads x head_size, total_q := \sum_{i=0}^{b} s_i
        const at::Tensor &k,   // total_k x num_heads x head_size, total_k := \sum_{i=0}^{b} s_i
        const at::Tensor &v,   // total_k x num_heads x head_size, total_k := \sum_{i=0}^{b} s_i
        const at::Tensor &out,   // total_q x num_heads x head_size
        const at::Tensor &softmax_lse_,     // b x h x s softmax logsumexp
        at::Tensor &dq,   // total_q x num_heads x head_size, total_q := \sum_{i=0}^{b} s_i
        at::Tensor &dk,   // total_k x num_heads x head_size, total_k := \sum_{i=0}^{b} s_i
        at::Tensor &dv,   // total_k x num_heads x head_size, total_k := \sum_{i=0}^{b} s_i
        const at::Tensor &cu_seqlens_q,  // b+1
        const at::Tensor &cu_seqlens_k,  // b+1
        const int max_seqlen_q_,
        const int max_seqlen_k_,          // max sequence length to choose the kernel
        const float p_dropout,         // probability to drop
        const float softmax_scale,
        const bool zero_tensors,
        const bool is_causal,
        const int num_splits,
        c10::optional<at::Generator> gen_
) {
    auto launch = &run_fmha_bwd;

    FMHA_dgrad_params params;

    set_params_dgrad(params, ...  num_splits);

    launch(params, stream, /*configure=*/true);

    launch(params, stream, /*configure=*/false);


    return { dq, dk, dv, softmax_d };
}

void run_fmha_bwd(FMHA_dgrad_params &params, cudaStream_t stream, const bool configure) {
  if (params.d <= 32) {
      run_fmha_bwd_hdim32(params, stream, configure);
  } else if (params.d <= 64) {
      run_fmha_bwd_hdim64(params, stream, configure);
  } else if (params.d <= 128) {
      run_fmha_bwd_hdim128(params, stream, configure);
  }
}

// csrc/flash_attn/src/fmha_fwd_hdim32.cu 

#include "fmha_fwd_launch_template.h"

void run_fmha_fwd_hdim32(Launch_params<FMHA_fprop_params> &launch_params) {
    FP16_SWITCH(launch_params.params.is_bf16, ({
        if (launch_params.params.seqlen_k == 128) {
            using Kernel_traits = FMHA_kernel_traits<128, 32, 16, 1, 4, 0x08u, elem_type>;
            run_fmha_fwd_loop<Kernel_traits>(launch_params);
        } else if (launch_params.params.seqlen_k >= 256) {
            using Kernel_traits = FMHA_kernel_traits<256, 32, 16, 1, 4, 0x08u, elem_type>;
            run_fmha_fwd_loop<Kernel_traits>(launch_params);
        }
    }));
}

run_fmha_fwd_loop -> fmha_fwd_loop_kernel

// csrc/flash_attn/src/fmha_fwd_launch_template.h

template<typename Kernel_traits>
void run_fmha_fwd_loop(Launch_params<FMHA_fprop_params> &launch_params) {
    BOOL_SWITCH(launch_params.is_dropout, IsDropoutConst, ({
        auto kernel = launch_params.params.is_causal
            ? (launch_params.return_softmax
               ? &fmha_fwd_loop_kernel<Kernel_traits, IsDropoutConst, true, true>
               : &fmha_fwd_loop_kernel<Kernel_traits, IsDropoutConst, true, false>)
            : (launch_params.return_softmax
               ? &fmha_fwd_loop_kernel<Kernel_traits, IsDropoutConst, false, true>
               : &fmha_fwd_loop_kernel<Kernel_traits, IsDropoutConst, false, false>);
        if( smem_size >= 48 * 1024 ) {
            FMHA_CHECK_CUDA(cudaFuncSetAttribute(
                kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, smem_size));
        }
        kernel<<<grid, Kernel_traits::THREADS, smem_size, launch_params.stream>>>(
            launch_params.params);
    }));
}

template<typename Kernel_traits, bool Is_dropout, bool Is_causal, bool Return_softmax>
__global__ void fmha_fwd_loop_kernel(FMHA_fprop_params params) {
    fmha::device_1xN_loop<Kernel_traits, Is_dropout, Is_causal, Return_softmax>(params);
}

device_1xN_loop

// csrc/flash_attn/src/fmha_fprop_kernel_1xN.h

#include "fmha_kernel.h"
#include <fmha/kernel_traits.h>
#include <fmha/gemm.h>
#include <fmha/utils.h>

template<typename Kernel_traits, bool Is_dropout, bool Is_causal, bool Return_softmax, typename Params>
inline __device__ void device_1xN_loop(const Params &params) {

    // The block index for the batch.
    const int bidb = blockIdx.x;
    // The block index for the head.
    const int bidh = blockIdx.y;
    // The thread index.
    const int tidx = threadIdx.x;

    fmha::device_1xN_<Kernel_traits, Is_dropout, Is_causal, Return_softmax, true, true>(params, bidb, bidh, STEPS, ph, 0);
}

template<typename Kernel_traits, bool Is_dropout, bool Is_causal, bool Return_softmax, bool Is_first, bool Is_last, typename Params, typename Prng>
inline __device__ void device_1xN_(const Params &params, const int bidb, const int bidh, int steps, Prng &ph, const int loop_step_idx) {
    ...
}

Dependencies

csrc/flash_attn/src/fmha
|-- gemm.h
|-- gmem_tile.h
|-- kernel_traits.h
|-- mask.h
|-- smem_tile.h
|-- softmax.h
`-- utils.h

Long sequences

WHY 原 flashattention 算法的并行依赖 bs * num_heads，A100 有 108 SMs，当 bs * num_heads > 80 时并行度利用率较高，但在长序列场景下，bs * num_heads 通常较小，无法充分利用 GPU 并行度。

HOW 前向：使用多个 thread blocks 并行处理同一个 attention head，head 按照 row 切分，可以无依赖并行 反向：多个 thread blocks 并行处理，head 按照 column 切分，thread 间需要聚合 query gradient。（如果按 row 切则需要聚合 key 和 value 的 gradient）

References

• github
• arxiv
• efficient
• long sequences
• triton

LoRA

arxiv

LoRA: Low-Rank Adaptation of Large Language Models

References

• microsoft LoRA
• huggingface PEFT

self.lora_A = nn.Parameter(self.weight.new_zeros((r, num_embeddings)))
self.lora_B = nn.Parameter(self.weight.new_zeros((embedding_dim, r)))
self.scaling = self.lora_alpha / self.r
# Freezing the pre-trained weight matrix
self.weight.requires_grad = False

result = F.linear(x, T(self.weight), bias=self.bias)
result += (self.lora_dropout(x) @ self.lora_A.T @ self.lora_B.T) * self.scaling

mathematics

Topics

• Pure Mathematics

  • Algebra
  • Analysis
  • Geometry
  • Topology

• Applied Mathematics

  • Modeling
  • Learning
  • Algorithms

• Physics

  • Classical mechanics
  • Statistical mechanics
  • Electromagnetism
  • Quantum mechanics

Definition

Scalar $\alpha \in \mathbf{R}$,

Vector $x\in {\mathbf{R}}^n,n\in \mathbf{N}$

Linearity $f:{\mathbf{R}}^n\mapsto {\mathbf{R}}^m$,

$$f(\alpha \cdot x+\beta \cdot y)=\alpha \cdot f(x)+\beta \cdot f(y)$$

Estimator

Bias

$$\mathrm{bias}(\hat{y})=E(\hat{y})-y$$

$\hat{y}$ unbiased if $\mathrm{bias}(\hat{y})=0$

Mean Squared Error

$$MSE=E((\hat{y}-y{)}^2)=\mathrm{bias}(\hat{y}{)}^2+\mathrm{Var}(\hat{y})$$

minimizing the mean squared error = maximum likelihood estimator

Moore-Penrose pseudoinverse

$$X^{+}=\underset{\alpha \to 0}{\lim }(X^{T}X+\alpha I{)}^{-1}{X}^T$$

Probability

$P(y,x)$ represent the probability of $x$ and $y$ hanppend.

$P(y|x)$ represent the probability of $y$ while $x$ already known been happend.

Conditional Probability

$$P(y|x)=\frac{P(y,x)}{P(x)}$$

If $x$ independent with $y$, then

$$P(y|x)=P(y),P(y|x)=P(y,x)$$

Bayes's Rule

$$P(x|y)=\frac{P(x)P(y|x)}{P(y)}$$

Experience

$$E(x_i)=\sum _{i}P(x_i)x_i$$

$$E(x)={\int }_{I}P(x)xdx$$

Variance

$$\mathrm{Var}(x)=E((x-E(x){)}^2)$$

Covariance

$$\mathrm{Cov}(x,y)=E((x-E(x))(y-E(y)))$$

Normal Distribution

$$\mathcal{N}(\mu ,{\sigma }^2)=\sqrt{\frac{1}{2\pi {\sigma }^2}}{e}^{-\frac{1}{2{\sigma }^2}(\cdot -\mu {)}^2}$$

If $x\sim \mathcal{N}(\mu ,{\sigma }^2)$, then $E(x)=\mu$, $\mathrm{Var}(x)={\sigma }^2$.

Logistic sigmoid

$$s(x)=\frac{1}{1+e^{-x}}$$

Softplus function

$$\zeta (x)=\log (1+e^x)$$

Softmax

$$\sigma (x{)}_i=\frac{e^{x_i}}{{\sum }_j{e}^{x_j}}$$

$L^p$ norm

$$||x|{|}_p={\left(\sum _i|x_i{|}^p\right)}^{\frac{1}{p}}$$

XOR

Linear model cannot perform XOR operation

$$f(x,W,c,w,b)=w^T\cdot max(0,W^{T}x+c)+b$$

$$f\left(\left[\begin{array}{cc}0& 0\\ 0& 1\\ 1& 0\\ 1& 1\end{array}\right]\right)=\left[\begin{array}{c}0\\ 1\\ 1\\ 0\end{array}\right]$$

$$W=\left[\begin{array}{cc}1& 1\\ 1& 1\end{array}\right],c=\left[\begin{array}{c}0\\ -1\end{array}\right],w=\left[\begin{array}{c}1\\ -2\end{array}\right],b=0$$

General Problem

For a data-set with $N$ samples

$$(x_i,y_i)\in {\mathbf{R}}^n\times {\mathbf{R}}^m,i\in \mathbf{N},$$

$$f:{\mathbf{R}}^n\mapsto {\mathbf{R}}^m$$

$\hat{y}=f(x)$

$$||\hat{y}-y|{|}_2^2=\sum _i({\hat{y}}_i-y_i{)}^2$$

Linear Case

Module $x\in {\mathbf{R}}^n,$ $w\in {\mathbf{R}}^n\times {\mathbf{R}}^m,$ $y\in {\mathbf{R}}^m,$ $b\in {\mathbf{R}}^m,$

$$\hat{y}=f(x)=x\cdot w+b$$

Likelihood function

Likelihood is the probability that a particular outcome $x$ is observed when the true value of the parameter is $\theta$,

$$\mathcal{L}(\theta \mid x)=p_{\theta }(x)=P_{\theta }(X=x)$$

Unlike probabilities, likelihood function do not have to integrate (or sum) to 1.

Quadratic Problem

$$F(x)=\frac{1}{2}xw{x}^T-bx$$

$$\nabla F(x)=xw-b$$

$${\nabla }^{2}F(x)=w$$

$w$ is invertable,

$$x_1=x_0-\nabla F(x_0)({\nabla }^{2}F(x_0){)}^{-1}=x_0-(x_{0}w-b)w^{-1}=b{w}^{-1}=x^{\ast }$$

$$\hat{w}=\underset{w}{\mathrm{argmin}}S(w)$$

summary

• maximum likelihood estimation
• least square regression
• minimum cross-entropy between the distributions
• minimum the KL divergence

prevent overfitting

• maximum a posteriori estimation
• regularized least square

cross-entropy $\ne$ negative log-likelihood of a Bernoulli or softmax distribution

PCA

• maximize variance of data after projection

maximize variance -- lagrangian multiplier --> eigenvector of covariance matrix

Theorem (Bochner’s theorem)

A continuous function of the form 𝑘(𝑥,𝑦)=𝑘(𝑥−𝑦) is positive definite if and only if 𝑘(𝛿) is the Fourier transform of a non-negative measure.

Random Fourier features — Random walks

Entropy

Definition

For a discrete probability distribution $p$ on the finite set $\{x_1,x_2,\dots ,x_N\}$ with $p_i=p(x_i)$, the entropy of $p$ is defined as

$$h(p)=-\sum _{i=1}^N{p}_i\log p_i.$$

For a continuous probability density function $p$ on an interval $[a,b]$, the entropy of $p$ is defined as

$$h(p)=-{\int }_a^{b}p(x)\log p(x)dx.$$

Theorem

For a probability density function $p$ on a finite set $\{x_1,x_2,\dots ,x_N\}$, then

$$h(p)\le \log n,$$

with equality iff. $p$ is uniform, i.e. $\forall i\le N,p(x_i)=1/n$.

Uniform probability yields maximum uncertainty and therefore maximum entropy.

Theorem

For a continuous probability density function $p$ on $\mathbb{R}$ with variance ${\sigma }^2$, then

$$h(p)\le \frac{1}{2}(1+\log (2\pi {\sigma }^2))$$

with equality iff. $p$ if Gaussian with variance ${\sigma }^2$, i.e. for some $\mu$ we have

$$p(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}.$$

Theorem

For a continuous probability density function $p$ on $(0,\infty )$ with mean $\lambda$, then

$$h(p)\le 1+\lambda ,$$

with equality iff. $p$ is exponential with mean $\lambda$, i.e.

$$p(x)=\frac{1}{\lambda }{e}^{-\frac{1}{\lambda }x}.$$

Cross entropy

The cross entropy of the distribution $q$ relative to a distribution $p$ over a given set is defined as follows:

$$H(p,q)=-E_p[\log q]=-\sum p\log q=H(p)+D_{\mathrm{K}\mathrm{L}}(p\Vert q)$$

Kullback-Leibler divergence

The Kullback-Leibler divergence (relative entropy) was introduced as the directed divergence between two distributions

$$D_{\text{KL}}(P\parallel Q)=-\sum _{x\in \mathcal{X}}P(x)\log \left(\frac{Q(x)}{P(x)}\right)$$

The Kullback-Leibler divergence is then interpreted as the average difference of the number of bits required for encoding samples of $P$ using a code optimized for $Q$ rather than one optimized for $P$.

Jensen-Shannon divergence

$$D_{\text{JS}}(P\parallel Q)=\frac{1}{2}{D}_{\text{KL}}(P\parallel M)+\frac{1}{2}{D}_{\text{KL}}(Q\parallel M)$$

where $M=\frac{1}{2}(P+Q)$

Mutual Information

Let $(X,Y)$ be a pair of random variables with values over the space $\mathcal{X}\times \mathcal{Y}$. If their joint distribution is $P_{(X,Y)}$ and the marginal distributions are $P_X$ and $P_Y$, the mutual information is defined as

$$I(X;Y)=D_{\mathrm{K}\mathrm{L}}(P_{(X,Y)}\Vert P_X\otimes P_Y)$$

where $D_{\mathrm{K}\mathrm{L}}$ is the Kullback–Leibler divergence.

PMFs for discrete distributions

The mutual information of two jointly discrete random variables $X$ and $Y$ is calculated as a double sum:

$$I(X;Y)=\sum _{y\in \mathcal{Y}}\sum _{x\in \mathcal{X}}{P}_{(X,Y)}(x,y)\log \left(\frac{P_{(X,Y)}(x,y)}{P_X(x)P_Y(y)}\right),$$

where $P_{(X,Y)}$ is the joint probability mass function of $X$ and $Y$, and $P_X$ and $P_Y$ are the marginal probability mass functions of $X$ and $Y$ respectively.

PDFs for continuous distributions

In the case of jointly continuous random variables, the double sum is replaced by a double integral:

$$I(X;Y)={\int }_{\mathcal{Y}}{\int }_{\mathcal{X}}{P}_{(X,Y)}(x,y)\log \left(\frac{P_{(X,Y)}(x,y)}{P_X(x)P_Y(y)}\right)dxdy,$$

where $P_{(X,Y)}$ is now the joint probability density function of $X$ and $Y$, and $P_X$ and $P_Y$ are the marginal probability density functions of $X$ and $Y$ respectively.

Mutual information and Kullback–Leibler divergence

Mutual information is the Kullback–Leibler divergence from the product of the marginal distributions

$$I(X;Y)=D_{\text{KL}}\left(p_{(X,Y)}\parallel p_X{p}_Y\right)$$

$$I(X;Y)={\mathbb{E}}_Y\left[D_{\text{KL}}\left(p_{X\mid Y}\parallel p_X\right)\right]$$

Newton's Method

Given a twice differentiable function $f:\mathbb{R}\to \mathbb{R}$, solve the optimization problem

$$\underset{x\in \mathbb{R}}{\min }f(x).$$

Consider the second-order Taylor expansion of $f$ around $x_k$

$$f(x_k+t)\approx f(x_k)+f^{\prime }(x_k)t+\frac{1}{2}{f}^{\prime \prime }(x_k)t^2.$$

Setting the derivative to zero,

$$0=\frac{\mathrm{d}}{\mathrm{d}t}\left(f(x_k)+f^{\prime }(x_k)t+\frac{1}{2}{f}^{\prime \prime }(x_k)t^2\right)=f^{\prime }(x_k)+f^{\prime \prime }(x_k)t,$$

the minimum is achieved for

$$t=-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)},$$

then

$$x_{k+1}=x_k+t=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}.$$

Regression

Linear Regression

Objective function

$$\underset{w}{\min }(\Vert y-x\cdot w{\Vert }_2^2)$$

Ordinary least square Lagrangian function

$$L(w)=||y-x\cdot w|{|}_2^2=(y-x\cdot w{)}^T(y-x\cdot w)=y^{T}y-w^T{x}^{T}y-y^{T}xw+w^T{x}^{T}xw$$

$$\frac{\partial L(w)}{\partial w}=2x^{T}xw-2x^{T}y$$

Let $\partial L(w)/\partial w=0$ leads to

$$w=(x^{T}x{)}^{-1}{x}^{T}y$$

With respect to different constraints,

• $||w|{|}_1<t$ for lasso,
• $||w|{|}_2<t$ for ridge.

Ridge Regression : L1 Regularization

Least absolute shrinkage and selection operator

$$\underset{w}{\min }(\Vert y-x\cdot w{\Vert }_2^2+\lambda \Vert w{\Vert }_1).$$

• When $\lambda =0$ , no parameters are eliminated. The estimate is equal to the one found with linear regression.
• As $\lambda$ increases, more and more coefficients are set to zero and eliminated (theoretically, when $\lambda =\infty$, all coefficients are eliminated).
• As $\lambda$ increases, bias increases.
• As $\lambda$ decreases, variance increases.

Lasso Regression : L2 Regularization

$$\underset{w}{\min }(\Vert y-x\cdot w{\Vert }_2^2+\lambda \Vert w{\Vert }_2^2).$$

$$L(w)=w^T{x}^{T}xw-w^T{x}^{T}y-y^{T}xw+y^{T}y+\lambda w^{T}w$$

$$\frac{\partial L(w)}{\partial w}=2x^{T}xw-2x^{T}y+2\lambda w$$

Let $\partial L(w)/\partial w=0$ leads to

$$w=(x^{T}x+\lambda I{)}^{-1}{x}^{T}y$$

Elastic net Regularization

$$\underset{w}{\min }(\Vert y-x\cdot w{\Vert }_2^2+{\lambda }_2\Vert w{\Vert }_2^2+{\lambda }_1\Vert w{\Vert }_1).$$

Logistic Regression

For $p\in (0,1)$,

$$\mathrm{logit}p=\ln (\mathrm{odds})=\ln \frac{p}{1-p}$$

Logistic function is inverse logit, sigmoid function defined as follows, for $q\in (-\infty ,\infty )$,

$$p=s(q)={\mathrm{logit}}^{-1}q=\frac{e^q}{e^q+1}=\frac{1}{1+e^{-q}}$$

Suppose $\mathrm{logit}y$ linear with $x$, i.e.

$$\ln \frac{y}{1-y}={\beta }_0+{\beta }_{1}x$$

then

$$y=s(x)=\frac{1}{1+e^{-({\beta }_0+{\beta }_{1}x)}}$$

Conjugate Descent

Definition

Let $Q\in {\mathbb{R}}^{n\times n}$ be a symmetric matrix. $v_i,v_j\in {\mathbb{R}}^n,v_i\ne v_j\ne 0$ are Q-orthogonal w.r.t. $Q$ if

$$v_i^{T}Q{v}_j=0.$$

Proposition

Let $Q$ be positive definite. If $v_1,v_2,\dots ,v_k$ are Q-orthogonal, then they are linear independent.

Proof

If $v_k={\alpha }_1{v}_1+\cdots +{\alpha }_{k-1}{v}_{k-1}$, then $$v_k^{T}Q{v}_k=v_k^{T}Q({\alpha }_1{v}_1+\cdots +{\alpha }_{k-1}{v}_{k-1})={\alpha }_1{v}_k^{T}Q{v}_1+\cdots +{\alpha }_{k-1}{v}_k^{T}Q{v}_{k-1}=0$$ controversy with postive definite.

Proposition

Let $v_1,v_2,\dots ,v_n\in {\mathbb{R}}^n$ vectors Q-orthogonal, $x_0\in {\mathbb{R}}^n$ be arbitrary starting point.

Assume

$${\alpha }_k=-\frac{g_k^T{v}_k}{v_k^{T}Q{v}_k}$$

$$x_{k+1}=x_k+{\alpha }_k{v}_k$$

$$g_k=Q{x}_k-b$$

then $x_n$ is the solution of $Qx=b$.

$$\underset{x\in {\mathbb{R}}^n}{\min }\frac{1}{2}{x}^{T}Qx-b^{T}x$$

Matrix $Q\in {\mathbb{R}}^{n\times n}$ is positive definite.

$$Qx=b,x\in {\mathbb{R}}^n$$

Let $v_1,v_2,\dots ,v_n\in {\mathbb{R}}^n$ vectors Q-orthogonal, which is linear independent, then $\forall x^{\ast }\in {\mathbb{R}}^n$ there exist ${\alpha }_1,{\alpha }_2,\dots ,{\alpha }_n\in \mathbb{R}$,

$$x^{\ast }=\sum _{i=1}^n{\alpha }_i{v}_i={\alpha }_1{v}_1+{\alpha }_2{v}_2+\cdots +{\alpha }_n{v}_n$$

$$v_i^{T}Q{x}^{\ast }=v_i^{T}Q({\alpha }_1{v}_1+{\alpha }_2{v}_2+\cdots +{\alpha }_n{v}_n)={\alpha }_i{v}_i^{T}Q{v}_i$$

$${\alpha }_i=\frac{v_i^{T}Q{x}^{\ast }}{v_i^{T}Q{v}_i}=\frac{v_i^{T}b}{v_i^{T}Q{v}_i}$$

$$x^{\ast }=\sum _{i=1}^n{\alpha }_i{v}_i=\sum _{i=1}^n\frac{v_i^{T}b}{v_i^{T}Q{v}_i}{v}_i$$

Proposition

Let $x_0\in {\mathbb{R}}^n$ arbitrary,

$$\begin{array}{rl}{v}_0& =b-Q{x}_0\\ {\alpha }_k& =-\frac{g_k^T{v}_k}{v_k^{T}Q{v}_k}\\ x_{k+1}& =x_k+{\alpha }_k{v}_k\\ g_k& =Q{x}_k-b\\ v_{k+1}& =-g_{k+1}+{\beta }_k{v}_k\\ {\beta }_k& =\frac{g_{k+1}^{T}Q{v}_k}{v_k^{T}Q{v}_k}\end{array}$$

Gradient Descent

Multi-variable function $f:{\mathbb{R}}^n\mapsto \mathbb{R}$, defined differentiable in a neighborhood of a point $x\in {\mathbb{R}}^n$, for $\lambda \in {\mathbb{R}}_{+}$ small enough,

$$x_{k+1}=x_k-\lambda \nabla f(x_k)$$

leads to $f(x_{k+1})\le f(x_k)$.

If $f$ convex and $\nabla f$ Lipschitz, $f_k$ converge to a local mimimum.

Optimization : Momentum

Let $\gamma <1,v_0=0$,

$$\{\begin{array}{ll}{x}_{k+1}& =x_k-v_k\\ v_k& =\gamma v_{k-1}+\lambda \nabla f(x_k)\end{array}$$

The momentum term increases for dimensions whose gradients point in the same directions and reduces updates for dimensions whose gradients change directions. As a result, we gain faster convergence and reduced oscillation.

Optimization : Nesterov

Version with correction,

$$\{\begin{array}{ll}{x}_{k+1}& =x_k-v_k\\ v_k& =\gamma v_{k-1}+\lambda \nabla f(x_k-\gamma v_{k-1})\end{array}$$

This anticipatory update prevents us from going too fast and results in increased responsiveness, which has significantly increased the performance of RNNs on a number of tasks.

Adagrad

It adapts the learning rate to the parameters, performing larger updates for infrequent and smaller updates for frequent parameters. For this reason, it is well-suited for dealing with sparse data.

$$x_{k+1}=x_k-\frac{\lambda }{\sqrt{G_k+ϵ}}\nabla f(x_k)$$

Application : learned to recognize cats in Youtube videos; GloVe word embeddings.

Adadelta

Adadelta is an extension of Adagrad that seeks to reduce its aggressive, monotonically decreasing learning rate. Instead of accumulating all past squared gradients, Adadelta restricts the window of accumulated past gradients to some fixed size w.

RMSprop

RMSprop and Adadelta have both been developed independently around the same time stemming from the need to resolve Adagrad’s radically diminishing learning rates. RMSprop in fact is identical to the first update vector of Adadelta.

$$v(x_k,t):=\gamma v(x_k,t-1)+(1-\gamma )(\nabla f(x_k){)}^2$$

$$x_{k+1}=x_k-\frac{\eta }{\sqrt{v(x_k,t)}}\nabla f(x_k)$$

Adam

Adaptive Moment Estimation (Adam) is another method that computes adaptive learning rates for each parameter. In addition to storing an exponentially decaying average of past squared gradients $v_k$ like Adadelta and RMSprop, Adam also keeps an exponentially decaying average of past gradients $u_k$, similar to momentum.

$$\begin{array}{rl}{u}_k& ={\beta }_1{u}_{k-1}+(1-{\beta }_1)\nabla f(x_k)\\ v_k& ={\beta }_2{v}_{k-1}+(1-{\beta }_2){\nabla }^{2}f(x_k)\end{array}$$

$$\begin{array}{rl}{\hat{u}}_k& =\frac{u_k}{1-{\beta }_1}\\ {\hat{v}}_k& =\frac{v_k}{1-{\beta }_2}\end{array}$$

$$x_{k+1}=x_k-\frac{\eta }{\sqrt{{\hat{v}}_k}+ϵ}\nabla f(x_k)$$

Principal Component Analysis

In order to maximize variance, the first weight vector $w_1$ thus has to satisfy

$${\mathbf{w}}_1=\underset{\Vert \mathbf{w}\Vert =1}{\mathrm{argmax}}\left\{\sum _i{\left({\mathbf{x}}_{(i)}\cdot \mathbf{w}\right)}^2\right\}=\underset{\Vert \mathbf{w}\Vert =1}{\mathrm{argmax}}\{\Vert \mathbf{X}\mathbf{w}{\Vert }^2\}=\underset{\Vert \mathbf{w}\Vert =1}{\mathrm{argmax}}\left\{{\mathbf{w}}^T{\mathbf{X}}^{\mathbf{T}}\mathbf{X}\mathbf{w}\right\}$$

Since $w_1$ has been defined to be a unit vector, it equivalently also satisfies

$${\mathbf{w}}_1=\mathrm{argmax}\left\{\frac{{\mathbf{w}}^T{\mathbf{X}}^{\mathbf{T}}\mathbf{X}\mathbf{w}}{{\mathbf{w}}^T\mathbf{w}}\right\}$$

The quantity to be maximised can be recognised as a Rayleigh quotient. A standard result for a positive semidefinite matrix such as $X^{T}X$ is that the quotient's maximum possible value is the largest eigenvalue of the matrix, which occurs when w is the corresponding eigenvector.

Let $x\in {\mathbb{R}}^n$ be origin samples, projection vector $w$ satisfies

$$w^{T}w=1.$$

Denote $\Sigma =\frac{1}{N}(X-\bar{X})(X-\bar{X}{)}^T$ the variance of origin samples, variance of projected samples,

$$\sigma (X,w)=\frac{1}{N}(w^{T}X-w^T\bar{X})(w^{T}X-w^T\bar{X}{)}^T=w^T\Sigma w$$

Lagrangian multiplier form,

$$\max w^T\Sigma w+\lambda (1-w^{T}w),$$

leads to

$$\Sigma w=\lambda w.$$

the variance $\sigma (X,w)$ is maximized with $w$ is the eigenvector of $\Sigma$, the maximum value is the related eigenvalue $\lambda$.

PCA is either done by

• singular value decomposition of a design matrix;
• eigenvalue decomposition on the covariance matrix.

Eigenvalue Decomposition

A (non-zero) vector $v\in {\mathbb{R}}^n$ is an eigenvector of a square matrix $A\in {\mathbb{R}}^n\times {\mathbb{R}}^n$ if it satisfies the linear equation

$$Av=\lambda v$$

where $\lambda \in \mathbb{R}$ termed the eigenvalue corresponding to $v$.

Since $\forall i\in \mathbb{N}$,

$$A{v}_i=\lambda v_i,$$

$$A\left[\begin{array}{cccc}{v}_1& v_2& \cdots & v_n\end{array}\right]=\left[\begin{array}{cccc}{v}_1& v_2& \cdots & v_n\end{array}\right]\left[\begin{array}{ccccc}{\lambda }_1& & & & \\ & {\lambda }_2& & & \\ & & \ddots & & \\ & & & {\lambda }_n& \end{array}\right]$$

$$AV=V\Lambda$$

Diagonalisable matrice $A$ can be factorized as

$$A=V\Lambda V^{-1}$$

where ${\Lambda }_{ii}={\lambda }_i$ and $v_i$ is the $i$-th column of $V$.