指数関数を含む連立方程式：効率的な数値解法

X : 未知の n×m 実行列 (2 <= m < n)
A : 既知の m×m 実正則行列（m次正則行列)
と置くとき，
X = e(X) A …(1)
が成り立つことがわかっています．但し，e(X)は，Xの各要素xをその指数関数exp(x)で置き換えたn×m 実行列を表すこととします．式(1)の右辺はe(X)とAの積です．

このとき，式(1)をXについて解きたいと考えています．恐らく，代数的に解くことは無理で，数値解法（数値アルゴリズム）を利用するほかないと思います．

この問題の場合，どのような数値解法が効率的でしょうか？
数値解法に疎いので，アドバイスを頂ければ嬉しいです．

No.1 ベストアンサー 回答者： sugakusya 回答日時： 2009/08/29 01:41

Y(X) = X － e(X) A

とYを定義すると、Yが0行列になるようなXが解と言えるので、

f：R^(mn)∋Y→R
f(Y)= Σ{(yij)^2}
(Yの各成分を２乗したものの、総和)

を評価関数として最急降下法は使えますし、簡単です。もっと効率的なものだと非線形共役勾配法(←アルゴリズムは知らないので使えるか不明。名前的に使えそう(笑))、ニュートン法、などが使えると思います。ただ、これより効率のいい方法があるのかも知れません。あくまで、最適化は少しかじりましたが、非線形方程式の解法を勉強したわけではない私の知りうる範囲です。

勉強するところから始めるなら最急降下法がもっとも効率いいです。

この回答へのお礼

早速，貴重なアドバイスを頂き，ありがとうございます。大変参考になりました。

関数f(Y)の最小化によって式(1)を解く方法を，実装できました（勉強を兼ねて最急降下法を使いたかったのですが，今回は時間がないため，統計ソフトRに実装されていた準ニュートン法を使いました）。

そこで，Aを乱数で作り，式(1)を解く数値実験を繰り返し行ったところ，以下のように観察されました。

○ 方程式(1)には，Xの解がある場合とない場合がある。
○ Xが複数の解をもつ場合（初期値によって数値解が異なる場合）もある。
○ Xの解には，（下記の例のように）全ての行ベクトルが同じものと，2種類の行ベクトルから成るものがある。

なぜこのような結果になったのか，理解することが，次の課題です。以下は数値解の一例です。

＜例＞
N = 5, M = 3，
A =
-0.7145957 0.1991183 0.4051587
-0.3763429 -0.6322557 0.3075774
-0.7828631 -0.6872320 -0.1921799
の場合，Xの数値解は初期値を変えても1つしか見つからず，
X =
-1.188371 -0.9125817 0.04576733
-1.188370 -0.9125825 0.04576654
-1.188371 -0.9125815 0.04576673
-1.188371 -0.9125817 0.04576708
-1.188371 -0.9125827 0.04576738
でした．このXに対するf(Y) = 6e-12

以上，報告させて頂きます。長くなってすみません。
重ねて，ありがとうございました。

お礼日時：2009/08/29 03:54

No.7 回答者： noname#101087 回答日時： 2009/08/31 16:11

2×2 行列 A の蛇足。

「逐次代入」だと複数解をうまくつかめず。
Newton ならアッという間に収束するので、初期値を振れば見つかります。

(EXP(X1) を e1 などと略記)
　a11e1 + a12e2 - x1 = y1
　a21e1 + a22e2 - x2 = y2
　　　　↓　微分
　(a11-1)e1*dx1 + a12e2*dx2 = dy1
　a21e1*dx1 + (a22-1)e2*dx2 = dy2
で勘定。(スプレッドシート / シート関数のみ使用。y1, y2 →0 )

一例だけでも。
|0.15　　0.10|
|0.10　　0.20| = A
　　　　↓　数回で収束。
{x1, x2} = {0.3779, 0.4640} & {1.514, 2.119}
　

この回答へのお礼

さらに具体例を示して頂き，ありがとうございます。
ニュートン法，準ニュートン法のように，微分を用いた最適化法でも，初期値を振れば，複数の解を見つけることができるのですね。
数回で収束とは速いですね！

お礼日時：2009/09/09 05:58

No.6 回答者： warumx 回答日時： 2009/08/31 11:13

#2の者です。

>いくつ含まれるかは，準ニュートン法における初期値に依存すると
>考えられます。

そうですね。解に関する情報があれば初期値の設定が楽になり
希望の解が得られやすくなります。

>"X()=零行列"という表現は，"X=零行列"という意味でしょうか。

そのとおりです。括弧()は無い方が分かりやすかったですね。

>また，この例でe(X)*Aは5×3行列なので，行ベクトルは5個かと
>存じます。3つとされているのは，列ベクトルのことでしょうか。

すいません、記入ミスです。5つの行ベクトルですね。

統計ソフトRで答えを出しておられるのでどうかと思いましたが、
Excel Solverも結構簡単に最適解を求めることが出来ますので
それのご紹介のつもりでした。

この回答へのお礼

ご丁寧に補足を頂き，ありがとうございます。
疑問点について，理解できました。
今後はExcel Solverも使ってみます！

お礼日時：2009/09/09 05:51

No.5 回答者： noname#101087 回答日時： 2009/08/30 12:54

とても多次元化まで見通せませんけど。

まずは、2×2 行列 A から。

>A : 既知の m×m 実正則行列（m次正則行列)

・非正則でも解が存在する。
　一例。
　　|0.1　0.1|
　　|0.1　0.1|
　これは、一次元における x = 0.2e^(x) と同じ解をもつ。
　実際、2×2 で「逐次代入」していくと収束し、一次元と一致。
　(自明だが)
　

この回答へのお礼

2次元の場合に拡張して頂き，ありがとうございます．
Aが正則でなくても解は存在するのですね．
Aが正則であるという条件は，私が直面した場面では満たされるので，念のためと思って質問の文章に含めた条件です．

挙げて頂いた例を，質問の文章における，n = 1, m = 2の場合の一例と捉えて，代数的に確認させて頂きました。
X = (x1, x2) と置くと，
e(X) A = (e^x1, e^x2) A = 0.1 (e^x1 + e^x2, e^x1 + e^x2)
より，方程式(1)⇔
(x1, x2) = 0.1 (e^x1 + e^x2, e^x1 + e^x2) …(1')

方程式 : x = 0.2 e^(x)の片方の解をxsと置くとき，
xs = 0.2 e^(xs) = 0.1 (e^xs + e^xs) だから，
X = (x1, x2) = (xs, xs) は確かに式(1')を満たす．

お礼日時：2009/08/30 21:56

No.4 回答者： noname#101087 回答日時： 2009/08/29 17:30

1 次元の蛇足。

>ただし解があれば 2点(二重根はダブルカウント)なので、両方向(x→y, y→x)から攻めます。

これは、0＜a≦1/e の場合でしたね。a＜0 なら解は一つ。
　

この回答へのお礼

ご丁寧にありがとうございます．以下のように，確認できました．

g(x) = x - a*e^x と置くと，
g'(x) = 1 - a*e^x
g''(x) = - a*e^x

I) a > 0 の場合
○ g'(x)は単調減少関数．x = log(1/a)のときg'(x) = 0となり，g(x)は極大かつ最大の値log(1/a) - 1を取る．
○ x→-∞のときg(x)→-∞
○ x→+∞のときg(x)→-∞
以上より，方程式(2) : g(x) = 0の解の存在条件はlog(1/a) - 1 >=0 つまり 0 < a <= 1/e（等号成立時のみ二重根）

II) a <= 0 の場合
○ x > 0 ⇒ g(x) > 0より，区間(0,-∞]にg(x)=0の解はない．g'(x) > 0よりg(x)は単調増加関数であり，区間(-∞,0]におけるg(x)の最大値はg(0) = - a (>= 0)．
○ x→-∞のときg(x)→-∞
以上より，方程式(2) : g(x) = 0の解は1つ．

お礼日時：2009/08/30 21:34

No.3 回答者： noname#101087 回答日時： 2009/08/29 14:26

迂遠ながら、1 次元の場合から。

　y = a*e^x　　　…(2)
の不動点 y = x を求めるには、「逐次代入法」が安直。
収束は Newton 法より遅いけど、計算式は単純。

ただし解があれば 2点(二重根はダブルカウント)なので、両方向(x→y, y→x)から攻めます。
多次元に拡張するには？
　

この回答へのお礼

貴重な回答を頂き，ありがとうございます。
方程式(1)は，右辺が与える写像の不動点を求める方程式と捉えられるのですね！気づきませんでした。
逐次代入法についても，思い出すことができました。

お礼日時：2009/08/30 20:48

No.2 回答者： warumx 回答日時： 2009/08/29 12:48

>代数的に解くことは無理で，数値解法（数値アルゴリズム）を利用す

>るほかないと思います．

そのとおりと思います。
非線形の連立方程式なので、既に皆さんが言われるような解法で
解くのが普通でしょう。物理的な意味を考慮すると効率的な解法が
あるのかも知れませんが当方には分かりませんね。

ということで、提示された例をExcel Solverで解いて見ました。
見掛けは15個の変数の連立方程式ですが、実際は3変数の
連立方程式が5組あるのと同等なので、行列X()の1行の3個の変数
x11,x12,x13について解いた結果は以下のとおりです。
　x11=-1.188371126
　x12=-0.912582511
　x13=0.045767108
質問者様の結果に一致してます。
因みに、ソルバーの解法は準ニュートン法に指定してます。

尚、質問者様の考察について簡単に言及させていただきます。

>○ Xの解には，（下記の例のように）全ての行ベクトルが同じもの
> と，2種類の行ベクトルから成るものがある。

この例では、X()=零行列のとき、e(X)*Aの3つの行ベクトルは同一に
なりますので、同一の解を持つことは容易に予想されますね。
物理的な意味は当方には分かりませんが偶然とは思えません。

Excel Solverは非線形方程式でもちょっと解くには手軽で便利です。
以上ご参考まで。

この回答へのお礼

貴重なアドバイスを頂き,また計算の再現までして頂き,ありがとうございます。とても参考になりました。Excel Solverは，便利そうですね。

> 見掛けは15個の変数の連立方程式ですが、実際は3変数の
> 連立方程式が5組あるのと同等なので
私が見落としていた，重要なポイントかと存じます。方程式(1)は，m次元行ベクトルxに関する方程式：
x = e(x) A …(2)
を，n組まとめて表現したものに過ぎないですね。
従って，方程式(2)がs個の解を持つならば，方程式(1)の解Xには，高々s個の異なる行ベクトルが含まれることになります。いくつ含まれるかは，準ニュートン法における初期値に依存すると考えられます。

ところで，頂いたコメントうちの一文：
> この例では、X()=零行列のとき、e(X)*Aの3つの行ベクトルは同一に
> なりますので、同一の解を持つことは容易に予想されますね。
について，"X()=零行列"という表現は，"X=零行列"という意味でしょうか。また，この例でe(X)*Aは5×3行列なので，行ベクトルは5個かと存じます。3つとされているのは，列ベクトルのことでしょうか。

お礼日時：2009/08/30 20:40