2

Haskell で行列を解くための Gauss Seidel と Conjugate Gradient 反復アルゴリズムの両方を作成しました (ただし、この質問は言語ではなくメソッドに関連しています)。私の理解では、これらのアルゴリズムは両方とも同様の収束特性を持つべきであり、ほとんどの場合 CG 法の方が高速であるはずです。http://math.nist.gov/MatrixMarket/の対称正定行列で多くのテストを実行しました。CG algを取得することはほとんどありません。GS はほぼ常に収束しますが、収束します。オンラインでテストする目的で右辺ベクトルが付随する対称正定行列が見つからないため、独自の RHS を任意に作成しています (これが問題の一部である可能性があります)。Ax = b で A の代わりに (transpose 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などのフローティングを備えた正確なライブラリを使用することで、2番目の質問に答えることができます: http://hackage.haskell.org/package/numbersまたはログを取り除く (これがフローティング制約を導入するものだと思います) Data.Ratio の合理性を使用するだけです。

もちろん、これは非常に遅くなります。しかし、収束に対する浮動小数点近似の影響を調査できるはずです。

于 2013-02-03T20:06:06.527 に答える