ベクトル値関数の最適化 (6)

入力もベクトル、出力もベクトルの関数、

f(x) --- (1)

において、近似することをせず、エラーEを、以下で定義します。

E = fT(x)f(x) --- (2)

式(2)を、xで微分してみましょう。

∂E/∂x = 2JTf --- (3)

式(3)をさらに微分します。

2E/∂x2 = 2(∂JT/∂x)f + 2JTJ --- (4)

式(4)の右辺第1項はややこしいのですが、これはf = 0の付近では、ゼロに近くなりますから、第1項を無視してしまいます。すると、式(4)は、

2E/∂x2 = 2JTJ --- (5)

となります。これが、ガウス-ニュートン近似と言われるもので、LM法で使われるものです。
スポンサーサイト

ベクトル値関数の最適化 (5)

今回は、前回紹介した式、

(JTJ + λI) h = -JTf --- (1)

を実際にどう解くか、という話です。

これは連立一次方程式なので、逆行列を求めれば解けますが(正則化項があるので、逆行列は必ずある)、実際には逆行列を求めることはしません。なぜかというと、時間がかかるからです。一般に、変数が多いときは、逆行列をまともにもとめることは推奨されません。

では、どうするか。よくやられるやり方は、Cholesky分解と、共役勾配法ですね。これは専門的になりますので、適当な書籍をご覧ください。最適化を扱っているものであれば、大抵は載っているはずです。

ベクトル値関数の最適化 (4)

前々回で紹介した、一次方程式系、

JTJh + JTf = 0 --- (1)

ですが、ここから、Levenberg–Marquardt法(LM法)に移行することができます。式(1)を移項して、さらに正則化項を加えてやると、

(JTJ + λI) h = -JTf --- (2)

という式が得られます。λは調整パラメタで、これが大きくなると、hは小さくなり、逆もなりたちます。まさに、機械学習における正則化です。

ベクトル値関数の最適化 (3)

前回の続きます。

あ、ここで終わっていてはいけませんでした。これは一次近似解なので、これを繰り返すわけですね。いわゆるニュートン法です。

つまり、hが求まりますから、

xx + h --- (1)

としてやって、同じ計算をして、収束したら(または決められた回数をおこなって)、終了します。

ベクトル値関数の最適化 (2)

ベクトル値を出力に持つ多変数関数の最適化について、以前まとめました。

http://kanouy.blog9.fc2.com/blog-entry-1913.html

ちょっと舌足らずなところがあったので、少し書き換えます。入力もベクトル、出力もベクトルの関数、

f(x) --- (1)

において、このベクトルの長さを最小化することが目的です。このままではどうしようもないので、式(1)をテイラー展開します。一次の項までとると、

f(x + h) = f(x) + Jh --- (2)

Jはヤコビアン(行列)で、以下です。

J = ∂f/∂x --- (3)

エラーE(これはスカラ)を、以下で定義します。これを最小化したいわけです。

E = f(x + h)Tf(x + h) --- (4)

式(4)に式(2)を代入すると、

E = (f(x) + Jh)T(f(x) + Jh) = hTJTJh + 2hTJTf(x) + f(x)Tf(x) --- (5)

式(5)をhで微分して、それをゼロとおきます。

JTJh + JTf(x) = 0 --- (6)

式(6)はhに関する一次方程式系なので、線型代数の方法により、解くことができます。あまり変わらなかった。

湘南工科大学非常勤講師 (10)

湘南工科大学の非常勤講師、線形代数のほうは如何に?

こちらは、ベクトルと行列です。応用上大切な内積をまずやって、そのあと、行列を使って、ベクトルを変換するやりかたを教えます。

計算はプログラミング(Java)でやるのですが、ここでの課題は、次元が低いので(2次元または3次元)、手計算でもできてしまうこと。これだとプログラミングの重要性が伝わりにくいです。実際のところ、手計算でやったような解答を出してくる学生さんもいます。

これは私の出す課題にも問題があるので、ちょっと思案のしどころ。すでに手は打ってありますが。

湘南工科大学非常勤講師 (9)

湘南工科大学の非常勤講師、ガイダンスを含め、5回が終わりました。だんだんペースが掴めてきた。

まず、微積のほうです。プログラミング(Java)にかけるので、結局のところ、級数をやっています。級数計算は、微分が絡みますし(テイラー展開)、実用的、計算もプログラミングでないとできないので(たくさんの項を足し合わせるので)、結構フィットしている気がします。

Chain rule

微分でおなじみの、chain ruleですが、chainという単語に、これまで引っかかっていました。

それが、"Deep Learning (2016)"を読んでいて、この呼び名の意味するところが、なんとなくわかりました。

同書の第6章で、各変数がノードで表されたグラフの微分を扱っているのですが(back propagationで必要)、まさにchain ruleです。ノードをたどるたびに、各ノードの微分を掛けていけばよいし、別の経路があれば、経路ごとに足し合わせればよいわけです。まさにchainです。なるほど。

余因子行列 (3)

逆行列は、掃きだし法で計算するのが基本です。

ただ、例外的に、逆行列=(数式)、として書き下したいことがあります。そのときは、以下の公式があります。

A-1 = adj(A) / det(A) --- (1)

adj(A)は、行列Aの「余因子行列」です。余因子行列は定義の難しいものです。先日、3x3行列の逆行列を計算するプログラムを、式(1)で作ったのですが、コーディングの間違いで、デバッグに時間がかかりました。

湘南工科大学非常勤講師 (8)

湘南工科大学の非常勤講師、今年は本日(2018年4月13日)が最初です。まずはガイダンス。

科目は数学で、線形代数と微積をやるのですが、それにプログラミングをかけるという、多少とも贅沢なカリキュラムです。プログラミング言語はJavaです。

いちおう主担当のN教授と進めかたは協議しましたが、今年からのカリキュラムなので、試行錯誤をしながら進めていくことになるでしょう。試行錯誤は不得手ではありませんし、悪いことでもありません。もしや、教育のアジャイル?学生さんの反応を見ながらやっていくことになるでしょう。

ベクトルでの微分 (3)

機械学習、特にニューラルネットをやっていると、ベクトルでの微分が出てきます。

でも、通常の微積教科書では、こういうのは載っていません。通常は、たとえば、

z = z(x, y) --- (1)
x = x(u, v) --- (2)
y = y(u, v) --- (3)

などとしたときに、

∂z/∂u = ∂z/∂x・∂x/∂u + ∂z/∂y・∂y/∂u --- (4)
∂z/∂v = ∂z/∂x・∂x/∂v + ∂z/∂y・∂y/∂v --- (5)

などと書かれてあります。でも、これだと変数が増えたときに対応できません。そこで、

x = x(x, y) --- (6)
u = u(u, v) --- (7)

と変数をベクトルとしてパッケージ化してやります。そうすると、式(4)(5)は、

dz/du = dz/dx・dx/du --- (8)

と簡潔に書くことができます。変数がいくつ増えても式(8)は変わりません。こういう記法をどんどん増やしたほうがいいと思います。

湘南工科大学非常勤講師 (7)

湘南工科大学の非常勤講師、来年度(2018年4月~)も継続することになったので、先日(3月26日)、研修会に出向きました。

大学の現状などをお聴きしました。大学もいろいろとたいへんですね。いちおう方針はわかったので、その線に沿って協力できればよいな、と思いました。

そのあとは懇親会です。すでに何人の先生がたとは面識ができました。

シラバスは作りましたが、ルーブリックを授業開始前に作らなければなりません...

対称正定値行列

対称正定値行列、日本語での呼び名は多々あれど、英語では、symmetric positive-definite、です。

この行列は、いろいろとよい性質があるのですが(たとえば、LU分解が必ずできる)、では具体的にはどのような行列?

私が知っているのは、共分散行列です。共分散行列Vは、

V = (1/n)ΣppT --- (1)

と書けます。でもΣはややこしいので、design matrixを使います。すなわち、

V = (1/n)ΦTΦ --- (2)

つまり、Φpを行ベクトルとして「パッケージ化」したものです。転置が逆になることに注意。

さて、Vは式(2)から対称行列です。なぜならば、

(ΦTΦ)T = ΦTΦ --- (3)

では、式(2)が正定値であることを確かめましょう。

(x, Vx) = (x, (1/n)ΦTΦx) = (1/n)(Φx, Φx) = (1/n)|Φx|2 --- (4)

式(4)は正ですから、Vは正定値であることが確認できました。

テイラー展開の謎 (3)

テイラー展開、必要があって、復習していました。

応用上は実関数しか必要ないのですが、きちんと理解しようとすると、やはり複素関数にいかなければなりません。つまりは正則関数の理論。ここは実は、あまり得意ではない、というか苦手なところです。

というわけで、志賀浩二先生「複素数30講(1989)」を再読することにしました。四半世紀前に購入したものです。

abc予想

「abc予想」を、日本人数学者(M先生)が証明した、というニュースが駆け巡っています。

これは2012年に発表されたものですが、まったく新しい数学を構築し、それを適用することで、abc予想を証明した、と主張するものです。なので、これが正しいかを検証するには、まず、まったく新しい数学(= Inter-universal Teichmüller Theory: ITU)を習得しなければならないわけです。

暫し前、ある英語ジャーナルに、この検証がなぜ困難を極めているのか、書かれてありました。つまり、もしもこの証明が間違っていたとすれば、せっかく新しい数学を習得したとしても、それは浮かばれないかもしれず、大量の時間を費やすので自身のキャリアを棒に振るリスクがある、ということです。その通りだと思います。

なので、M先生はもちろん素晴らしいに決まっていますが、それを検証しようとした方々も、同等に素晴らしいと思います。

対数の歴史 (3)

先日(2017年12月15日)の、湘南工科大学・非常勤にて、対数の説明をいたしました。

対数関数は、数学書だと、指数関数の逆関数として説明されます。あるいは、1/xの積分として定義されることもある。

ただ、実際のところ、対数というのは、実用的な「計算ツール」として発明されました。16世紀末のことです。これは、航海をするにあたり、大きな数の掛け算をする必要があったことによります。当然、計算機はないので、人手でやった。これがたいへんだということで、掛け算を足し算で代用できないか、と考え出されたのが、対数というわけです。

...などという話を、導入でしたところ、多少興味を持ってもらえた?

平面幾何学 (2)

以前ネタにした、平面幾何学の問題ですが、

http://kanouy.blog9.fc2.com/blog-entry-1516.html

先日(2017年12月4日)、数学関連の飲み会にて、同じ問題を出してみました。

すると、とあるK氏が、「わかった!」と言って、正解を出しました。私の図をじっと見ていただけです。5分くらいでしょうか。ちょっと驚きました。K氏は、数学科のご出身ではありますが、弁護士をされているという、異色の方です。

eの説明 (2)

eの説明を考えるにあたって、eに関する私の知識が、いかにいい加減であるかを再認識いたしました。以下、備忘録。

まず、eの定義ですが、これは、

lim ((eh - 1)/h) = 1 (h→0) --- (1)

を満たす数として定義されます。つまり、x=0における傾きが1である指数関数の底。

そうすると、y=exの微分が、微分の定義式を計算すると、y'=exとなることが導ける、というわけです。合ってますか?

Levenberg–Marquardt法 (2)

Levenberg–Marquardt法(LM法)について、おバカな記事を以前書きました。

http://kanouy.blog9.fc2.com/blog-entry-1904.html

かなりヤバい記事です。すぐにでも削除したいところですが、いたしません。私はこのとき、このようなおバカなことを考えていた、という記録です。そして、それがその後、改善されたということであれば、進歩が確認できるのです。なので、おバカなほうが、進歩は容易。いくつになっても進歩はできる。

eの説明

代表的な無理数に、eとπがあります。無理数はたくさんありますが、つねにこの二つしか登場しないのが、私の長年の疑問です。

ところで、このeとπ、説明するには、πはやさしいですが(円周率なので)、eはどうすればよいのでしょうか?eは極限や微積に絡んで出てくるので、説明が難しいのです。

私のいまの知識では、指数関数y=a^xにおいて、x=0における傾きが1となるようなaをeと定める、というものですが...これより簡単な説明があるのでしょうか?あれば教えてください!

ベクトル値関数の最適化

ベクトル値を出力に持つ多変数関数の最適化をまとめました!あってるかな...

ベクトル値関数を、

f(x) --- (1)

とします。この長さをゼロにすることが、最適化の目標です。

式(1)をテイラー展開します。

f(xx) = f(x) + JΔx --- (2)

ここで、Jはヤコビアンで、以下です。

J = (∂f/∂x) --- (3)

エラーEを、以下で定義します。これがゼロに近づけばよろしい。

E = f(xx)Tf(xx) --- (4)

式(4)に式(2)を代入すると、

E = (f(x) + JΔx)T(f(x) + JΔx) = ΔxTJTJΔx + 2ΔxTJTf(x) + f(x)Tf(x) --- (5)

式(5)をΔxで微分して、それをゼロとおきます。

JTJΔx + JTf(x) = 0 --- (6)

式(6)はΔxに関する一次方程式系なので、線型代数の方法により、解くことができます。オシマイ!

オイラーの定理

多作のオイラー、彼の定理はたくさんありますが、おそらく一番有名なのが、こちら。

e = conθ+ i sinθ --- (1)

神秘的な式です。θ=πとすると、これも有名な(小説の題材にもなった)、

e + 1 = 0 --- (2)

が導けます。基本的な5つの数が、ひとつの方程式で結ばれているという、驚くべきもの。

さて、式(1)の導出ですが、教科書では、テイラー展開を使うものが多いですね。確かにそのとおりなのですが、あまり直観的な説明ではないですね。数式をいじるとそのようになる、という感じなので、いまいち納得感がありません。もっと直観的にわかる導出って、ないんですかね?そうなると、物理的な導出?

Hidden Figures

映画「ドリーム」を観ました。

これは、アメリカでは2016年に公開されたものです。原題は、"Hidden Figures"。表舞台に出なかった人たちの話、という意味でしょうか。

あらすじや評価などは、ネットにいろいろとありますので省きます。私の興味は、コンピュータ黎明期の、ヒトによる手計算の時代です。要するに、エンジニアリング用途の数値計算。ファインマン自伝の時代と被ります。

最初の場面に出てくる4次方程式の解ですが、2次方程式ふたつに因数分解されていて、解がふたつでてくるので、もうひとつの方程式は、解ナシ(複素解)ということでしょうね、おそらく...(よく見ることができなかった)

最後の場面の、「オイラー法」というのは、これは数値積分のやり方のことだと思います。これも、手計算でやっていたとは...

ジャイロ効果 (2)

広瀬茂男先生著「ロボット工学―機械システムのベクトル解析」、「4.7 ジャイロ効果の再考」、に進みます。改訂版で追加されたところです。以下の複雑な式を使います。

N = I (ω'* + ω0×ω) + ω× --- (1)

角速度ωxで回る円板を、角速度ωyで回すと、z軸周りに回転する、というのがジャイロ効果です。ここで、

ω = ωx + ωy --- (2)
ω0 = ωy --- (3)

とおくのがミソ。これで式(1)を計算すると、z軸周りのトルクが計算できます。

でも、結果だけ得たいのであれば、

N = ωy×x --- (4)

としても、同じ結果が得られます。実際のところ、通常は式(4)で求めると思います。敢えて式(1)を使う理由ですが、ジャイロ効果をより明確に理解するためだそうです(改訂版序文より)。

ジャイロ効果

広瀬茂男先生著「ロボット工学―機械システムのベクトル解析」、改訂版を読んでいます。

http://kanouy.blog9.fc2.com/blog-entry-1864.html

「4.6 伸縮回転する角運動量ベクトルの微分」、のところで、複雑な式が出てきます。すなわち、

N = I (ω'* + ω0×ω) + ω× --- (1)

同書では、式(4.42)です。これは複雑です。通常では、

N = Iω' + ω× --- (2)

程度ですけどね(これも複雑だが)。では、なぜ式(1)なのか?これは、次の「4.7 ジャイロ効果の再考」、で使われます。初版にはなかったところです。

湘南工科大学非常勤講師 (6)

湘南工科大学の非常勤講師、金曜日の2コマめです。この日は、8時にオープンする、湘南T-SITE内、スタバ湘南蔦屋書店、に入り、準備をするのが慣例となりました。楽しい時間です。

前回(2017年10月20日)は雨。以下、スタバ店員さん(かわいい女性)との、注文時の会話。

女性:今日は寒いですね~
私 :雨ですからね...
女性:このあたり(T-SITEのこと)にお住みですか~
私 :いえ、もう少し先です。
女性:どちらに行かれるんですか~
私 :湘南工科大というところです。
女性:えっ、なにを教えられてるんですか~
私 :数学です。
女性:....

やはり、「数学」というのは、雑談ではタブーであることを、再認識いたしました。この単語を出すと、円滑な会話がストップ。変人とも思われますから、気をつけよ!

Levenberg–Marquardt法

Levenberg–Marquardt法(LM法)を使っているのですが、よくわからないところがあります。

まず、ベクトルでのコスト関数を設定します。そして、それをパラメタで偏微分します。偏微分の方向にパラメタを変化させてやればよいはず。

しかし、もしも、パラメタが本来と逆の方向を指していたらどうでしょう?これはあり得ます。なぜかというと、ベクトルでのコスト関数は、符号を反転させてやっても、これの二乗和は変わらないからです。対して、偏微分は符号が反転する。

というわけで、もしも、偏微分の方向でコストが増大すれば、逆方向も探索しないといけないと思うのですが、通常の実装ではどうもそうはしないみたいですね。なぜならば、正しい方向に行っても、コストが増大することがあるからです。その場合は、いまより少しだけ進む。それでも増大すれば、更に少し進む。こうしていくと、最初の状態となり、決して逆方向には進みません。この理屈はわかりました。

う~ん、困ったな...おそらく、私がどこか、勘違いをしているのだと思います。

関数の近似 (2)

関数の近似ですが、やはり参照するのは、Numerical Recipes。私は、第二版はハードカバーにて、第三版はオンラインにて。

第5章が、"Evaluation of Functions"として、このテーマに当てられています。とにかくこれを読めばよろしい。

最初に、Lanczosの、"Applied Analysis"が参考文献として挙げられています。これは1956年に書かれた古典です。Lanczonの本は、"The Variational Principles of Mechanics"を持っていて、これを読めば、力学の全てが解るわけですが、読めていないです。残念...

"Applied Analysis"は購入しようと思います。古いものですが、このあたりの内容は、現代でもさほどの変化はないと思います。

関数の近似

いまさらながら、関数の近似技術が重要であると、再認識中。

たとえば、三角関数、sinθですが、スマホ電卓でも正確に計算できます。でも、

sinθ = x - x3 / 6 --- (1)

と近似してやると、数値的な洞察を得ることができるわけです。エンジニアリングでは、具体的な数値、すなわち数感が大切。

式(1)はテイラー展開ですが、関数近似の一般理論みたいなのはないのでしょうか?数学ではなく、応用技術でしょうかね、このあたりは...

KKT条件

金谷健一先生の「これなら分かる最適化数学」で、KKT条件の復習をしています。

7.2 ラグランジュ乗数、のところです。制約式が等式であれば、通常のラグランジュ未定乗数法を使えばよいのですが、制約式が不等式の場合は、難しくなります。要するに、

λigi(x) = 0 --- (1)

が成り立つということですが、これは解釈が難しいですね。サポートベクタのような、具体的な事例を勉強したほうが、よいかもしれませんね。
プロフィール

加納裕(かのうゆたか)

Author:加納裕(かのうゆたか)


[略歴]
1983年3月東京工業大学工学部機械物理工学科卒業
1983年4月(株)図研入社
1987年1月同社退社
1987年2月(株)ソリッドレイ研究所を6名で設立、取締役
1994年3月同社退社
1994年4月(株)スリーディー入社
1996年10月同社取締役
1999年12月上海大学兼務教授
2002年10月同社代表取締役
2009年9月ものつくり大学非常勤講師~現在
2009年10月同社代表退任/退社
2010年1月ソフトキューブ(株)入社~現在(横浜オフィス)
2011年11月甲南大学特別講師
2011年11月関西大学特別講師
2012年11月東京理科大学特別講師
2017年4月湘南工科大学非常勤講師~現在


[業界団体・学会活動]
電気学会・第三期次世代インタラクティブディスプレイ協同研究委員会(幹事)/最先端表現技術利用推進協会・アカデミック部会(旧:三次元映像のフォーラム)(副部会長)/日本バーチャルリアリティ学会ハプティクス研究委員会(委員)/ACM・SIGGRAPH(Professional Member)/情報処理学会(正会員、CVIM会員)/3Dコンソーシアム(賛助会員)/3DBiz研究会(個人賛助会員)/URCF(特別会員)

----------------

前職:立体映像産業推進協議会(幹事)/日本バーチャルリアリティ学会・論文委員会(委員)


[資格]
TOEIC805点
数学検定1級(数理技能)
中型・普自二免許
サッカー4級審判員

最新記事
最新コメント
最新トラックバック
月別アーカイブ
カテゴリ
検索フォーム
RSSリンクの表示
リンク
ブロとも申請フォーム

この人とブロともになる

QRコード
QRコード