LAPACK 是否落后于密集线性代数的前沿?

计算科学 线性代数 特征值 scipy 拉帕克 朱莉娅
2021-12-04 13:04:11

我最近一直在研究一些数值线性代数,并特别阅读LAPACK 如何解决对称特征值问题我注意到*stegr计算对称三对角矩阵的特征值和特征向量的例程(更一般的情况被简化)是O(n2)及时,如参考文献中所述(有趣的是,最后一个 URL 具有后缀“#holygrail”)。LAPACK 文档(与第一个链接相同)指出该算法(相对稳健的表示算法)“除了少数情况外,比所有其他例程都快,并且使用的工作空间最少。”

然而,我在三对角矩阵的维基百科文章中注意到存在计算实对称三对角矩阵特征值的O(nlogn)

给出的算法的参考来自 2012 年。根据这个表,看起来 LAPACK 的最后一个版本是在 2000 年。那么 LAPACK 是不是明显不理想的情况在其性能方面?如果是这样,为什么这么多包(如Julia 的 LinearAlgebraPython 的 SciPy)在后台使用 LAPACK?O(nlogn)

(请注意,我看到了这篇文章,这表明对称三对角矩阵的算法无助于密集矩阵的整体渐近性能,因为三对角化步骤更昂贵。但我上面的内容仍然促使我想知道:LAPACK 是否有已知的方式不能在很大程度上表现得那么好?)O(nlogn)

3个回答

当有人说算法的顺序为时,这可能意味着复杂度由以下公式给出:添加的每个新元素都会增加运行时间(有效地)。有数学头脑的人经常忘记的是,这些陈述不包括常数的大小。这当然会延续到等。我无法绝对回答您的问题,是否存在其他库更快的特殊情况,但遗憾的是复杂性顺序只是故事的一半。当你设计一个算法时,你经常需要权衡,允许你将工作从线性项推到常数项,即在进行元素计算之前进行某些预计算。O(n)c+bnO(n²)针对给定问题优化订单并不总是会改善您的运行时间。盈亏平衡点可以任意大。这里唯一正确的方法是测量。

LAPACK 已经走在最前沿大约三十年了,而且可能仍然是为了它的利基市场。然而,鉴于 LAPACK 传统上构建的更简单的 BLAS 类型矩阵运算库的最新发展我们可以想象在未来几年内我们会看到传统的基于 FORTRAN 的 LAPACK 的严重竞争对手的出现。

换句话说,BLIS (C)、MKL (C/C++) 和Octavian.jl (Julia) 等库似乎几十年来第一次击败了传统的 BLAS,这表明确实仍有性能在桌面上对于这样的基本操作。

为了比较(通过 Octavian.jl 存储库): BLAS 类型的库基准

参考 LAPACK 库优先考虑可移植性而不是性能,因此它缺少许多优化,例如矢量化和线程化。(此外,它的大部分性能都依赖于 BLAS,因此使用参考 BLAS 会显着限制性能。)

但是,大多数 LAPACK 用户不使用参考库,而是各种优化的 BLAS+LAPACK 库,提供相同的 API,但优化各种例程。示例包括 Intel 的 MKL、IBM 的 ESSL、OpenBLAS、Cray 的 libsci、ARM 的性能库和 CUDA 的 cuBLAS。此共享 API 允许编写一次代码,同时能够使用特定系统上可用的最佳 LA 实现。

Julia 默认使用 OpenBLAS,而 SciPy 似乎让用户控制他们对实现的选择。