疑似逆行列について質問があります．

疑似逆行列について質問があります．
（もしかしたら疑似逆行列の問題で無いかもしれませんが．．．）

例えばF =　Ａτという関係が成り立っているとします．
ただし
τ= [t1,t2,t3,t4,t5,t6,t7,t8]
F = [f1，f2，f3]
A は3行8列の行列式．

ここで0<t1,t2,t3,...,t8 となるような条件のもと
tn(n=1,2,3...8)のノルムの2乗和が最小となるようなτを求めたいとき
どのようにすればよいのでしょうか．



個人的に考えたのは
疑似逆行列を使うと
τ　=Ａ^†F + (I - A^†A)k
のkの部分をいじること何とかなると思ったのですが
どうにも解くことができません．




もしよい方法を知っている方がいましたら教えてください．
また，解法までは分からなくても参考になりそうな考え等が
ありましたらそちらも教えていただけると
大変ありがたいです．





よろしくお願いします

No.1 回答者： ur2c 回答日時： 2010/07/06 14:10

F = Ａτ

0 <τ

の両者を満たす解の存在を確かめます。解はないかもしれません。もしあれば、それらの制約条件のもとで 、「 '」を転置として

τ' τ

を最小化します。計算はたとえば R でできます。

この回答へのお礼

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

τ' τ を最小にするということはつまり
τの各要素の2乗和を最小にするということであると思うのですが，
τ＞0を常に満たすという条件も含まれているのでしょうか．

お礼日時：2010/07/07 12:59

No.2 回答者： alice_44 回答日時： 2010/07/07 00:58

一般逆行列にまつわる用語は、揺らぎが大きくて、

単語を聞いただけでは意味が掴みにくいのですが…
A^+ という記号を使っているところからすると、
貴方の言う「疑似逆行列」とは、私方の方言では
「ムーア・ペンローズ型一般逆行列」のことでしょうか。
最小二乗型 かつ ノルム最小型 かつ 反射型の
一般逆行列のことを、ムーア・ペンローズ型と言います。
任意の行列に対して唯一つ存在するので、よく使われます。

ムーア・ペンローズ型一般逆行列も、ノルム最小型
一般逆行列の一種なので、一次方程式に解が在る場合には、
その中でノルム最小のものを与えます。
k をいじらなくても、τ = (A^+)F だけで ok です。

ただし、F = Aτ が常に解を持つとは限らないので、
注意が必要です。解が存在するための必要十分条件は、
A と F の列を並べた３行９列の行列の rank が
rank A と一致していることです。

この回答へのお礼

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

ただτ = (A^+)F のみだと
τの要素[t1,t2,t3,t4...t8]が正の値だけでなく負の値をとる可能性が
残ってしまいます．
そのような場合は各要素を正の値あるいは0にするにはどのようにすればよいのでしょうか？

お礼日時：2010/07/07 12:55

No.3 回答者： ur2c 回答日時： 2010/07/07 17:17

> τ＞0を常に満たすという条件も含まれているのでしょうか．

はい。ANo.1 にある「それらの制約条件」は「F = Ａτ, 0<τ」を指しますから。実際に計算するには、それらの制約を満たす、いわゆる実行可能解を 1 つ、初期値にとれることが必要です。

Moore-Penrose 一般化逆行列は 0<τ を考慮しませんので、この問題には不向きです。

質問は数学の問題と言うより、現象のモデル化の問題のように見えます。今の定式化では実行可能解の存在が保証されないので、不自然な感じを受けます。現象自体の説明があれば、もう少し適切な助言ができるかもしれません。

この回答へのお礼

丁寧な回答ありがとうございます．
非常に助かります．


現象については説明が難しいのですが，力の再分配を考えています．
力ＦをＡ＝[a1,a2,a3...a8]の方向の力の合力で表わすために必要な各方向anに加える力tnで表わそうというものです．


回答を観ているうちになんとなく今更ながら論点がつかめてきました．
私が考えていた問題とはつまり
目的関数　z = [F-Aτ]^T[F-Aτ]
を最小にするtを求めよ．
ただし
τ＞0
f1=a11*t1+a21*t2+a31*t3...+a81*t8
f2=a12*t1+a22*t2+a32*t3...+a82*t8
f3=a13*t1+a23*t2+a33*t3...+a83*t8
（この3行はF=Aτを表わす）
※このときaijは既知の値

という非線形計画問題ということです．

ただ，私は実は計画問題というものには触れたこともなくどのようにすればよいか全くわかりません．
このような問題は解けるものなのでしょうか．

お礼日時：2010/07/09 19:42

No.4 回答者： ur2c 回答日時： 2010/07/09 22:28

> 目的関数　z = [F-At]^T[F-At] を最小にするtを求めよ．

F = A t なので z = 0 です。

> 力の再分配

だとすると、定式化の第 1 の要点は、0 = t_j を許すか許さないかだと思います。つまり a_j 方向に張った支線が、たるんでも良いかどうか。

たるんで良いなら簡単で、Σ t_j を最小化、つまり a_j 方向の力の和を最小にすることにすれば、問題は線形計画法に落ちます。そして解は、どれか 3 本が張って残りはたるんでいるというものになるか、あるいは存在しません。t' t を最小化でも、大差ないです。

Q1. たるんで良いですか、だめですか？

たるんではいけないなら、話はちょい難しくなります。

> このような問題は解けるものなのでしょうか．

解があれば、解けます。第 2 の要点はそこ、つまり、解が存在しない場合にどうするか、ということです。「ないよ」で済むのか。済まなくて、なくても、できるだけ F に近いものを A t で実現しなければならないのか。

Q2. なければないですみますか、だめですか？

> どのようにすればよいか

一番簡単な場合、つまり「0 = t_j を許し、解がなければないで済む」なら、すぐです。

Q3. 数値計算 software は何を使ってますか？

この回答へのお礼

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

ur2cさんの質問に対する返答ですが
まず

＞Q1. たるんで良いですか、だめですか？
については”たるんでもよい”つまり条件式はti>=0の方が適切でした．
ありがとうございます．

>Q2. なければないですみますか、だめですか？
”できるだけ F に近いものを A t で実現しなければなりません”
ただ，今回の場合，幾何学的に必ず解が存在するように力を分配するので
原則的には解が存在するのではないかと考えています．


>Q3. 数値計算 software は何を使ってますか？
C言語によるプログラミングによって実現しようと考えています．
（数学的に特殊な関数の入ったヘッダーファイルを読み込む予定もありません．）
なので基本的に自ら立式し，解の形まで計算する必要があります．









最後に質問ですが，
目標関数をΣ t_j の最小化問題に置き換えることができるとありますが，
(Σti)^2<Σ(t_i)^2
となりΣ t_j解が常にΣ(t_i)^2の最小であることを示すことは難しいように感じました．
（自分も計算してみましたが，できませんでした．）
これは証明できうるものなのでしょうか．
もし証明できるのならそちらについても教えていただけると嬉しいです．

お礼日時：2010/07/10 11:01

No.5 回答者： ur2c 回答日時： 2010/07/10 11:46

では、F = A t, 0 =< t の解は存在する、とします。

> 目標関数をΣ t_j の最小化問題に置き換えることができるとありますが，

いえ、置き換えることができるわけではありません。言いたいのは「目的関数を Σ t_i^2 に固定する理由はまだない」ということです。

各支線を同じ太さでなるべく細くしたいなら、各支線にかかる力の最大値を最小化すれば良いので、目的関数は Σ max_i t_i みたいなものでしょう。作用点にかかる力の和をなるべく小さくしたいなら Σ t_i になります。目的関数 Σ t_i^2 はどういう（設計上の？）要請から出て来るのでしょうか？

> C言語によるプログラミングによって実現しようと考えています．
>（数学的に特殊な関数の入ったヘッダーファイルを読み込む予定もありません．）

いきなり C で書くのは労力が無駄です。普通は、まずはもっと高級な言語で楽に prototype を作り、満足なものができた後に必要なら、しかたなく C に書き換えます。

この回答へのお礼

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

>各支線を同じ太さでなるべく細くしたいなら、各支線にかかる力の最大値を最小化すれば良いので、
>目的関数は Σ max_i t_i みたいなものでしょう。作用点にかかる力の和をなるべく小さくしたいなら
>Σ t_i になります。目的関数 Σ t_i^2 はどういう（設計上の？）要請から出て来るのでしょうか？

なるほど，全くその通りだと思います．
私は目的関数 Σ t_i^2を最小化することはエネルギー的に小さくなることと同義と思っていたので
最もスマートな力の再分配方法だと勘違いしていました．
言われてみればΣ t_i^2では作用点に加わる力の方向と反対方向を向く支線にも力が分配されるためあまりよい方法ではありません．

今回は作用点にかかる力の和をなるべく小さくしたいのでΣ t_i にしようと思います.




なお，Ｃ言語で書くのは同じくＣ言語を使って計算結果の描画を行うプログラムを作成したためそちらを利用しようと考えているからです．（私がＣ言語とＨＴＭLしか使えないというのもありますが．．．）

お礼日時：2010/07/12 15:10

No.6 回答者： ur2c 回答日時： 2010/07/10 19:30

訂正です。

誤: Σ max_i t_i
正: max_i t_i

No.7 ベストアンサー 回答者： ur2c 回答日時： 2010/07/12 20:04

A x = y, 0 =< x の解は存在する。

[... 1 ...] x を最小化する。A は (3,8)。R でやります。

各列の Euclidean norm が 1 である行列 A を作る。

> set.seed(123)
> A0 <- matrix(runif(3*8,-1,1),nrow=3,ncol=8)
> A <- A0 %*% diag((diag((t(A0) %*% A0)))^(-1/2))

A =
-0.575 0.518 0.071 -0.094 0.403 0.607 -0.276 0.351
0.780 0.595 0.989 0.990 0.165 -0.385 0.730 0.256
-0.246 -0.614 0.130 -0.101 -0.900 -0.695 0.625 0.901

0 =< y として良い。

> (y <- c(round(abs(runif(3,0,9)))))
[1] 6 6 5

解く。

> library(boot)
> x <- simplex(rep(1,8),A3=A,b3=y)
> x$solved
[1] 1

なので解けてる。0 以外の解は

> x$soln[x$soln != 0]
x2 x3 x8
5.3335022 0.4952125 9.1188452

目的関数の値は

> x$value
b
14.94756

A x = y を確認する。

> as.vector(A %*% x$soln)
[1] 6 6 5


高水準言語なら簡単にできるので、1 つは使えるようにしておくべきです。数式を扱うなら Sage とか Maxima とか。数値だけなら R とか Scilab とか。

> Ｃ言語を使って計算結果の描画を行うプログラムを作成したためそちらを利用しようと考えている

描画も高水準言語で作り直した方が、あとあとメンテが楽です。

> 私がＣ言語とＨＴＭLしか使えないというのもありますが

C で書く正当な理由は「計算はマイコン」とか「一瞬を争う実行」しか思いつきません。

この回答へのお礼

丁寧な回答ありがとうございます．


これを参考に実際にプログラムを書いていこうと思います．

また，アドバイスを参考に高水準の言語も勉強していきたいと思います．
本当にありがとうございました．

お礼日時：2010/07/21 19:01

No.8 回答者： ur2c 回答日時： 2010/07/13 22:36

訂正です。

機能は同じですけど。

誤: (y <- c(round(abs(runif(3,0,9)))))
正: (y <- c(round(runif(3,0,9))))