Matlab[Q,R,E]=qr(A)后面对应的LAPACK函数是什么？

计算科学 matlab 拉帕克
2021-12-21 02:25:37

我目前正在尝试廉价地计算矩阵的良好秩估计A. 因此，我使用计算柱旋转 QR 分解

[Q,R,E]=qr(A)

在 Matlab 中。我估计排名A使用

tol = size(A,n)*eps*norm(A,'fro'); 
r = sum(abs(diag(R))>tol)

这工作正常，并且 R 的所有对角线条目的图如下所示： 情节（排序（abs（diag（R）），1，'下降'），'r *'）

如果将整个算法移植到 C/Fortran，我会使用 LAPACK 中的 DGEQP3 替换 [Q,R,E]=qr(A)，它还计算旋转 QR 分解的列。但是，如果我对排名使用相同的估计，我大多会出错。相同的情节R从 DGEQP3 产生的看起来像

两个实验的输入矩阵完全相同。

我现在的问题是 Matlab 中旋转 QR 分解的列依赖于哪个 LAPACK 函数？

感谢您的帮助，格里苏

编辑： DGEQPF 给出了同样的错误结果。

编辑2：

  • 输入矩阵A是密集的并且构建为E+sign(E,F)
  • A可在此处获得：http ://www-e.uni-magdeburg.de/makoehle/A.mtx.gz（MatrixMarket 格式）
  • 错误的 R: http://www-e.uni-magdeburg.de/makoehle/R_wrong.mtx.gz
  • 我将 LAPACK 3.4.0 与 OpenBlas/GotoBLAS（64 位）一起使用
  • Matlab 7、2007b、2010b Linux 32 位

Edit3: - 使用 GDB 我发现，Matlab 2010b 调用 DGEQP3: #3 0xaa46ce2f in dgeqp3_ () from /usr/ubuntu10.04/matlabr2010b/bin/glnx86/../../bin/glnx86/../。 ./bin/glnx86/mllapack.so 为什么我使用 LAPACK（3.4.0 包含工作说明 176 中提到的修复）得到错误的结果？

2个回答

这里有两个问题：

  • 是A密集还是稀疏？
  • 您是否拥有与 MATLAB 的内部库相同的软件堆栈？

密集还是稀疏？

MATLAB 不再明确提及它调用的 LAPACK 例程以获得 QR 分解，如果A是稠密的。如果MATLAB R2008b 文档中的信息也适用于以后的版本，则DGEQP3当您调用[Q,R,E] = qr(A). 如果A是稀疏的，然后 MATLAB 调用Tim Davis 组中的SuiteSparseQR，该组与 SuiteSparse 库中的 UMFPACK 捆绑在一起。

您是否拥有与 MATLAB 的内部库相同的软件堆栈？

可能不是，这可能是您得到不同结果的原因之一。

我在对我正在编写的使用 QR 分解的库进行单元测试时遇到了这个问题。我使用 MATLAB 对我的工作进行原型设计，得到的结果与使用 LAPACK 或 NumPy 不同。据我所知，由于 MathWorks 无法轻松找到这些信息，MATLAB 使用不早于 3.1.1 版的 LAPACK 版本，以及英特尔的 MKL BLAS 库（适用于 Windows、英特尔 Mac 和 Linux）9.1 版或更高（见这里）。我找不到有关 SuiteSparse MATLAB 使用的版本的任何信息。通过在线挖掘或查看系统的库文件，您可能能够收集到更多信息。您可以尝试更改 MATLAB 链接到的库，以便能够跨软件包与相同的库进行比较；Eric Chu 提供了一篇不错的文章这至少显示了如何用自己的替换 MATLAB 的 BLAS 库（当然，这样做的风险由您自己承担）。他建议你也可以用 LAPACK 做同样的事情。甚至可以将 MATLAB 使用的 SuiteSparse 版本替换为您自己的版本。

LAPACK 版本会发生变化，BLAS 版本也会发生变化；他们可能使用不同版本的算法或不同的排序约定，即使A=QR和Q正交和R是上三角形，与版本无关。这些变化使再现性成为挑战。甚至您用于确定数字排名的容差也是一种判断。您似乎正在使用标准公差。

我最终使用 NumPy 对我的 QR 分解结果进行原型制作，因为它使用系统 BLAS 和 LAPACK 库。NumPy 和 SciPy 不是 MATLAB 的直接替代品，因为这两个库结合起来缺少 MATLAB 的一些功能，但是对于这个特定的线性代数任务，Python + NumPy + SciPy + Matplotlib 应该可以很好地工作。

请参阅 Leslie Foster关于排名显示软件的页面。另请参阅此LAPACK 工作笔记分析排名显示 QR 的失败xGEQP3。

您应该能够通过在调试器中设置断点并检查堆栈来找出 MATLAB 使用的例程。最后，我承认，几年前，MATLAB 使用了共享库，在这种情况下，符号名称不能被剥离，所以你会在调用堆栈上看到函数名称（但看不到参数，因为它肯定没有保留调试信息）。

其它你可能感兴趣的问题
上一篇查找成本与计算成本 下一篇混合整数线性规划的有效解