共轭梯度法收敛

我已经编写了Gauss Seidel和Conjugate梯度迭代算法来解决Haskell中的基元问题（但这个问题与方法不是很相关）。我的理解是，这两种算法应该具有相似的收敛特性，并且CG方法在大多数情况下应该更快。我已经在http://math.nist.gov/MatrixMarket/的对称正定矩阵上运行过很多测试，而且我几乎可以不用CG阿尔格。融合，而GS几乎总是这样。我无法找到任何对称的正定矩阵和相应的右手边矢量用于在线测试目的，所以我一直在任意创建我自己的RHS（也许这是问题的一部分？）。如果我在Ax = b中使用（转置A）* A而不是A，则可以使CG方法收敛，这只是迫使矩阵对称。我在这里包含了CG代码。它显然不会按原样编译。如果有人需要它的功能来帮助，我会发布一切。它正在为来自(Pseudocode and example)的简单示例(Similar question)正确工作。关于共轭梯度与高斯赛德尔会聚标准有什么不足之处？任何人都可以指出我正确的方向来实现这个目标吗？谢谢。共轭梯度法收敛

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算法时病态很大。简单的预处理在某些情况下有所帮助。

来源

2012-12-12 MFlamer

“如果我只是强制矩阵是对称的，我可以使CG方法收敛。”不是'A'已经对称了吗？ –

我是这么认为的，那是我混乱的一部分。这是矩阵http://math.nist.gov/MatrixMarket/data/Harwell-Boeing/bcsstruc2/bcsstk14.html。我只注意到它不是对角线的主宰。我认为这是高斯赛德尔集中要求，但不是CG。也许这就是问题所在。如果是这样，CG似乎只适用于特定类别的矩阵。 – MFlamer

CG需要对称正定（spd）矩阵进行收敛。 Gauss-Seidel仅保证对角占优和spd矩阵收敛。 –

您可以通过使用浮动的确切库（如来自此处的CReal）来回答第二个问题：http://hackage.haskell.org/package/numbers或删除您的日志记录（我认为引入了浮动约束），并仅使用Data.Ratio 。

这当然会非常慢。但它应该让你调查浮点近似对收敛的影响。

来源

2013-02-03 20:06:06 sclv

谢谢，我没有想过这样做。 – MFlamer

共轭梯度法收敛

回答

相关问题