現在固有値問題を解いているのですが、ヤコビ法とQR法とでは
どちらがどれくらい計算が速いでしょうか?
解こうとしている行列は最小で1000×1000で、精度を良くする
ために5000×5000程度の行列で解こうと思っています。
どなたか教えていただければ幸いです。お願いします。

このQ&Aに関連する最新のQ&A

A 回答 (2件)

n×n行列Aの固有値問題。


困り度3ということですから、中途半端ながらres付けてしまいます。

ヤコビ法は大きいnには手間が大きくて辛いでしょうね。
QR法が便利だけれど、一手間掛けまして、原点を移動する。つまり狙いの固有値を絶対値≒0になるように座標変換しておいてからQR法を使うのが旨いやり方だと思います。つまり予想される値μを使って、Aの代わりにA-μI の固有値を求めて、得られた答えにμを加えればよい。

反復計算で固有値と固有ベクトルを求める段階では、固有値の絶対値同士がどのぐらい散らばっているかによって計算の手間が全然違う。絶対値の順に固有値を並べた時、隣り合う固有値の比が1に近いと手間が掛かります。固有値が一個ずつ出てくるのを使って問題を減次するのも、nが大きいときには結構な手間です。なかなか一概には論じきれないな。

また固有値と固有ベクトルの近似値μ,vを改良するという仕上げは欠かせません。つまり
(A-μI)x = v
を解いてxを求め、改良された固有ベクトルの近似値xと、固有値の近似値(Ax)[j]/x[j]を求める。これを繰り返します。(これは最初の近似が良ければ凄く速く収束します)。QR法はこの計算にも向いてます。

スパースな行列(つまり零要素が殆ど、という場合)だと、また話が違うようです。

いずれにせよ、誤差の累積が怖いので、計算の精度(有効桁数)を変えて2度計算し、両者の結果が一致することを確認すべきです。もし食い違うようなら、有効桁数が足りない訳です。
    • good
    • 0
この回答へのお礼

とても有用なご回答ありがとうございます。
もう少しいろいろ調べながらやっていこうと思います。
ありがとうございました。

お礼日時:2001/12/27 23:17

大きい固有値だけに限って求められれば十分、という応用が沢山あります。

また、制御工学では不安定な固有値成分だけ取り出したい。そういう場合どうするのか、ちょっと調べてみたら
「ランチョス法に、等角写像を組み合わせる。」
というテクニックがあるみたいです。
等角写像
y = (x+c) / (x-c) (x,yは複素数、cは適当な正の定数)
によって、不安定なやつを単位円の外に、安定な奴を単位円の中に写像しておいて、ランチョス法(絶対値が大きい固有値から順に出てきます)を適用するらしい。
 まだ調査中ですが、お急ぎのようなのでとりあえず。
    • good
    • 0

このQ&Aに関連する人気のQ&A

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人はこんなQ&Aも見ています

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Q固有値問題

現在、因子分析や主成分分析などを行うために対称行列の固有値及び固有ベクトルを求めるプログラムを作っています。
現在、ハウスホルダー法を用いて行列を3重対角化し、QR法にて固有値が求まるところまでは、できたのですが、この後の固有ベクトルを求める方法がわかりません。
文献などを調べておりますと、逆反復法を用いて固有ベクトルを求めれば良いとあるのですが、この逆反復法のアルゴリズムが良く分かりません。
また、この逆反復法以外で効率的に求める手法がありましたら、教えてください。
よろしくお願いします。

Aベストアンサー

逆反復法は固有ベクトルを計算するのによく用いられます.

問題の行列を A (既知),近似的固有値を Ea (既知) とします.
任意の初期ベクトル v_0 を出発点にとり,
(A-Ea)^(-1) を繰り返し掛けようというのが逆反復法のアルゴリズム
(というか,原理)です.
v_0 を A の固有ベクトル φ_j (未知)で展開します.
(1)  v_0 = Σ_j a_(j0) φ_j
係数 a_(j0) は未知.
これに (A-Ea)^(-1) を k 回掛けます.
左辺は (A-Ea)^(-k) v_0 = v_k と書くことにします.
φ_j は A の固有ベクトル(固有値を E_j とします)なので,
(2)  (A-Ea)^(-1) φ_j = (E_j - Ea)^(-1) φ_j
で,k 回繰り返して
(3)  (A-Ea)^(-k) φ_j = (E_j - Ea)^(-k) φ_j
したがって
(4)  v_k = Σ_j a_(j0) (E_j - Ea)^(-k) φ_j
です.
(E_j - Ea)^(-k) がミソで,Ea に近い E_j をもつ j 成分が他の成分に比べて
格段に大きくなります(繰り返し回数 k に対して指数関数的に増大).

(A-Ea)^(-1) を掛ける操作をそのままやろうとすると逆行列を求める必要があります.
しかし,
(5)  v_k = (A-Ea)^(-1) v_(k-1)
ですから
(6)  v_(k-1) = (A-Ea) v_k
というわけで,これなら単に連立方程式を解くだけですからできます.

不幸にして最初の v_0 の選び方が悪くて,求める固有ベクトルと直交していると,
うまく行きません.
極端な話,v_0 を A の真の固有ベクトルに選んでしまったら,
いつまで経ってもその固有ベクトルから抜け出せません.
こういうことが起きるのは,対称性の問題に起因することがほとんどです.
まあ,適当にごちゃごちゃした A で,v_0 を対称性の違うベクトルの線形結合に
とっておけばたいてい大丈夫です.

実際のプログラミングにはいろいろ注意が要ります.
固有値が縮退している可能性なども考えておかないといけません.
森正武「FORTRAN77 数値計算プログラミング」(岩波書店,1988)
など参考にされてはいかがでしょうか.

逆反復法は固有ベクトルを計算するのによく用いられます.

問題の行列を A (既知),近似的固有値を Ea (既知) とします.
任意の初期ベクトル v_0 を出発点にとり,
(A-Ea)^(-1) を繰り返し掛けようというのが逆反復法のアルゴリズム
(というか,原理)です.
v_0 を A の固有ベクトル φ_j (未知)で展開します.
(1)  v_0 = Σ_j a_(j0) φ_j
係数 a_(j0) は未知.
これに (A-Ea)^(-1) を k 回掛けます.
左辺は (A-Ea)^(-k) v_0 = v_k と書くことにします.
φ_j は A の固有ベクトル(固有値を E_j と...続きを読む

Q行列演算: 固有ベクトルの解法

現在、対称行列の固有値、固有ベクトルを求めるプログラムを作成し、つい最近完成しました。

しかし、とても使い物にならないプログラムになってしまいました。
理由はとても遅いのです。

解法の手順として、まず固有値を求めてから固有ベクトルを求めるようと考え、入力の対称行列をHouseHolder法により三重対角行列に変換し、それをQR法により対角化してまず固有値を求めました。

固有値を求めることができたので、次に固有ベクトルを求めます。手順として、固有値ごとに入力対称行列の対角成分から固有値を減算した行列をLU分解し、連立一次方程式を解くように固有ベクトルを求めていきます。

この一連の手順で、対称行列の固有値、固有ベクトルを求めることができたのですが、とても時間がかかってしまいます。

ただし、対称行列の固有値を求めるまでの時間はとても高速です。

500×500の行列の固有値、固有ベクトルを求めるのに30分はかかってしまいますが、その中で固有値を求める時間は2秒しかかかりません。

つまり今固有値がわかっている状態で、固有ベクトルを高速に求めたいと考えています。

なにか高速に固有ベクトルを求める方法(アルゴリズム)はあるでしょうか?

現在、対称行列の固有値、固有ベクトルを求めるプログラムを作成し、つい最近完成しました。

しかし、とても使い物にならないプログラムになってしまいました。
理由はとても遅いのです。

解法の手順として、まず固有値を求めてから固有ベクトルを求めるようと考え、入力の対称行列をHouseHolder法により三重対角行列に変換し、それをQR法により対角化してまず固有値を求めました。

固有値を求めることができたので、次に固有ベクトルを求めます。手順として、固有値ごとに入力対称行列の対角成分から固...続きを読む

Aベストアンサー

 一つ気づいたのですが、もしかして500個の固有値を全部求めていませんか?。そうすると全体で30分(1800秒)だとしても、1個あたりなら、ほんの数秒で、かなり優秀なコードと思えます。
 ふつうは、低次の固有値・固有ベクトルを重視するので、せいぜい10個,多くて20個と思います。さっきの話が正しいとして、それなら数10秒です。
 もし500個の固有値という事であれば、共役勾配法を利用した原点移動付逆ベキ乗法を試すしかないかな?、と思いました。

Q固有値、固有ベクトル、対角化...何のため?

私は文系出身の32歳会社員です。

ふとしたきっかけで数学を学び直そうかなと
独学で最近始めました。

そこで...
本当に素朴で基本的な疑問で恐縮なのですが...

(1)何のために固有値を求めるのでしょうか?
(2)何のために固有ベクトルを求めるのでしょうか?
(3)何のために行列の対角化を行うのでしょうか?

回答は歴史的背景、学術的背景、感情...etc、なんでも結構です。
例)
・特定の法則で計算すると固有値が求められるので求めた。
・固有ベクトルは縦に並べてベクトルとしてみた方がすっきりするから「数列」ではなく「ベクトル」と呼んでみた。
・意味はない!目的はない!ただ数学として突き詰めているだけだ!
...などなど

あっ、でも急を要している訳ではないので
もしご存知の方、もしくは自論をお持ちの方は
お時間のある方はご回答いただければ幸いです。

ちなみにテキストは共立出版の『やさしく学べる基礎数学~線形代数・微分積分~です。
やっと線形代数が終わって、微分積分に入ろうというところで、ふと疑問を持ってしまいました...(~~;

本当に漠然とした質問で恐縮ですが、どうぞよろしくお願いいたしますm(__)m

私は文系出身の32歳会社員です。

ふとしたきっかけで数学を学び直そうかなと
独学で最近始めました。

そこで...
本当に素朴で基本的な疑問で恐縮なのですが...

(1)何のために固有値を求めるのでしょうか?
(2)何のために固有ベクトルを求めるのでしょうか?
(3)何のために行列の対角化を行うのでしょうか?

回答は歴史的背景、学術的背景、感情...etc、なんでも結構です。
例)
・特定の法則で計算すると固有値が求められるので求めた。
・固有ベクトルは縦に並べてベクトルとしてみた方がすっ...続きを読む

Aベストアンサー

行列を 1個固定して考えてみます.
この行列は何かよくわからないんですが線形変換を表します.
たいていのベクトルはこの線形変換によって変な方向を向いてしまうんですが, まれに方向が変わらず長さだけが変わるベクトルがあります.
このように「長さだけが変わるベクトル」がこの線形変換 (ひいては行列) の固有ベクトルとなります. で, 長さの変化率が固有値.

Q行列の正定・半正定・負定

行列の正定・半正定・負定について自分なりに調べてみたのですが、
イマイチ良くわかりません。。。
どなたか上手く説明していただけないでしょうか?
過去の質問の回答に

>cを列ベクトル、Aを行列とする。
>(cの転置)Ac>0
>となればAは正定値といいます。
>Aの固有値が全て正であることとも同値です。

とあったのですが、このcの列ベクトルというのは
任意なのでしょうか?
また、半正定は固有値に+と-が交じっていて、
負定は固有値が-のみなのですか?

どなたかお願いしますorz

Aベストアンサー

まず、行列の正定・半正定・負定値性を考えるときは、
行列は対称行列であることを仮定しています。
なので、正確な定義は、

定義 n次正方 "対称" 行列 A が正定値行列であるとは、
『ゼロベクトルではない任意の』n次元(列)ベクトル c に対して、
(cの転置)Ac>0
となることである。

です。

対称行列Aが正定値なら、その固有値はすべて正です。
(cとして固有ベクトルをとってみればよいでしょう。)
逆に、対称行列Aの固有値がすべて正なら、Aは正定値行列です。

ただし、対称行列ではないAの固有値がすべて正だからといって、
(cの転置)Ac>0とは限りません。
例えば、
A =
[ 1 4 ]
[ 0 1 ]
とすると、Aは対称行列ではなく、固有値は1です。
しかし、
(cの転置) = [ 1, -2]
とすると、
(cの転置)Ac = -3 < 0
となってしまいます。(実際に計算して確かめてください。)
なので、行列Aが対称行列であるという条件はとても重要です。

また、半正定値の定義は、上の定義で
『ゼロベクトルではない任意の』 --> 『任意の』
と書き直したものです。
このとき、半正定値行列の固有値はすべて0以上です。(つまり0も許します。)
逆に、対称行列の固有値がすべて0以上なら、その行列は半正定値です。

さらに、負定値の定義は、『ゼロではない任意の』ベクトルcに対して
(cの転置)Ac<0
となることです。
固有値についてはもうわかりますね。

まず、行列の正定・半正定・負定値性を考えるときは、
行列は対称行列であることを仮定しています。
なので、正確な定義は、

定義 n次正方 "対称" 行列 A が正定値行列であるとは、
『ゼロベクトルではない任意の』n次元(列)ベクトル c に対して、
(cの転置)Ac>0
となることである。

です。

対称行列Aが正定値なら、その固有値はすべて正です。
(cとして固有ベクトルをとってみればよいでしょう。)
逆に、対称行列Aの固有値がすべて正なら、Aは正定値行列です。

ただし、対称行列...続きを読む

QLNK2019: 未解決の外部シンボルのエラーが出る

Microsoft Visual Studio 2008
Version 9.0.21022.8 RTM
Microsoft .NET Framework
Version 3.5 SP1
----------------------------------------------------------------
新しいプリジェクト→Win32 コンソール アプリケーション(ソリューションのディレクトリを作成 チェック外す)→Windows アプリケーション(空のプロジェクト チェック外す)
----------------------------------------------------------------
 プログラム

 mymain.cpp
#include "myhelper.h"
#include "mymain.h"

//自キャラのデータ
Point2D g_jikipos = {40, 400};//自キャラの座標

//画像ハンドル
int g_jikiimage[11];

//色々なファイルの読み込み
int LoadFiles(){
//画像ファイル読み込み
if(LoadDivGraph("media\\player01.bmp",
11,11,1,64,64,g_jikiimage) == -1) return -1;

return 1;
}


 mymain.h
//他から呼び出させるMyMainの関数
void MyMain();
int LoadFiles();


 myhelper.h(サンプルなので打ちミスはない)
#include "DxLib.h"
#include <limits.h>
#include <math.h>

//構造体宣言
//座標またはベクトルを記録する構造体
struct Vector{
float x,y;
};
typedef Vector Point2D;
//線を記録する構造体
struct Line2D{
Point2D startpos, endpos;
float katamuki;//傾きをラジアン値で記録
Vector speed;//移動している場合は速度をセット
};
//球体を記録する構造体
struct Ball2D{
Point2D position;
float hankei;//半径
};
//四角形を記録する構造体
struct Rect2D{
Point2D lefttop;
Point2D rightbottom;
float width;
float height;
};


//ライブラリ関数
Point2D PosInView(Point2D in);
int XInView(float inx);
int YInView(float iny);
void ScrollToLeft(float jikiposx);
void ScrollToRight(float jikiposx);
void ScrollToUp(float jikiposy);
void ScrollToDown(float jikiposy);
void DrawLineInView(float x1, float y1, float x2, float y2, int Color, int Thickness);
void DrawCircleInView(float x, float y, float r, int Color, int FillFlag);
void DrawAnimation(float x, float y, double ExtRate, double Angle,int TurnFlag,
int *imgarray, int allframe, float fps);
//ベクトル関数
Vector CreateVector(Vector in, float veclen);
Vector AddVector(Vector v1, Vector v2);
Vector SubVector(Vector v1, Vector v2);
Vector AddVectorInFrameTime(Vector pos, Vector speed);
Vector AddVectorInFrameTime2(Vector pos, Vector speed, Vector accel);
Vector Normalize(Vector in);
Vector RotateVector(Vector in, float radian);
float VectorLengthSquare(Vector in);
float DotProduct(Vector v1, Vector v2);
float CrossProduct(Vector v1, Vector v2);
void SetLine2DKatamuki(Line2D *in);
void DrawLine2D(Line2D in, int Color, int Thickness);
void DrawBall2D(Ball2D in, int Color, int Fill);
//当たり判定関数
bool HitTestLineAndBall(Line2D linein, Ball2D ballin);
bool IsPointAtLineFace(Line2D linein, Point2D ptin);
bool HitTestLineAndLine(Line2D line1, Line2D line2);
bool HitTestBallAndBall(Ball2D a, Ball2D b);
bool HitTestPointAndBox(Rect2D rect, Point2D pt);
//タイマー関数
void SetSimpleTimer(int idx, int time);
int GetPassedTime(int idx);


//グローバル変数
extern float g_frametime;
extern Rect2D g_framerect;//画面領域(当たり判定)
extern Point2D g_current_field_pos;//現在の左上座標
extern Rect2D g_stagesize;//ステージサイズ

//定数宣言
const float ZEROVALUE = 1e-10f;
const float PIE = 3.1415926f;
const int SCROLL_LIMIT = 200;
----------------------------------------------------------------
 エラー内容
1>myhelper.obj : error LNK2019: 未解決の外部シンボル "void __cdecl MyMain(void)" (?MyMain@@YAXXZ) が関数 _WinMain@16 で参照されました
1>C:\Documents and Settings\Owner\My Documents\Visual Studio 2008\Projects\my\Debug\my.exe : fatal error LNK1120: 外部参照 1 が未解決です
1>my - エラー 2、警告 0
ビルド: 0 正常終了、1 失敗、0 更新不要、0 スキップ
----------------------------------------------------------------
画像を貼り付けときます
(見えにくい場合→http://www.dotup.org/uploda/www.dotup.org154142.jpg.html)
初心者なのでわかりやすくお願いします

Microsoft Visual Studio 2008
Version 9.0.21022.8 RTM
Microsoft .NET Framework
Version 3.5 SP1
----------------------------------------------------------------
新しいプリジェクト→Win32 コンソール アプリケーション(ソリューションのディレクトリを作成 チェック外す)→Windows アプリケーション(空のプロジェクト チェック外す)
----------------------------------------------------------------
 プログラム

 mymain.cpp
#include "myhelper.h"
#include "mymain.h"

//自...続きを読む

Aベストアンサー

ファイル構成から推測するに
mymain.cpp というファイルに
void MyMain(void) {
// ここに処理を書く
}
という関数が必要なようです。

Qエルミート行列の固有値

エルミート行列の固有値は必ず実数になることはどうやって示せますか。

Aベストアンサー

どこの教科書にも載ってるような話だけど、
本で調べなかったの?

行列 A の転置共役を A* と書くことにする。
行列 H がエルミートであるとは、H* = H のこと。

H の固有値 λ に属する固有ベクトルを x と置く。
(x*)Hx = (x*)(Hx) = (x*)λx = λ(x*)x.
また、
(x*)Hx = (x*)(H*)x = ((Hx)*)x = ((λx)*)x = (λ*)(x*)x.
固有ベクトル x は零ベクトルではないから、
λ(x*)x = (λ*)(x*)x の両辺を (x*)x で割って、
λ = λ*. これは、λ の共役が λ と等しいこと、
つまり、λ が実数であることを示している。


人気Q&Aランキング