2

我已经编写了 Gauss Seidel 和 Conjugate Gradient 迭代算法来解决 Haskell 中的矩阵（但这个问题与方法有关，而不是语言）。我的理解是，这两种算法都应该具有相似的收敛特性，并且在大多数情况下，CG 方法应该更快。我已经对来自http://math.nist.gov/MatrixMarket/的对称正定矩阵进行了许多测试而且我几乎永远无法获得CG算法。收敛，而 GS 几乎总是这样。我找不到任何带有右侧向量的对称正定矩阵用于在线测试目的，所以我一直在任意创建自己的 RHS（也许这是问题的一部分？）。如果我在 Ax = b 中使用 (transpose A) * A 而不是 A，我可以让 CG 方法收敛，这只是迫使矩阵对称。我在这里包含了 CG 代码。它显然不会按原样编译。如果有人需要它来提供帮助，我会全部发布。对于来自（伪代码和示例）的简单示例（类似问题） ，它可以正常工作. 关于共轭梯度与高斯赛德尔收敛标准，我有什么遗漏吗？谁能指出我正确的方向来完成这项工作？谢谢。 

conjGrad :: (Floating a, Ord a, Show a) => a -> SpMCR a -> SpVCR a -> SpVCR a -> (SpVCR a, Int)
conjGrad tol mA b x0 = loop x0 r0 r0 rs0 1
  where r0  = b - (mulMV mA x0)        
        rs0 = dot r0 r0
        loop x r p rs i
          | (varLog "residual = " $ sqrt rs') < tol = (x',i)          
          | otherwise                               = loop x' r' p' rs' (i+1)
          where mAp = mulMV mA p
                alpha = rs / (dot p mAp)
                x' = x + (alpha .* p)
                r' = r - (alpha .* mAp)
                rs' = dot r' r'
                beta = rs' / rs
                p'  = r' + (beta .* p)



(.*) :: (Num a) => a -> SpVCR a -> SpVCR a
(.*) s v = fmap (s *) v

编辑：果然，我没有考虑到 MM 文件格式只包含对称矩阵的下对角线这一事实。谢谢。现在算法收敛了，但似乎比它应该进行的迭代次数更多。我的理解是，当使用精确算术时，CG 应该始终以少于矩阵顺序的迭代次数收敛。使用浮点 (Double) 的事实是否会产生如此大的差异（1.5 - 2 x 矩阵阶数是合理收敛所需的迭代）？

跟进：对于任何可能偶然发现这一点的人，事实证明我的大部分问题都与我用于测试的矩阵有关。似乎它们对于使用 CG 算法进行求解的条件很差。在某些情况下，简单的预处理会有所帮助。

4

1 回答 1

1

您可以从这里使用带有浮动的精确库（例如 CReal）来回答您的第二个问题：http://hackage.haskell.org/package/numbers或摆脱您的日志记录（我认为这是引入浮动约束的原因）和仅使用 Data.Ratio 中的有理数。

这当然会非常缓慢。但它应该让您研究浮点近似对收敛的影响。

于 2013-02-03T20:06:06.527 回答