为什么Matlab的mldivide比dgels好得多? [英] Why is matlab's mldivide so much better than dgels?

问题描述

解决Ax = b.真正的双重. A用M >> 2超确定为Mx2.b是Mx1.我已经对mldivide运行了大量数据，结果非常好.我用MKL LAPACKE_dgels编写了一个mex例程，但效果却差强人意.结果会产生大量噪声，而基本信号几乎不存在.我首先对照MKL示例结果检查了例程.我已经搜索了mldivide doc(流程图)和SO问题.我发现的是Matlab使用QR分解处理过高的矩形.

Solve Ax = b. Real double. A is overdetermined Mx2 with M >> 2. b is Mx1. I've run a ton of data against mldivide, and the results are excellent. I wrote a mex routine with MKL LAPACKE_dgels and it's nowhere near as good. The results have a ton of noise and the underlying signal is barely there. I checked the routine against the MKL example results first. I've searched through the mldivide doc (flowchart) and the SO questions. All I found is Matlab uses QR factorization for overdetermined rectangular.

接下来我应该尝试什么?我使用了错误的LAPACK例程吗?请帮助指引我正确的方向.

What should I try next? Am I using the wrong LAPACK routine? Please help guide me in the right direction.

更新: 对于求解矢量上E-15浮点数之差，英特尔MKL LAPACKE_dgels对于真实的双重超定(矩形)问题具有与Matlab mldivide相同的结果.据我所知，这是使用的QR方法.

Update: To within E-15 floating point difference on the solution vector, Intel MKL LAPACKE_dgels has the same result as Matlab mldivide for real double overdetermined (rectangular) problems. As far as I can tell, this is the QR method used.

当心从这些dgel返回的残差.它们不等于b-Ax.他们中的许多人都接近这个值，而有些人则离它很远.

Beware the residuals returned from this dgels. They do not equate to b - Ax. Many of them are close to this value, whereas some are far from it.

推荐答案

问题不是解决方案x，而是问题DGELS返回的残差.该例程的输出在输入数组指针上就地修改. MKL文档表示输入数组b被输出矢量覆盖x用于前N行，然后是N+1至M中的残差.我用我的代码确认了这一点.

The problem was not the solution x, rather the returned residuals from DGELS. This routine's outputs are modify-in-place on the input array pointers. The MKL doc says the input array b is overwritten with the output vector x for the first N rows, then the residuals in N+1 to M. I confirmed this with my code.

错误在于将b[N+1]残差与原始输入b[1]对齐，并对此做出进一步的算法决策.残差与原始输入的正确对齐方式是b[1]至b[1].前N个残差不可用；您必须随后进行计算.

The mistake was in aligning the b[N+1] residuals to original inputs b[1], and making further algorithmic decisions on that. The correct alignment of residual to original input is b[1] to b[1]. The first N residuals are not available; you have to compute those afterwards.

文档并没有说它们本身就是残差，而是具体来说

The doc doesn't say they are residuals per se, rather specifically

每列中解决方案的剩余平方和由该列中元素n+1至m的模量平方和得出.

the residual sum of squares for the solution in each column is given by the sum of squares of modulus of elements n+1 to m in that column.

这篇关于为什么Matlab的mldivide比dgels好得多?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！