八度不准确地求解线性方程组？答案

【问题标题】：Octave inaccurately solving systems of linear equations?八度不准确地求解线性方程组？
【发布时间】：2017-06-14 09:35:43
【问题描述】：

我发现 Octave 无法准确地求解方程组。

k= [1.0000e+015,  -2.6667e-001,  3.3333e-002;-2.6667e-001,  5.3333e-001,  -2.6667e-001;3.3333e-002,  -2.6667e-001,  2.3333e-001];
f = [0;6.66667;41.66667];
x = (k\f)
x1 = inv(k)*f

x = (k\f)给了

-3.8893e-031
-1.9136e+001
 1.1192e+001

x1 = inv(k)*f给了

 4.8339e-014
 2.3752e+002
 4.5003e+002

两个答案之间存在差异。只是为了检查其他软件是否可以解决它，我在python numpy中尝试过。

import numpy as np
k = np.array ([[1.0000e+015,  -2.6667e-001,  3.3333e-002],[-2.6667e-001,  5.3333e-001,  -2.6667e-001],[ 3.3333e-002,  -2.6667e-001,  2.3333e-001]])
f = np.array([[0],[6.66667],[41.66667]])
x = np.linalg.solve(k, f)
print(x)

[[  4.83388655e-14]
 [  2.37521517e+02]
 [  4.50034428e+02]]

它可以正确解决。我想在倍频程中使用运算符\ 来求解方程（我想这是 MATLAB 中推荐的方法）。有什么建议吗？

在 Windows 上使用 Octave 4.2.x。

编辑：

为什么答案被否决？这是 GNU Octave 所说的 https://www.gnu.org/software/octave/:

与许多 Matlab 脚本兼容的插件

尝试在 MATLAB 中运行相同的脚本并自己思考。

评论者归咎于病态矩阵。一个好的求解器应该能够以足够的精度求解（就像 numpy 和 MATLAB 一样）。

【问题讨论】：

我不这么认为。它特定于八度运算符 \ 并且可能是该运算符的算法。其他方法有效，但不推荐。如果浮点数有问题，这两种方法都不行。
您在尝试这些时是否碰巧注意到警告？ “警告：机器精度的矩阵奇异，rcond = 6.66612e-17”
是的，不过numpy 给出了更准确的答案。事实上，八度音程inv 也有作用。查看 octave \ 和 inv 中的错误幅度，这两种方法都会发出警告。
这是一个条件非常差的矩阵。遇到重大错误也就不足为奇了。
我认为k\f 与inv(k)*f 不同，而是在数值上更稳定的伪逆，即inv(k'k)*(k'f)：这就是结果不同的原因。

标签： matlab numpy octave

【解决方案1】：

总结：Octave 开发人员选择返回（标记的）奇异方程组的“最小范数解”。

================

k的第一个元素比其他所有元素都大，而f的对应项为0。在x1解决方案中

4.8339e-014
2.3752e+002
4.5003e+002

第一项实际上是 0。

对剩余的 2x2 矩阵执行求解，产生相同的值：

>> inv(k(2:end,2:end))*f(2:end)
ans =    
   237.52
   450.03

>> k(2:end,2:end)\f(2:end)
ans =
   237.52
   450.03

在numpy:

In [518]: np.linalg.solve(k[1:,1:],f[1:])
Out[518]: 
array([[ 237.52151737],
       [ 450.0344278 ]])

==================

由于某种原因，\ 求解器对矩阵奇异性更敏感。

>> k\f
warning: matrix singular to machine precision, rcond = 6.66612e-17
warning: matrix singular to machine precision, rcond = 6.66612e-17

是的，测试这个解决方案表明它已经完成了

>> k*(k\f)
ans =
    5.4761
  -13.1904
    7.7144

>> rcond(k)
ans =    6.6661e-17
>> rcond(k(2:end,2:end))
ans =  0.083327

我不确定这有什么意义，但 numpy lstsqr 求解器返回与 k\f 相同的错误值：

In [525]: np.linalg.lstsq(k,f)
Out[525]: 
(array([[ -3.88925778e-31],
        [ -1.91360857e+01],
        [  1.11917602e+01]]),
 array([], dtype=float64),
 2,
 array([  1.00000000e+15,   6.89292234e-01,   7.73677656e-02]))

>> pinv(k)*f  #  Moore-Penrose pseudoinverse
ans =
  -3.8893e-31
  -1.9136e+01
   1.1192e+01

http://savannah.gnu.org/bugs/?42181 - 声称已记录在案。对于单数、欠定或超定方程，Octave 返回“最小范数解”，而不是“基本解”（MATLAB 的选择）。

http://wiki.octave.org/FAQ#Solvers_for_singular.2C_under-_and_over-determined_matrices

【讨论】：