はじめての親子ハイキングに挑戦!! >>

よろしくおねがいします。

タイトルの通り、
n✕n行列の固有値を求めるプログラムを作成しようと考えています。
ただし、行列は非対称行列とするためJacobi法等は使えません。

そこで、
このプログラムを作成する際の一般的なアルゴリズムを教えていただきたいです。
例えば、どういった法則を使うのか?などです。
具体的であればあるほどありがたいです。

しょぼい質問ですがお願いします。

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

A 回答 (5件)

4 x 4 でしたら、固有方程式を直接といたほうが簡単です。


ヤコビとか QR とかは巨大な行列に使う手法です。

ベアストウヒッチコックとかDKAとか、いろいろとやり方があります。
    • good
    • 0

肝心なところを書き忘れました。



原理的には QR だけで固有値は求まります。小さな行列ではこれで十分ですが、
数十、数百の大きさの行列では QR だけでは収束が遅すぎるので、
他の手法(加速処置)との組み合わせが必要です。とんでもなく速くなります。
    • good
    • 0
この回答へのお礼

具体的なところまでありがとうございます。

実際求めたい行列は4✕4の非対称行列なので
QR法のみで大丈夫な気がします。

QR法は初めて聞くアルゴリズムなので
勉強してから実装したいと思います!

実際のところ
Cで組むとしたら大変でしょうか??

お礼日時:2012/10/06 15:24

私が勉強した頃は、非対称行列では



Reduction(ヘッセンベルク化、3重対角化) + QR + 原点移動 + 減次

が定石だったと思います。最近は違っていたらすいません。

Reduction(ヘッセンベルク化、3重対角化) はハウスホルダー法を実装したことが
ありますが、他にもいろいろなやり方があるようです。
#双ランチョス法とか・・・実装経験なし。対称行列のように3重対角化できるようなので
#メチャ速いかも。

Reduction で計算量を減らしたあと、QR + 原点移動 + 減次 をセットで行うと
収束が速いようです。

逆にセットでやらないととても遅くなって実用的ではありませんでした。

拙い情報ですが・・・
    • good
    • 0

条件にもよるけど, べき乗法 (逆べき乗法) や LR法など, いくつかあるよ.



というか, 調査した?

この回答への補足

回答ありがとうございます。

条件は非対称行列である、
すべての固有値を求めたい、
C言語を使いたいぐらいです。

調べて見たところ、
フレーム法+ベアストウ法を用いるものは見つけましたが、
なかなか非対称行列の例題がないため他はわかりませんでした。

そこで
その他にもアルゴリズムがあるのかを知りたくて投稿しました。

言葉足らずですみません。

補足日時:2012/10/05 18:26
    • good
    • 0

企業で統計を扱う部門の者です。



固有値は対称行列しか求めることができません。

非対称行列の固有値を求めたいということは、
ランク落ちに近い、すなわち最小固有値がきわめて0に近い
状態を知りたいのでしょうか。

もし、そうであれば、

rank(A)=rank(A'A)=rank(AA')

という公式がありますので、
それを用いて対称行列にしてから固有値を求めてはどうでしょうか。
    • good
    • 0
この回答へのお礼

回答ありがとうございます。

簡単にいうと近似して求めるということですね。

参考にさせていただきます。

お礼日時:2012/10/05 19:01

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

計算 逆行列」に関するQ&A: 平面の計算方法

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

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

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

Q非対称行列の固有値と正定値性について

画像にあるように、非対称な行列について質問があります。

ある英語の本の問題で、画像に書いてある行列について以下のような問があります。
Is the matrix positive definite?
(これは正定値行列ですか?)

正定値行列がわからなかったので、調べたのですが、普通は対称行列に対して求めるもので、非対称な行列に対してそもそも取り扱われていないようでした。

また、正定値行列は固有値が全て正であるとのことだったのですが、この行列の固有値を求めたところ、複素数が出てきました。

これって正定値行列かどうか、判定することができるのでしょうか?

Aベストアンサー

No.3のエルミート行列の下りは間違い

nを2以上の整数とし
Aをn次エルミート行列とし
xをn次複素列ベクトルとする
(任意の行列XについてX^*をXの複素共役転置とする)
零ベクトルでない任意のxについて
S=x^*・A・x
が正であるときAを正値であるという
Aが正値⇄Aの実数部行列が正値
は成り立たない
従ってエルミート行列の正値性は対称行列の正値性の拡張として意味がある

従って修正版としては後半を除いて以下の通り

nを2以上の整数とし
Aをn次実行列とし
xをn次実列ベクトルとする
(任意の行列XについてX^TをXの転置とする)
零ベクトルでない任意のxについて
S=x^T・A・x
が正であるときAを正値であるという
S=S^Tなので
S=(S+S^T)/2=x^T・(A+A^T)・x/2
であるから
Aが正値⇄(A+A^T)/2が正値
M=(A+A^T)/2
は実対称行列であるから
行列が正値であるかどうかを判定するには
その行列を対称化した行列により判定すればよい
(行列Aの対称化行列とは(A+A^T)/2のことである)
従って実非対称正方行列の正値性は
実対称行列の正値性の拡張としてあまり役に立たない

今回の場合
√2・M=√2・(A+A^T)/2=
[1 0]
[0 1]
でありMの固有値は1/√2であり正であるから
与行列は正値である

No.3のエルミート行列の下りは間違い

nを2以上の整数とし
Aをn次エルミート行列とし
xをn次複素列ベクトルとする
(任意の行列XについてX^*をXの複素共役転置とする)
零ベクトルでない任意のxについて
S=x^*・A・x
が正であるときAを正値であるという
Aが正値⇄Aの実数部行列が正値
は成り立たない
従ってエルミート行列の正値性は対称行列の正値性の拡張として意味がある

従って修正版としては後半を除いて以下の通り

nを2以上の整数とし
Aをn次実行列とし
xをn次実列ベクトルとする
(任意の行列XについてX^TをX...続きを読む

QC言語のプログラムで質問です。

C言語のプログラムで質問です。
ヤコビ法で固有値・固有ベクトル絵を求めたいのですが、固有値は上手く出るのですが、
固有ベクトルがめちゃくちゃな値が出ます。何が間違っているのでしょうか教えてください。下のは一応自分で作ってみたプログラムです。

#include <stdio.h>
#include <math.h>
#define N 4
#define PI 3.14159
void jacobi(double a[N][N], double λ[N], double v[N][N]);
int main(void)
{
int i, j;
double a[N][N],
λ[N], v[N][N], x;

for(i=0;i<N;i++){
printf("\n a[%d][0] a[%d][1] a[%d][2] a[%d][3]=",i,i,i,i);
scanf("%lf %lf %lf %lf",&a[i][0],&a[i][1],&a[i][2],&a[i][3]);
}

jacobi(a, λ, v);

for(j=0; j<N; j++){
printf("\n固有値λ[%d]=%lf", j, λ[j]);
printf("\n固有ベクトル\n");
for(i=0; i<N; i++){
printf("%lf\n", v[i][j]);
}
}
return 0;
}
/* ヤコビ法による固有値計算 */
void jacobi(double a[N][N], double λ[N], double v[N][N])
{
int i, j, kmax=100, repeat, p, q;
double eps, c, s,theta, gmax;
double apq, app, aqq, apqmax, apj, aqj, vip, viq;

gmax=0.0;
for(i=0; i<N; i++){
s=0.0;
for(j=i+1; j<N; j++){
s += fabs(a[i][j]);
}
if(s>gmax) gmax=s;
}
eps=0.000001*gmax;
for(i=0; i<N; i++){
for(j=0; j<N; j++){
v[i][j]=0.0;
}
v[i][i]=1.0;
}

for(repeat=1; repeat<kmax; repeat++){
/* 収束判定 */
apqmax = 0.0;
for(p=0; p<N; p++){
for(q=0; q<N; q++){
if(p!=q){
apq=fabs(a[p][q]);
if(apq>apqmax) apqmax=apq;
}
}
}
if(apqmax<eps) break;
for(p=0; p<N-1; p++){
for(q=p+1; q<N; q++){
apq=a[p][q];
app=a[p][p];
aqq=a[q][q];
if(fabs(apqmax)<eps) break;
/*回転角計算*/
if(fabs(app-aqq)>=1.0e-15){
theta= 0.5*atan(2.0*apq/(app-aqq));
}else{
theta= PI/4.0;
}
c = cos(theta);
s = sin(theta);
a[p][p] = app*c*c + 2.0*apq*c*s + aqq*s*s;
a[q][q] = app*s*s - 2.0*apq*c*s + aqq*c*c;

a[p][q] = 0.0;
a[q][p] = 0.0;

for(j=0; j<N; j++){
if(j!=p && j!=q){
apj = a[p][j];
aqj = a[q][j];
a[p][j] = apj*c + aqj*s;
a[q][j] = -apj*s + aqj*c;
a[j][p] = a[p][j];
a[j][q] = a[q][j];
}
}
/*固有ベクトル*/
for(i=0; i<N; i++){
v[i][p] = v[i][p]*c+v[i][q]*s;
v[i][q] = -v[i][p]*s+v[i][q]*c;

}
}
}
eps=eps*1.05;
}
for(i=0; i<N; i++) λ[i] = a[i][i];
}

C言語のプログラムで質問です。
ヤコビ法で固有値・固有ベクトル絵を求めたいのですが、固有値は上手く出るのですが、
固有ベクトルがめちゃくちゃな値が出ます。何が間違っているのでしょうか教えてください。下のは一応自分で作ってみたプログラムです。

#include <stdio.h>
#include <math.h>
#define N 4
#define PI 3.14159
void jacobi(double a[N][N], double λ[N], double v[N][N]);
int main(void)
{
int i, j;
double a[N][N],
λ[N], v[N][N], x;

for(i=0;i<N;i++){
printf("\n a[%d][0] a[%d][1]...続きを読む

Aベストアンサー

> 何が間違っているのでしょうか

/*固有ベクトル*/
for(i=0; i<N; i++){
v[i][p] = v[i][p]*c+v[i][q]*s;
v[i][q] = -v[i][p]*s+v[i][q]*c;
}

forの中の1行目で値を変えた後のv[i][p]を2行目の右辺で使っている

Q固有値が複素数のときの固有ベクトルの求め方

固有値が複素数のときの固有ベクトルの求め方

( -7 -5 )
( 13 9 )

の2x2行列で固有値を求めると 1±2i になると思いますが

Av = λv の形で固有ベクトルを求めようとすると

( -8 + 2i ) x - 5 y = 0
13 x + ( 8 + 2i ) y = 0

の形になり、その先を求めることが出来ません。
何度も計算したので最後の2つの式は間違いは無いと思うのですが、
固有値が複素数の時は、Av = λv の方法で計算することは出来ないということでしょうか?
またどのように計算できるのでしょうか?
お知恵をお貸しいただければ幸いです。

Aベストアンサー

固有値は1±iになるかと…

そこから先の計算は普通に実数の時と同じ方法で計算できます.

Qヤコビ法とQR法について

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

Aベストアンサー

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度計算し、両者の結果が一致することを確認すべきです。もし食い違うようなら、有効桁数が足りない訳です。

n×n行列Aの固有値問題。
困り度3ということですから、中途半端ながらres付けてしまいます。

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

反復計算で固有値と固有ベクトルを求める段階では、固有値の絶対値同士がどのぐらい散...続きを読む

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ランキング