Lecture3 More MatMul and the Roofline Performance Model

Review

Computational Intensity: 非常重要

t m / t f t_m/t_f tm​/tf​ 被称为 Machine Balance，是机器效率的关键。

Blocked (Tiled) Matrix Multiply

可以发现按照如下流程计算，最终得到的 CI 差不多就是每个矩阵块的大小 b（如果实现正确的话）。

因此，为了让 CI 更大，我们需要更大的 b。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-9TrPyBj7-1678604023050)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230305205404264.png)]

我们假设快速内存大小为 M。则一定有 b ≤ M / 3 b\le\sqrt{M/3} b≤M/3 ​，因为 M 需要装上 3 个 b x b 的块。

Is there a more elegant approach?

Recursive Matrix Multiplication

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-rwCZMkX2-1678604023051)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230305210031389.png)]

Define C = RMM(A, B, n)
if (n = 1) {
	C = A * B;
} else {
	C00 = RMM(A00, B00, n/2) + RMM(A01, B10, n/2);
	C01 = RMM(A00, B01, n/2) + RMM(A01, B11, n/2);
	C10 = RMM(A10, B00, n/2) + RMM(A11, B10, n/2);
	C11 = RMM(A10, B01, n/2) + RMM(A11, B11, n/2);
}

定义 Arith(n) = number arithmetic operations in RMM(…n)

所以 Arith(n) = 8 x Arith(n/2) + 4(n/2)^2 if n > 1, else 1

= 2 n 3 − n 2 2n^3-n^2 2n3−n2

定义 W(n) = number of words moved between fast, slow memory by RMM(…n)

所以 W(n) = 8 x W(n / 2) + 4 x 3(n / 2)^2 if 3n^2 > M_fast else 3n^2

= O ( n 3 / ( M f a s t ) 1 / 2 + n 2 ) O(n^3/(M_{fast})^{1/2}+n^2) O(n3/(Mfast​)1/2+n2)

不论是 Arith 还是 W 都和 blocked 矩阵乘法相似，但是如果有多级 cache。

Cache Oblivious

在实践中，会在 1x1 的块运算之前剪枝。

  • 称为在小矩阵块上的 micro-kernel

Pingali et al 展示了2/3 左右的性能峰值

  • 使用 Recursive + optimized micro-kernel
  • PPT - A Comparison of Cache-conscious and Cache-oblivious Programs PowerPoint Presentation - ID:2922418 (slideserve)

有很多优化可以做。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-yaowQSZF-1678604023051)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230305213417516.png)]

Alternate Data Layouts

如果想要最小化内存移动，考虑到最快的方式获取连续内存，否则就会造成 latency 问题。如果想要最大化 bandwidth，最好从 slow memory 中连续地取内存。

所以有时候也要使用 copy optimization。就是在程序开始之前，我们对于内存重新规划，使其和原来的顺序不同。

有两种方法，一种是分块之后按照 b x b 矩阵块 row - major进行（适合两级内存结构，如果有多级内存结构呢？）

另外一种方法是递归的 Z-Morton 顺序。它适用于任何 cache 大小，但是 index calculation 非常 expensive。并且当逐渐递归下去到最后可能变成 row/col major。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-qbwisCRw-1678604023052)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230305214525558.png)]

Theory: Communication Lower bound

Theorem (Hong & Kung, 1981):

Any reorganization of matmul (using only commutativity and associativity) has computational intensity q = O ( ( M f a s t ) 1 / 2 ) q = O((M_{fast})^{1/2}) q=O((Mfast​)1/2), so # words moved between fast/slow memory = Ω ( n 3 / ( M f a s t ) 1 / 2 ) \Omega(n^3/(M_{fast})^{1/2}) Ω(n3/(Mfast​)1/2)

这里做了一个课程预告：

Extensions, both lower bounds and (some) optimal algorithms (later in course)

  • 之后会讨论 parallel matrix multiply，优化 latency 和 bandwidth。最终达到一个性能下界，还有 latency 下界。
  • 同样适用于 Rest of linear algebra (Gaussian elimination 高斯消元法, least squares, tensors contraction 有 lower bound 和 optimal 算法)
  • 之后可以将这些理论扩展到比线性代数更加通用的代码，比如 Nested loops accessing arrays (e.g. All-Pairs-Shortest-Paths, N-body, …) —— 访问数组的嵌套循环
  • 开放问题：
    • Small loop bounds (e.g. matrix-vector vs matrix-matrix multiply) 矩阵向量乘法是矩阵乘法的一个特例，当矩阵越来越细，算法如何进行拓展？
    • Dependencies, i.e. when only some reorganization are correct —— 在矩阵乘法中，交换 i j k 的顺序并不会导致结果发生错误，但是在一些问题中，不同的执行顺序可能导致错误结果。

Strassen’s Matrix Multiply

一个在 60 年代发现的结论。矩阵乘法可以不必是 O ( n 3 ) O(n^3) O(n3) 的，而是可以 O ( n 2.81 ) = O ( n log ⁡ 7 ) O(n^{2.81}) =O(n^{\log 7}) O(n2.81)=O(nlog7)

考虑一个 2 x 2 的矩阵乘法，传统的进行 8 次乘法，4 次加法

  • Strassen’s 进行了 7 次乘法，18 次加法

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-iKnzIqnk-1678604023052)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230305221555092.png)]

T ( n ) = 7 ∗ T ( n / 2 ) + 18 ∗ ( n / 2 ) 2 = O ( n log ⁡ 7 ) = O ( n 2.81 ) T(n) = 7*T(n/2) + 18*(n/2)^2 = O(n^{\log 7})=O(n^{2.81}) T(n)=7∗T(n/2)+18∗(n/2)2=O(nlog7)=O(n2.81)

当 n 很大时，这种方法会更快一些（也取决于机器）(“Tuning Strassen’s Matrix Multiplication for Memory Efficiency”)

  • Possible to extend communication lower bound to Stassen
    • #words moved between fast and slow memory = Ω ( n log ⁡ 7 / M log ⁡ 7 / 2 − 1 ) ∼ Ω ( n 2.81 / M 0.4 ) \Omega(n^{\log 7}/M^{\log 7 / 2-1})\sim \Omega(n^{2.81}/M^{0.4}) Ω(nlog7/Mlog7/2−1)∼Ω(n2.81/M0.4) (Ballard D, Holtz, Schwartz, 2011, SPAA Best Paper Prize)
    • Attainable too, more on parallel version later

Other Fast Matrix Multiplication Algorithms

  • 世界记录是 O ( n 2.37548... ) O(n^{2.37548...}) O(n2.37548...) Coppersmith & Winograd, 1987
  • 后来减少到了 2.3728642 (Virginia Vassilevska Williams, UC Berkeley & Stanford, 2011)
  • 从 2.37293 减少到了 2.3728639 (Francois Le Gall, 2014)
  • Latest Record! 2.3728639 reduced to 2.3728596 (Virginia Vassilevska Williams and Josh Alman, 2020)

Basic Linear Algebra Subroutine (BLAS)

  • 网站：BLAS (Basic Linear Algebra Subprograms) (netlib)，BLAS Techical Forum (netlib)

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-LMKiJL6I-1678604023052)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230306075715647.png)]

**Q：**我听说 Strassen 方法精度不高，他还能被广泛应用吗？

**A：**精度不是限制这个方法推广的问题，其实使用半精度浮点数已经表明不需要精确结果了。这个方法没有被广泛使用的原因可能是政治方面的，LINPACK 要求不使用 Strassen 方法，因此许多供应商主要优化立方矩阵方法，而不是 Strassen 方法。

Roofline Model

ChatGPT 介绍什么是 Roofline Model：

Roofline model是一种性能建模方法，用于评估计算机系统的性能限制。它是由David Patterson等人在2009年提出的。该模型提供了一种图形化方式来可视化内存带宽和计算性能之间的关系，并确定系统在特定工作负载下的性能瓶颈。

Roofline model使用两条线来表示系统性能的限制。第一条线是表示计算能力的峰值，通常是浮点运算速度的峰值。第二条线是表示内存带宽的峰值，通常是系统内存总带宽的一部分。

Roofline图形显示了工作负载的计算密集度和内存访问模式，以及这两条限制线。通过将工作负载的性能指标表示为数据点绘制到图形上，可以确定系统的性能瓶颈是计算密集度还是内存带宽。

Roofline model是一种有用的性能建模方法，可以帮助开发人员优化他们的代码以提高系统性能。

我们为什么有性能模型？

现在我们有很多的评价方法：Running time、Bandwidth、Memory footprint、Energy Use、Percent of Peak 等。

那么有这种性能模型有什么意义呢？这是因为当我们采用算法并进行优化时，我们想知道什么是正确的优化。我们希望一些指标帮助我们判断是否需要投入大量时间。

另外 Roofline Model 可以为程序每个部分提供不同的性能瓶颈（评估每个部分的用时，针对性的进行优化，知道程序的瓶颈改变）。

History of the Roofline Model

Samuel Williams, Andrew Waterman, David Patterson, “Roofline: an insightful visual performance model for multicore architectures.”

Roofline

Idea：应用同时被计算峰值和内存带宽限制。一些例子：

  • Bandwidth bound（matvec）
  • Compute bound（matmul）

What’s in the Roofline Model?

主要包含三个部分，其中两个是针对机器的，1 个针对应用。

针对机器：

  • 浮点数运算性能 (flops/sec)

    • 主要取决于时钟速度和并行度（ILP, SIMD, Multicore）

  • 内存带宽 (bytes/sec)

    • 不包括延迟 Latency

针对应用：

  • Computation Intensity
    • 将直接计算计算强度，而不是使用一些其他优化算法，我们将只计算你的算法的 flops 次数，然后计算输入输出大小，获得计算强度的上限。
    • 如果结果很高的话，可能告诉我们 tiling 算法可能有用。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-lS1kja64-1678604023052)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230306091006925.png)]

这是 Roofline 基础图表。横轴上是计算强度，以 DRAM 每字节的浮点数来衡量（这是对数刻度的）。如果是 1 表示缓存和 DRAM 之间移动一个字节所需的时间内完成一次 flops。

在横轴上将会根据我们评估的算法的计算强度选择合适的位置。

纵轴上是 attainable Gflops/s，同样是对数刻度的。表示性能峰值 roof

Different Roofs

接下来介绍四种不同的 roof。

假如你在没有任何其他类型的限制的情况下，每秒执行 64 Gflops，这是你的电脑的计算峰值。但是你不一定能达到。

如果乘法和加法的数量相同，大多数机器可以使用 FMA(Fuse Multiply Add) 获得计算峰值。而如果你不使用 FMA，就会导致性能的损失（例如一半）

假设你可以执行 SIMD，那么如果你不执行 SIMD，则性能 roof 又会下降

假如你机器可以使用 ILP，那么如果你不使用 ILP，则性能 roof 接着下降。

What is a better model?

  • 矩阵-矩阵乘法比矩阵-向量乘法更快，体现在 roofline 图上 dgemm 更高，但是仅仅知道计算性能峰值，不能解释为什么矩阵-矩阵乘法更快。
  • 对于矩阵-向量乘法，最好的情况是每移动一个 word，进行两次浮点运算（考虑了乘法和加法）。
  • 2 flops / word = 1/2 flops per byte (single precision) = 1/4 for double pricision

Data Movement Complexity

  • Assume run time ~= dat moved to/from DRAM

  • 假设我们性能收到了数据移动的限制。那么应该如何估计呢？

  • 强制性的数据移动（数据结构大小）are good first guess

  • Roofline 为我们提供的是一些问题的最大上限（一定无法超过，很多时候难以接近）

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-avi1HjgI-1678604023053)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230306100228646.png)]

Machine Balance and Computational Intensity

  • Machine balance is: balance = (Peak DP FLOP/s) / Peak Bandwidth
  • 一些机器的 Balance：
    • Haswell 10 Flops / Bytes
    • KNL is 34 Flops/Bytes to DRAM or 7 Flops/Bytes to HBM(High Bandwidth Memory)
  • Computational / Arithmetic Intensity (CI/AI) is
    • CI = FLOPs Performed / Data Moved
  • 一些算法的 CI 值

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Y61d0jx7-1678604023053)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230306100838525.png)]

（DRAM）Roofline

Roofline 的另外一个部分的限制

假设：

  • 理想的处理器和 cache
  • 冷启动（数据都在 DRAM 中）

接下来要计算 max time（分别 time spent on flops 和 time spent on data movement）：

为什么计算 Max 而不是 Sum time？

反映出这些可能并行运行的事实，我可以将内存交换和计算部分重叠

接下来通过将浮点运算数量除以时间计算每秒 Flops，由于分母都是 time，所以可以得到两者的最小值：

相当于：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-WwMtcVUP-1678604023054)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230306163902410.png)]

如果受到了运算峰值限制，则 GFlops/s 更小，否则受到了内存移动限制（下面的公式是 Computation Intensity 乘以 Peak GB/s）

这里的 #FP ops / Time 表示实际中总时间平摊下来进行的浮点数（浮点数/总时间）。这个值对应实际的纵轴。

这个 min 的式子表示实际的浮点数/总时间取决于两个部分：

  • 峰值的计算性能
  • 运算强度（总浮点数/总内存移动数）与峰值内存移动速率的乘机

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-ROpNRkN7-1678604023054)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230306164112925.png)]

画出图就和上面长得比较类似。如果我们计算一下两个部分的分界点，可以发现就是峰值运算速度/峰值内存移动速度，也就是之前说的机器的 balance。

当实际的计算强度处于 Balance 的左边时，表明整体受到了内存移动的限制，即 Bandwidth Bound；当实际计算强度处于 Balance 的右边时，受到了最大峰值运算速度的限制。

The Roofline Performance Model

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-vVGIQcPc-1678604023054)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230306165424240.png)]

从这张图可以看到，针对不同的机器，有不同的计算峰值（是否支持 ILP、是否支持 SIMD 等）；同时也有不同倾斜的部分，例如是否使用 NUMA（分布式内存架构，每个核有一个和自己靠的近的内存区域，当访问这块内存区域时速度很快，访问和其他核靠的近的内存区域时速度很慢）、是否支持预取内存等优化。

Roofline Across Algorithms

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xk1Bxtqs-1678604023054)(D:\zhang_kg\BUAA_undergraduate\TyporaImage\image-20230306170718677.png)]

这是一些算法在不同机器下的表现情况

  • 在左边 Xeon X5550 上面取决于是否使用单精度，可以发现峰值计算速度是使用双精度的两倍。
  • 最快的是 DGEMM，并且呈现出一个斑点（blob），因为计算强度有变化。

/gs/home/liaojianjin/Zhang_kg/lapack-3.11.0

export LIBRARY_PATH="$LIBRARY_PATH:/gs/home/liaojianjin/Zhang_kg/lapack-3.11.0"
export C_INCLUDE_PATH="$C_INCLUDE_PATH:/gs/home/liaojianjin/Zhang_kg/lapack-3.11.0/LAPACKE/include:/gs/home/liaojianjin/Zhang_kg/lapack-3.11.0/CBLAS/include"
/liaojianjin/Zhang_kg/lapack-3.11.0

```powershell
export LIBRARY_PATH="$LIBRARY_PATH:/gs/home/liaojianjin/Zhang_kg/lapack-3.11.0"
export C_INCLUDE_PATH="$C_INCLUDE_PATH:/gs/home/liaojianjin/Zhang_kg/lapack-3.11.0/LAPACKE/include:/gs/home/liaojianjin/Zhang_kg/lapack-3.11.0/CBLAS/include"

更多推荐

『UC Berkeley CS267』Lecture3: More MatMul and the Roofline Performance Model