为什么matlab的mldivide比dgels好那么多？答案

【问题标题】：Why is matlab's mldivide so much better than dgels?为什么matlab的mldivide比dgels好那么多？
【发布时间】：2019-09-30 16:06:53
【问题描述】：

解决Ax = b。真双。 A 用 M >> 2 超定 Mx2。b 是 Mx1。我针对mldivide 运行了大量数据，结果非常好。我用 MKL LAPACKE_dgels 编写了一个 mex 例程，但它远没有那么好。结果有大量噪音，基本信号几乎不存在。我首先根据 MKL 示例结果检查了例程。我搜索了mldivide doc（流程图）和 SO 问题。我发现的只是 Matlab 对超定矩形使用 QR 分解。

接下来我应该尝试什么？我是否使用了错误的 LAPACK 例程？请帮助指导我正确的方向。

更新： 对于解向量上的 E-15 浮点差，英特尔 MKL LAPACKE_dgels 具有与 Matlab mldivide 相同的结果，用于处理真正的双超定（矩形）问题。据我所知，这是使用的 QR 方法。

注意从这个 dgels 返回的残差。它们不等同于 b - Axe。其中许多接近这个值，而有些则远非如此。

【问题讨论】：

为了让我们检查您是否以错误的方式使用该函数，您最好提供一个最小的工作示例（即使使用小矩阵也不会重现该问题）。
作为替代方案，请尝试 LAPACK 的DGELSY()。它还利用 QR 分解来解决超定系统。最后，如果矩阵 A 是病态的，请尝试DGELSD()。见stackoverflow.com/questions/41637108/…和stackoverflow.com/questions/55367024/…
@francis。谢谢，这是个好建议。我读到犰狳使用 QR，如果它的条件不好，则回退到 SVD（dgelsd）。但是 dgels 的 int info 返回参数永远不会是正数（总是零）。
确实，文件include/armadillo_bits/glue_solve_meat.hpp 描述了对status = auxlib::solve_approx_fast() 的调用和对auxlib::solve_approx_svd() 的调用（如果状态为假）。 solve_approx_fast() 调用 lapack::gels 和 auxlib::solve_approx_svd() 调用 lapack::gelsd()。如果找到空枢轴，则 dgels 的返回参数变为正数，从而阻止系统被求解。它始终为零的事实与矩阵的良好调节一致。恭喜您自己解决了问题！

标签： matlab linear-algebra lapack intel-mkl

【解决方案1】：

问题不是x 的解，而是DGELS 返回的残差。该例程的输出在输入数组指针上就地修改。 MKL doc 表示输入数组 b 被第一行 N 的输出向量 x 覆盖，然后是 N+1 到 M 中的残差。我用我的代码确认了这一点。

错误在于将b[N+1] 残差与原始输入b[1] 对齐，并对此做出进一步的算法决策。残差与原始输入的正确对齐是b[1] 到b[1]。第一个N残差不可用；你必须在之后计算这些。

文档并没有说它们本身是残差，而是特别说明

每列中解的残差平方和由该列中元素n+1 到m 的模数平方和给出。

【讨论】：