2018年6月17日 (日)

最小二乗法の心(2)

前回のブログ を書いた後、最小二乗法について検索してみました。

正規方程式

前記事にも追記していますが、あの「デタラメな」解法には裏付けがあり、式 (1) は正規方程式と呼ばれるものとのことです。

\[ {}^t\!X\boldsymbol{y} = {}^t\!XX\boldsymbol{\alpha} \tag{1} \]

正規方程式は次のように求まります。式 (2) のように両辺が必ずしも一致していない時に近似解を求めることを考えます。

\[ \boldsymbol{y} \fallingdotseq X\boldsymbol{\alpha} \tag{2} \]

式 (2) の両辺の差(の二乗)は式 (3) になります。

\[\begin{align} \|X\boldsymbol{\alpha} - \boldsymbol{y}\|^2 &= ({}^t\!\boldsymbol{\alpha}{}^t\!X - {}^t\!\boldsymbol{y})(X\boldsymbol{\alpha} - \boldsymbol{y})\\ &= {}^t\!\boldsymbol{\alpha}{}^t\!XX\boldsymbol{\alpha} - 2\ {}^t\!\boldsymbol{\alpha}{}^t\!X\boldsymbol{y} +{}^t\!\boldsymbol{y}\boldsymbol{y} \tag{3} \end{align}\]

式 (3) が最小となる時、即ち式 (3) を \(\boldsymbol{\alpha}\) で微分したものが \(0\) となる時、正規方程式 (1) が導かれます。なお、 \(\boldsymbol{\alpha}\) での微分とは勾配( \(\alpha\) による \(\operatorname{grad}\) )のことです。

距離(の2乗)を最小にしたら?

最小二乗法は y 方向の誤差(の2乗)が最小になるような直線を決定するわけですが、データと直線との距離(の2乗)が最小になるように直線を決定してはどうかという件について、幾つか記事が有りました。特に 相関係数に関する一考察 に詳しい考察が有ります。その記事では \(x\),\(y\) の入れ替えやデータの回転に対する回帰直線の性質を考察しています。それらについては、重要性が理解できなかったのですが、終わりの方で、直線との距離(の2乗)が最小になるような直線は「多変量解析の主成分分析に対応しているようである」と有ったのに興味が惹かれました。

他の記事で、最小二乗法で求めた回帰直線の特徴として、以下のことが記されていました。

  • 誤差の平均は \(0\)
  • 分散が最小になる

この特徴が、直線との距離(の2乗)が最小になるような直線の場合はどうなるのか興味が有りますが、自力で解決するのはチョット・・・。

最後に、直線 \(y=ax+b\) とデータ \((x_i,\ y_i)\) との距離の2乗を書いておきます。これが最小になる直線が気になるのでした。

\[ \sum\frac{(y_i - ax_i - b)^2}{(a^2 + 1)} \]

2018年6月12日 (火)

最小二乗法の心

YouTube で見た数学動画からのネタ第2弾です。動画は下記に有ります。

【大学数学】最小二乗法(回帰分析)【確率統計】

\((x_i,\ y_i)\quad i=1,\dots,n\) をデータとしたとき、最小二乗法で直線 \(y=ax+b\) の傾き \(a\) と切片 \(b\) を求めよという問題です。動画では

\[ \sum (y_i - ax_i - b)^2 \tag{1} \]

が最小になるように \(a\), \(b\) を求めています。

ここでは、別の方法で \(a\), \(b\) を求めることにします。まず \(y_i = ax_i + b\quad i=1,\dots,n\) を行列を用いて式 (2) のように書きます。

\[ \begin{pmatrix} y_1\\ \vdots\\ y_n \end{pmatrix} = \begin{pmatrix} x_1 & 1\\ \vdots & \vdots\\ x_n & 1 \end{pmatrix} \begin{pmatrix} a\\ b \end{pmatrix} \tag{2} \]

これを \(\boldsymbol{y}=X\boldsymbol{\alpha}\) と置くことにします。式 (2) から \(a,\ b\) を求めたいのですが、行列 \(X\) が正方行列でないため一般に解は有りません。それでも無理矢理 \(a,\ b\) を求めるために、式 (2) の両辺に \({}^t\!X\) を掛けることで \(\boldsymbol{\alpha}\) に掛かる行列を正方行列にします。

\[\begin{align} {}^t\!X\boldsymbol{y}&={}^t\!XX\boldsymbol{\alpha}\\ \begin{pmatrix} \sum x_iy_i\\ \sum y_i \end{pmatrix} &= \begin{pmatrix} \sum x_ix_i & \sum x_i\\ \sum x_i & n \end{pmatrix} \begin{pmatrix} a\\ b \end{pmatrix} \tag{3} \end{align}\]

そして、式 (3) を解いて \(a,\ b\) を求めると、それが最小二乗法の解になります。

\[ \frac{1}{\sum{x_ix_i}n-\sum{x_i}\sum{x_i}} \begin{pmatrix} n & -\sum{x_i}\\ -\sum{x_i} & \sum{x_ix_i} \end{pmatrix} \begin{pmatrix} \sum{x_iy_i}\\ \sum{y_i} \end{pmatrix} = \begin{pmatrix} a\\ b \end{pmatrix} \]

このデタラメな方法でうまく行く根拠は説明できないのですが、何故か最小二乗法の答えに一致するのです。

追記1:下記に説明が有りました。

最小二乗法の行列表現(単回帰,多変数,多項式)

ところで、式 (1) で \(ax_i+b\) は直線上の点です。したがって、式 (1) は \(y\) 方向の誤差(距離)を考えていることになります。即ち、最小二乗法は \(y\) 方向の誤差が最小になる直線を求めているわけです。しかし、本当の距離はデータから直線に下した垂線の長さなのだから、その長さが最小になる直線を考えるべきだと思うのですが、どうなのでしょうか?

y 方向の誤差 直線とデータとの距離

2018年5月 5日 (土)

ある漸化式の一般項

YouTube で数学の大学入試問題を解説している動画が幾つも有り、そのうちの1つで扱っている問題に興味を持ったので、それについて書きます。動画は次の所に有ります。

YouTube の動画:横浜市立(医)漸化式

問題は、次の漸化式の一般項を求めよというものです。

\[ a_{n+2} - 5a_{n+1} + 6a_n - 6n = 0 \tag{1}\\ a_1 = a_2 = 1 \]

おそらく、動画でやっている方法が普通の解法だと思うのですが、ここでは、行列を用いてこの問題を解いてみることにします。

式 (1) において、 \(6n\) が無ければ、次の式 (2) で行列のn乗を計算すれば良いのですが、 \(6n\) が有るのでそう単純ではありません。一工夫必要になります。

\[ \begin{pmatrix} a_{n+2}\\ a_{n+1} \end{pmatrix} = \begin{pmatrix} 5 & -6\\ 1 & 0 \end{pmatrix}^n \begin{pmatrix} a_2\\ a_1 \end{pmatrix} \tag{2} \]

式 (1) を行列で表すと式 (3) になります。

\[\require{color} \begin{pmatrix} a_{n+2}\\ a_{n+1}\\ b_{n+1}\\ 1 \end{pmatrix} = \begin{pmatrix} 5 & -6 & \textcolor{red}6 & 0\\ 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 1 \end{pmatrix}^n \begin{pmatrix} a_2\\ a_1\\ b_1\\ 1 \end{pmatrix} \tag{3} \]

この行列の右下部分

\[\begin{array}{cc} 1 & 1\\ 0 & 1 \end{array}\]

が工夫した箇所です。この部分は

\[ \begin{pmatrix} n+1\\ 1 \end{pmatrix} = \begin{pmatrix} 1 & 1\\ 0 & 1 \end{pmatrix} \begin{pmatrix} n\\ 1 \end{pmatrix} \]

のように、特定の要素を 1 大きくする働きが有ります。したがって、初期値を適当に選べば、その要素は項番号と同じものになります。これと式 (3) の第1行3列の \(6\) とで \(6n\) の働きをさせることが出来ます。

さて、式 (3) の行列を \(A\) と置きます。 \(A\) のn乗を計算するために、 \(A\) を式 (4) のように分解します。

\[ A = PJP^{-1} \tag{4} \]

\(J\) は対角行列・・・にしたいところなのですが、今回の \(A\) は対角化出来ないのが悩ましいところです。とりあえず、なるべく対角行列に近いものに変形します( Jordan 標準形)。 \(A\) の特性方程式

\[ \det(x - A) = (x - 3)(x - 2)(x - 1)^2 = 0 \]

より、固有値は \(x = 3, 2, 1\) となります。これらの固有値に対応するベクトル \(\boldsymbol{p}_1, \boldsymbol{p}_2, \boldsymbol{p}_3, \boldsymbol{p}_4\) を並べると \(P\) が求まり、式 (5) になります(この辺りの計算が面倒です)。固有値 \(1\) は重複しているので、それに対応するベクトル(独立なもの)も2つ並べています。

\[ P = (\boldsymbol{p}_1 \boldsymbol{p}_2 \boldsymbol{p}_3 \boldsymbol{p}_4) = \begin{pmatrix} 3 & 2 & 6 & 18\\ 1 & 1 & 6 & 12\\ 0 & 0 & 2 & 1\\ 0 & 0 & 0 & 2 \end{pmatrix} \tag{5} \]

この時 \(J\) は式 (6) になります。

\[ J = \begin{pmatrix} 3 & 0 & 0 & 0\\ 0 & 2 & 0 & 0\\ 0 & 0 & 1 & 1\\ 0 & 0 & 0 & 1 \end{pmatrix} \tag{6} \]

以上より \(A^n\) が求まります。

\[\begin{align} A^n &= PJ^nP^{-1}\\ &= P \begin{pmatrix} 3^n & 0 & 0 & 0\\ 0 & 2^n & 0 & 0\\ 0 & 0 & 1 & n\\ 0 & 0 & 0 & 1 \end{pmatrix} P^{-1} \end{align} \tag{7} \]

ここに \(P^{-1}\) は式 (8) となっています。

\[ P^{-1} = \frac{1}{4} \begin{pmatrix} 4 & -8 & 12 & 6\\ -4 & 12 & -24 & -24\\ 0 & 0 & 2 & -1\\ 0 & 0 & 0 & 2 \end{pmatrix} \tag{8} \]

そして、式 (3) で \(a_2 = a_1 = 1, b_1 = 1\) とおいて計算すると一般項は式 (9) となります。ただし、式 (9) で求まるのは \(a_{n+2}\) なので \(a_n\) に書き換えた方が良いかも知れません。

\[ a_{n+2} = \frac{1}{2}(7\cdot 3^{n+1} - 20\cdot 2^{n+1} + 6n + 21) \tag{9} \]

以上の方法は、連立方程式を解いたり、逆行列を求めたりと手間がかかるので、入試問題の解法としては適さないでしょう。その点で YouTube で示された方法が良いと思います。そもそも、ここで行った行列に関する操作は、高校の範囲を越えていますし。しかし、遷移行列という単純なアイデアで解いて行くところは気に入っています。

2018年2月17日 (土)

ドレスト光子の本買いました

ドレスト光子についての本「ドレスト光子」 (ISBN978-4-254-21040-8) を買いました。「ドレスト光子」で検索すると、種々の書籍通販サイトの当該書籍ページが幾つもヒットします。極め付けは当該書籍の閲覧ページです。
「ドレスト光子」の閲覧
ただし、ここで見れるのは第2章のみです。これを見て、面白そうだったので購入しました(印刷の普通の本)。

2章、4章そして付録Aがドレスト光子の質量獲得を理解するための参考になりそうなので、そこを重点的に読んでいます。計算は丁寧に説明されています。丁寧過ぎて、木を見て森を見ずという感は有ります--読む側の能力の問題なのですが。この本を読む目的は、ドレスト光子の質量獲得の理解ですが、そこに到達するには暫く掛かりそうです。取り敢えず、気付いたことを書いてみます。

ドレスト光子の意味

「ドレスト光子」は dressed phton 即ち、ドレスを纏った光子のことです。では、ドレスを纏った光子とはどういう意味でしょうか。それは、以下のようになっています。

ナノ寸法領域での光と電子・正孔のハミルトニアンは式 (1) になります。

\[\begin{align} \hat{H} = &\sum_{\boldsymbol{k}\lambda}\hbar\omega_{\boldsymbol{k}}\hat{a}^\dagger_{\boldsymbol{k}\lambda}\hat{a}_{\boldsymbol{k}\lambda}\\ &+ \sum_{\alpha\gt F,\beta\lt F}(E_\alpha - E_\beta)\hat{b}^\dagger_{\alpha\beta}\hat{b}_{\alpha\beta} + \hat{H}_\mathrm{int} \tag{1} \end{align}\]

ここに \(\hat{a}^\dagger_{\boldsymbol{k}\lambda}\) (\(\hat{a}_{\boldsymbol{k}\lambda}\)) は光子の生成(消滅)演算子、 \(\hat{b}^\dagger_{\alpha\beta}\) (\(\hat{b}_{\alpha\beta}\)) は電子・正孔対の生成(消滅)演算子です。 \(\alpha, \beta\) は電子・正孔対のエネルギー準位、 \(F\) はフェルミ準位です。 \(\hat{H}_\mathrm{int}\) は光子と電子・正孔対の相互作用です。さまざまな状態の光子や電子・正孔対について和を取っていることが、ナノ寸法領域の特徴(だそう)です。

式 (1) をあれこれ変形して行くと、式 (2) のハミルトニアンが得られます。

\[\begin{align} \tilde{H} &= \hat{U}^{-1}\hat{H}\hat{U}\\ &= \sum_{\boldsymbol{k}\lambda}\sum_{\alpha\gt F,\beta\lt F}\left[ \hbar\omega'_\boldsymbol{k}\tilde{a}^\dagger_{\boldsymbol{k}\lambda}\tilde{a}_{\boldsymbol{k}\lambda} + (E'_\alpha - E'_\beta)\tilde{b}^\dagger_{\alpha\beta}\tilde{b}_{\alpha\beta} \right] \tag{2} \end{align}\]

式 (2) では、見かけ上相互作用が無くなっています。 \(a\) や \(b\) の上に付いているのが \(\hat{}\) から \(\tilde{}\) に変わっているに注意して下さい。その \(\tilde{a}\), \(\tilde{b}\) は式 (3), (4) になっています。

\[\begin{align} \tilde{a}_{\boldsymbol{k}\lambda} &= \hat{a}_{\boldsymbol{k}\lambda} - iN_\boldsymbol{k}\sum_{\alpha\gt F,\beta\lt F}\left( \rho^*_{\alpha\beta\lambda}(\boldsymbol{k})\hat{b}_{\alpha\beta} + \rho^*_{\beta\alpha\lambda}(\boldsymbol{k})\hat{b}^\dagger_{\alpha\beta} \right) \tag{3}\\ \tilde{b}_{\alpha\beta} &= \hat{b}_{\alpha\beta} - i\sum_{\boldsymbol{k}\lambda}\left( \rho_{\alpha\beta\lambda}(\boldsymbol{k})\hat{a}_{\boldsymbol{k}\lambda} + \rho^*_{\beta\alpha\lambda}(\boldsymbol{k})\hat{a}^\dagger_{\boldsymbol{k}\lambda} \tag{4} \right) \end{align}\]

式 (3) は右辺第1項の光子が第2項の電子・正孔対の衣を纏っていることを表しています。そして、 \(\boldsymbol{k}\), \(\lambda\) について和を取った \(\sum_{\boldsymbol{k}\lambda}\tilde{a}_{\boldsymbol{k}\lambda}\), \(\sum_{\boldsymbol{k}\lambda}\tilde{a}^\dagger_{\boldsymbol{k}\lambda}\) によって表される光子がドレスト光子です。また、式 (4) は光子と電子・正孔対が入れ替わってほぼ同じ形になっています。即ち、電子・正孔対も光子の衣を纏っていることになります。

グラスファイバー先端のモデル化

グラスファイバー先端のドレスト光子はどのように扱うのかと思っていたら、なんと、単純に原子を1次元格子振動子として計算を進めていました。式 (5) がその1次元格子振動子のハミルトニアンです。

\[ H = \sum_{i=1}^N\frac{\boldsymbol{p}_i^2}{2m_i} + \sum_{i=1}^{N-1}\frac{k}{2}(\boldsymbol{x}_{i+1} - \boldsymbol{x}_i)^2 + \sum_{i=1,N}\frac{k}{2}\boldsymbol{x}_i^2 \tag{5} \]

ここで \(\boldsymbol{x}_i\), \(\boldsymbol{p}_i\), \(m_i\) は \(i\) 番目の原子の(平衡位置からの)変位、運動量、質量、 \(k\) はバネ定数です。

式 (5) にドレスト光子のハミルトニアンやドレスト光子とフォノンとの相互作用などを付加して考察対象となるハミルトニアンが得られるのですが、今回は割愛します。

ところで、私が興味あるのは、 \(A_\mu\) を電磁ポテンシャルとしたとき、式 (6) を導くことです。そのためにはラグランジアン密度を決定する必要が有るのですが、この本からそれを読み解いて行きたいと思っています。

\[ (\Box + m^2)A_\mu = 0 \tag{6} \]

2018年1月14日 (日)

ドレスト光子

光に質量を与えたい -3- で超伝導体内で光が質量を獲得するカラクリが多少分かったわけですが、そもそもの発端は近接場光が質量を獲得することを理解したいということでした。その近接場光に関して新たに資料を見つけました。それが次の3点です(PDF ファイル)。

近接場光による光技術の質的変革
ドレスト光子によるバルク結晶シリコン発光素子
大津教授の講演資料

この中に出て来るキーワードで次のものが目に入りました。

  • ドレスト光子
  • 物質中の励起と光子が結合した場

ドレスト光子とは物質励起の衣をまとった光子という意味です。光子と電子とが結合した状態を表す準粒子がドレスト光子であるとの説明も有りました。

「物質中の励起と光子が結合した場」は、真空が相転移してゲージ場が質量を獲得するということを連想させます。たぶん、近接場光(ドレスト光子)を場の理論として説明するということは既に出来ているのでしょう。近接場光に対して持っていた疑問が、いよいよ解消できるかも知れません。

«転がっている球はいつ止まる?