C言語で N行*M列 の逆行列を求める方法は

大体５０×１００００のような大きい配列の行列と
５０＊５０ほどの配列と、２つの配列の
逆行列を求めたいと思っています。

この場合どのような方法を使って求めればいいでしょうか？

C言語で求めたいのですが、
フリーのライブラリなどあるのならば、それでもいいですし、
可能なら、ここに関数を書いていただいてもいいです。

とにかく求められればいいので、アドバイス頂きたいです。
よろしくお願いします。

No.5 ベストアンサー 回答者： echoes_x86 回答日時： 2008/09/29 21:48

こんばんは．

特異値分解はm*n行列Xの転置であるn*m行列X^Tを分解しても一般性を失わないので，
とりあえず(m>n)と書いたのがまずかったようですね．
Xが横長でも縦長でも
X = U*S*V^T
と分解したならば，その擬似逆行列は
X^+ = V*S^(-1)*U^T
です．
(m>n)の場合は以下のようにすれば n*n単位行列Eが計算できます．
(X^+)*X = V*S^(-1)*U^T*U*S*V^T = E
(m<n)ならば次で m*m単位行列Eが計算できます．
X*(X^+) = U*S*V^T*V*S^(-1)*U^T = E
UとVの列が正規直交であることに注意して確認してください．

当初のご質問とずれてきていますので，サイトの主旨に照らしてこのくらいにしたいと思います．
これ以上は別カテゴリで専門家にお聞きになるのが適切でしょう．
頑張ってください．

この回答へのお礼

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

とても参考になりました。
教えていただいた分だけで十分やっていけそうです！

どうもありがとうございました。

お礼日時：2008/09/30 15:19

No.4 回答者： echoes_x86 回答日時： 2008/09/28 18:01

まず，m*n行列(m>nとします)行列XをX=U*S*V^Tと特異値分解します(^Tは転置記号とします)．

ここで, 行列のサイズはUがm*n，Sがn*n，Vがn*n行列です．
UとVは直交行列の一部であり，それぞれXの列空間と行空間の直交基を成します．
直交行列の逆行列は自らの転置ですから，Xの疑似逆行列X^+は
X^+ = V*S^(-1)*U^T
と計算できます．
問題となるのはS^(-1)ですが，これはXの特異値を対角要素とする行列ですから，
対角要素の逆数を取ることで容易に逆行列を計算できます．
ただし，Sの対角要素で零に近いものがある場合は除算を行ってはいけません．
この場合は除算に対応するUとVの列を計算から除外してください．

なお，特異値分解が必要ならばATLASでは対応できません
(ATLASは線型方程式ソルバのみを提供しますので)．
GSLや「Numerical Recipes in C」を用いるのでなければLAPACKのdgesvd_()等が必要になります．

この回答へのお礼

丁寧な解説ありがとうございます。
とても参考になります。

質問なのですが、 (m<n)の場合ならば、
Xの疑似逆行列 X^+ = U^T * S^(-1) * V で求まるのでしょうか？

それと、正方行列の場合は、(m<n)の式でも(m>n)の式でも
どちらからでも求められる、といった感じでしょうか？

すみませんが、教えて下さい。

お礼日時：2008/09/29 15:26

No.3 回答者： MASA_H 回答日時： 2008/09/28 14:36

VC++6で動くよう解説付です。

GSLについての説明ですがなにをしているんのか理解できればATLASにも応用できます。

参考URL：http://www.tmps.org/index.php?GNU%20Scientific%2 …

この回答へのお礼

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

ＶＣ＋＋６も持っているので、参考ＵＲＬの方、読んでみたいと思います。
ありがとうございます！

お礼日時：2008/09/28 17:23

No.2 回答者： echoes_x86 回答日時： 2008/09/28 13:17

こんにちは．

回答No.1のライブラリはWindows環境での実行を前提としていないようです
(利用不能という意味ではありません)．

「逆行列を計算したい」のであれば，LU分解と呼ばれる方法を用いれば良いと思います．
ただし，これは正方行列(50*50の行列)の場合で，ご質問の50*10000の行列の場合には使えません．
この場合はQR分解や特異値分解によって疑似逆行列を計算することになります．
さて，手段の方ですが「Numerical Recipes in C 日本語版」(参考URLの書籍)のライブラリを用いるのが簡単であるようです．

最後になりますが，この種の問題は逆行列そのものが必要かどうかで少し話が変わってきます．
解きたい問題について補足がいただければもう少し詳細な回答を差し上げることができるかもしれません．

参考URL：http://gihyo.jp/book/1993/4-87408-560-1

この回答へのお礼

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

色々と勉強不足ですみません。
助かります。

特異値分解によって擬似逆行列が求まるとのことですが、
Ｘ＝ＵＳＶが基本の形として、Ｘの逆行列はどうやって求めるのでしょうか？
実は、特異値分解の関数は前に頂いたことがあり持っているので
Ｘという行列さえ与えてやれば　Ｕ，Ｓ，Ｖは求められます。

これを使って、長方系の行列の擬似逆行列も、正方行列の逆行列も
求まるということですよね？

Ｕ、Ｓ、Ｖの掛け方でＸの逆行列が求まると思うのですが
その掛け方を是非教えて下さい。
よろしくお願いします。

（ちなみに今回は逆行列そのものを必要としています。）

お礼日時：2008/09/28 17:22

No.1 回答者： MASA_H 回答日時： 2008/09/28 00:03

オープンな実装ならATLASを使ってみるのが一番かも。

ベンダーからBLASの実装を手に入れることがそちらのほうが大抵早いです。スパコンなら大抵用意されてますし、intel系ならIntel Math Kernel Library、AMD系ならAMD Core Math Libraryに含まれています。

参考URL：http://www.mlab.ice.uec.ac.jp/~ej-sib/numerical/ …

この回答へのお礼

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

早速サイトの方みて、導入してみようと思ったのですが
関数ライブラリ使うのが初めてでして導入方法が理解できませんでした。

申し訳ありませんがアドバイスお願いします。

今、ＵＲＬにあるライブラリ構成方法の部分を試してみました。
RANDLIB_V90.tar.gzというファイルを保存したので、
その後のコマンドプロンプトだと思うのですが、それを打つ部分が
わかりませんでした。これってwindowsでも使えるんですよね？

ちなみに私の環境はWinXP
VisualStudio2005を使用しています。

VisualStudio2005で関数使えるようにするまでのアドバイスが頂けると助かります。
よろしくお願いします。

お礼日時：2008/09/28 09:42