• ベストアンサー

共役勾配法 (Conjugate Gradient)の失敗?

共役勾配法で線形方程式を解く時に、誤差を10の-8乗 として計算をしています。 たいていの計算はこの誤差で解けるのですが、ある計算 では、誤差が10の-3乗くらいまで減った後に、誤差が増加 していってしまいます。 共役勾配法を使っていると、こういうことが起きるのでしょうか? 数値計算プログラムの作成ミスだろうかと考えてチェックしていますが、 今の所おかしな部分は見つかっておりません。

noname#56787
noname#56787

質問者が選んだベストアンサー

  • ベストアンサー
noname#221368
noname#221368
回答No.1

　共役勾配法は、桁落ちや丸め誤差の影響を受けやすいので、「前処理付き」で行えというのが標準だったと思いますが、「前処理」が面倒で、私はじつはやった事がありません。なので、最初に考えられる事は、 　　(1)前処理をなさっているか？ です。次に考えられる事は、 　　(2)係数行列が特異でないか？ です。前処理していない経験で言うと、係数行列が特異に近くなっても、共役勾配法はかなりもつ、という印象を受けました。しかし本当の特異行列を扱うなら、何らかの安全装置は必要かもしれません。 　行列が特異なら、コンピュータの内部精度ぎりぎりの小さな数での除算が、どこかに入るはずなので、そこで制御し、計算を打ち切るとか・・・。

noname#56787
質問者

お礼

詳細な回答ありがとうございました。 どうもワークステーションの種類かコンパイラの種類によってか、 上記の桁落ちか丸め誤差がでているようです。 前処理に関して調べてみます。 ありがとうございました。

  1. カテゴリ
  2. 学問・教育
  3. 数学・算数

関連するQ&A

  • 共役勾配法を用いた最適化について

    弾性体の変形を見るために，停留ポテンシャル法を用いようとしているのですが，その際に目的関数にひずみエネルギーをとり 独立な24変数の関数の極小値を共役勾配法で求めようとしています． Hesse行列を手計算でだし，最適化を行うプログラムを書いたのですが どうも上手くいきません． 下にプログラムの一部分をのせます． 共役勾配法について勘違いしている可能性もあるので 問題点を教えて頂ければ幸いです． x_vectorは24変数ベクトルで，dは探索方向，HはHesse行列です． for(i=0;i<10;i++) { lambda = (x_vector*d)/(d*H*d); x_save_vector = x_vector; x_vector=x_vector - lambda*H*d;//変数べクトル再設定 if(i==0) { // 初期値 gamma = 1; } else { gamma =x_vector*x_vector/(x_save_vector*x_save_vector); } d = x_vector + gamma*d;//探索方向再設定 }

  • 不完全ＬＵ分解前処理つき双共役勾配法についておしえてください。

    連立方程式を解くために不完全ＬＵ分解前処理つき双共役勾配法 について勉強しています。 前処理の際に、行列Ａを不完全ＬＵ分解しその逆行列(LU)^(-1)というのを使用します。LU分解まではできたのですが、この逆行列は普通にＬＵ分解＋直接法という形でもとめるのでしょうか。だとしたら、直接法をつかっていてあまり高速化が期待できない様な気がしました。 不完全コレスキー分解つき共役勾配法(ICCG)のときは、不完全コレスキー分解後、間接的にAの逆行列をもとめて使用する方法がありましたのでなにかいい方法があるのかと思い質問しました。 はじめてのプログラミングで見当違いなことをいっているかもしれませんがよろしくおねがいします。

  • ラプラス方程式の解法（情報処理として）と共役勾配法

    情報処理の数学はカテゴリ違いかもしれませんが、適当なところを見出せませんでしたので。 （以下、長文で申し訳ありません。） ラプラス方程式やポアソン方程式の境界値問題を解く方法にＳＯＲ法が あります。これは連立一次方程式の反復解法の１つだとも言えると思い ます。連立一次方程式の解法は通常AX=BとなってマトリックスAと既知 ベクトルBに対して未知ベクトルＸを解くわけですが、ラプラス方程式 をＳＯＲ法で解く場合、AX=Bという形にはしません（そういう風に表現 はできますが）。例えば2次元を考えると100×100の2次元の領域を解く 場合、未知数は10000個です。これをAX=Bという風に考えたら、10000元 連立一次方程式を解くことになり、A(10000,10000)のようにメモリ確保 も大変ですが、実際にはそのように解きませんね。マトリックスＡとい うメモリの取り方をしないのです。 100×100の格子領域でSOR法で解くとプログラムステップとしては数十 行ぐらいとなり、大したことはありませんね。 しかし、情報処理の連立一次方程式の解法の本を見ると、それはあくま でも連立１次方程式があった場合ということを前提としているのでその ようなメモリの確保から解説がスタートするようです。 そこで質問ですが、ラプラスやポアソン方程式を（前処理付）共役勾配 法などで解く場合、メモリ確保の前提条件としては連立一次方程式の形 ではないはずです。本を見るとマトリックス解法としての共役勾配法が 出ていますが、ラプラス方程式、ポアソン方程式の特化して解説してい る本などないでしょうか。また、簡単に説明できるのであれば教えて頂 きたいのですが。

  • 「セカント法」（割線法）についてです.

    非線形方程式の解を求めるアルゴリズム「セカント法」（割線法）についてです。 yahoo知恵袋の方でも質問しているのですが、 非線形方程式の解を求めるアルゴリズム「セカント法」（割線法）についてです。 セカント法で解を求める際、その収束次数は黄金比{（1+√5）/2}になると言われていますが、それを確認する方法はないものでしょうか。 セカント法で求めるプログラムはjavaの方で既に組んでおります。 なにかいいアイデアがあれば是非回答よろしくお願いします。 補足 組んだプログラムは与えられた関数の解と反復回数、それぞれの計算における関数値、真の解との誤差をかえすものです。 お願いします。

  • 最小二乗法における誤差の求め方

    こんばんは。 皆様よろしくお願いいたします。 あるデータにフィットさせる関数の係数を 最小二乗法を用いて、自分でプログラムを作って 見つけようと考えているのですが、 係数の誤差をどのように求めればいいかが わかりません。方法を探しても、直線の 場合はきれいに連立方程式が解けて、 誤差の伝播から計算するべき式が求まりますが、 一般の曲線の場合は解けないと思います。 gnuplotなどでフィッティングするとerrorが 出てきますが、あれはどのようにして 計算しているのでしょうか。 よろしくお願いいたします。

  • オイラー法、２次ルンゲクッタ法、４次ルンゲクッタ法のＣ言語プログラムに

    オイラー法、２次ルンゲクッタ法、４次ルンゲクッタ法のＣ言語プログラムについて教えてください！ 課題なのですが、まったくわからず困ってます＞＜ 1　常微分方程式　dy/dx=f(x,y),y(0)=1　の数値解をオイラー法を用いて計算するプログラムを作為せよ。ただし、f(x,y)=3-6x^2-4x+2xyとする。 2　α=1,β=1,γ=1/2,σ=1/2　の場合の２次ルンゲクッタ法を考える。1と同じ常微分方程式(f(x,y)も同じ)を考え、その数値解を求めるプログラムを作成せよ。また、オイラー法と２次ルンゲクッタ法の実行結果を示して、2つの近似精度を比較せよ。 3　1と同じ常微分方程式(f(x,y)も同じ)を考え、その数値解を４次ルンゲクッタ法を使って求めるプログラムを作成せよ。また、オイラー法、２次ルンゲクッタ法、４次ルンゲクッタ法の実行結果を示して、3つの近似精度を比較せよ。 以上の３つです。 休んでいた自分が悪いのですが、ネットで調べてもよくわからなくて… わかる方、よろしくおねがいします…

  • 非線形計画法について

    非線形計画法を現在勉強しています。 １． どういうときに線形でどういうときに非線形となるのか良く分かりません。 例えば、ある従属変数ｙを線形関数ｆ＝Σcx　で表したいときにパラメータｃの絶対値の和が定数ｂより小さくなるという制約のもとで、ｙとｆの二乗誤差を最小化するパラメータｃを求める問題を考えます。 この場合、制約条件はｃについて線形ですが、最小化したいのは、ｙとｆの二乗誤差なのでこの場合は非線形ということになるのでしょうか？それとも関数ｆはｃに関して線形関数なので、線形計画法で解くことになるのでしょうか？ ２． 以下のサイトで勉強しているのですが、このサイトにある楕円型の等高線はおそらく、従属変数ｙと目的関数ｆの誤差を表しているのだと思うのですが、なぜ「楕円」になるのですか？二乗誤差を考えるのならば、「円」になるのではと疑問で仕方ありません。 http://www.sist.ac.jp/~suganuma/kougi/other_lecture/SE/opt/nonlinear/nonlinear.htm#2.2 疑問が晴れずにもやもやしています。 回答もしくはアドバイス、よろしくお願いします。

  • ニュートン法に関して

    数値計算初心者です。数値計算で分からないことがあるので質問します。よろしくお願いします。 y=f(x,a)という関数があってパラメータaを非線形最小２乗法のニュートン法やマルカート法を使って求めたいのですが、計算過程でf(x)を各パラメータで偏微分してヤコビ行列を求める必要があると思われます。 例えばf(x)が複雑な関数で偏微分するのに困難な関数であった場合、 偏微分をしなくてΔxを決定するにはどのような方法があるのでしょうか？

  • 連立1次方程式の数値解法の使用条件

    数値計算で物理現象を解いていく問題では多くの場合、最終的には大規模な連立1次方程式を解くというところに行きつきます（行列を解く）。その場合、解く物理問題によって行列の形式が決まってくる場合とか、あるいはその性質は期待できない問題とがあると思います。例えば、行列が必ず（正値の？）対称行列になることがわかっている問題とか、です。では行列の性質を期待しない場合の高速解法にはどのようなものがあるでしょうか。ガウスの消去法はまるで紙と鉛筆で解いていくことをプログラミングしたようなところがあります。ピボットの問題とかはありますが。でもたぶん遅いんだろうなあとは思いますが。 教科書を読んでみるとだいたい、細かい解説が書いてありますが、ユーザという立場からは内容はどうでもいいから、自分の問題に対して早くて正確な解を得る方法だけ教えて欲しいと思うのですが。共役勾配法を使った計算例があったので、勉強してプログラムを作って走らせてみたら全く結果がおかしいのですが、よく読んでみたら、共役勾配法は正値対象行列に限定だそうで、がっかりしました。 行列の性質を限定しない高速解法にはどのようなものがあるのでしょうか。なお、よく問題になる条件数の問題とか10^10と10^(-10)が係数に含まれるとかいろいろ問題がありますが、今回の問題としてはそのような極端なことが起こらないということではあります。共役勾配法はいろんなファミリーがありますが、その中で正値対象行列でなくてもいいというものあるでしょうか。解説書には解法の冒頭に書いてもらいたいものですが。 よろしくお願いします。

  • ニュートン法

    こんにちは。 aを解とするある非線形方程式に,aの近傍の値x0を初期値としてニュートン法を一反復だけ適用したところ、誤差は |x0-a|=10^-4から|x1-a|＝10^-7に減少した。このとき、ニュートン法をもう一反復適用したら誤差|x2-a|はどの程度になると推定されるのでしょうか？

注目のQ&A

カテゴリ

専門家に質問してみよう

職業から探して質問する

あなたにピッタリな商品が見つかる！ OKWAVEセレクト

質問する