# Math Methods for Physicists <br> Final report : mutual information for single neuron model
> 109022122呂秉宬

---

這篇報告裡主要想討論的是在什麼樣的條件下神經元的輸入端$x_i$和輸出端$y$之間的 mutual information 會有最大值，這個條件會和Hebb learning得到一樣的條件嗎？最後我們會用電腦模擬Hebb learnig，看看用Hebb learning長出來的結果是怎麼樣。

---
## mutual information
1. Mutual information (entropy maximize)
說到相互資訊我們就要先說到Shannon entropy，Shannon entropy和我們在熱統計一開始學到的entropy有一點不一樣，他的定義是$S=\sum_i P_i\log(P_i)$
當裡面的對數是以二當作底數的時候Shannon entropy還有一個特別的意義，那就是：我平均最少需要多少位元才能夠表示系統處在的狀態。我們可以想一個簡單的例子，想像一個箱子裡有四顆白球、兩顆紅球及兩顆黑球，若是我們想紀錄抽出來的球是什麼顏色我們要怎麼做？因為我們已經知道白球最有可能出現，所以第一個位元就可以用來記錄結果是不是白球（1是白球0則不是），如果不是白球的話木就需要用第二個位元來紀錄結果 ex(1是紅球0是黑球)，所以在這個例子中我們平均需要1.5個位元來儲存抽出來的結果，如果我們說抽出來的球是一個系統的話那麼這個系統的$S=\sum_i P_i\log(P_i)=\frac{1}{2}\times1+\frac{1}{4}\times2+\frac{1}{4}\times2=1.5$
在上面的這個例子我們可以看到我們需要用的位元確實$\geq S$。
現在我們可以來看一下相互資訊 $I$ 到底是什麼了。
兩個不同的隨機變數之間的$I(x;y)$ 是
:::info
\begin{gather*}
I(x;y)=h(x)-h(x\mid y)
\end{gather*}
我們可以用需要的最少位元來理解$I(x;y)$，等號右邊的第一項是紀錄 $x$ 所需要用到的位元數，第二項則是在知道 $y$ 的情況下紀錄 $x$ 需要用到的位元，若是我們透過 $y$ 就有辦法完全得知 $x$ 的結果，那麼 $h(x\mid y)=0$ 因此 $x$ 和 $y$ 之間的相互資訊就是 $h(x)$，相反的如果已知 $y$ 對知道 $x$ 完全沒有幫助 $h(x\mid y)=h(x)$，那麼$x$ 和 $y$ 之間的相互資訊就是$I(x;y)=0$
:::

其中$h(x)$是$x$這個變數自己的Shannon entropy，$h(x\mid y)$ 則是在我們已經知道$y$的情況下$x$的Shannon entropy，。若這兩個變數都是連續的話，那麼原本的加總就會變成積分
\begin{gather*}
h(x)=-\int p(x,y)\log(p(x))dxdy\\
h(x\mid y)=-\int p(x,y)\log(\frac{p(x,y)}{p(y)})dxdy
\end{gather*}
:::success
接著我們也將$I(x;y)$ 寫成積分的形式
\begin{gather*}
I(x;y)=\int p(x,y)\log(\frac{p(x,y)}{p(x)p(y)})dxdy
\end{gather*}
可以注意到mutual information 在交換 $x$ 跟 $y$ 之後還是一樣$I(x;y)=I(y;x)$，和KL divergence不一樣。
:::
---

2. Hebb learning and oja rule
Hebb learning 和 oja rule 都是在描述神經在成長的過程會發生什麼事
（圖片仿自林秀豪老師的講義[mutual information](https://elearn.nthu.edu.tw/pluginfile.php/515878/mod_resource/content/1/2023-0410-mutual-information.pdf)）
圖片中左邊的 $x_i$ 是神經元的輸入訊號， $W_i$ 是輸入端的權重 synaptic weight， $y$ 是輸出端，N則是雜訊。輸出端和輸入端的關係可以用這條式子來表示 $\tau\frac{dy}{dt}=-y+G(\sum_ix_iW_i)+N$
因為神經反應的速度應該要非常快，所以這條式子會快速地達成平衡 $\frac{dy}{dt}=0$ 達到平衡的時候（假裝G=1）
\begin{gather*}
y=\left(\sum_ix_iW_i\right)+N
\end{gather*}
:::info
除此之外每一個輸入端的權重$W_i$也是會隨著時間改變
\begin{gather*}
\tau_w\frac{dW_i}{dt}=x_i y
\end{gather*}
這條式子就是Hebb learning，因為神經反應的時間比神經成長的時間短很多，通常 $\tau\ll\tau_w$，所以權重的變化要考慮一段時間的平均
:::
現在我們可以把平衡後的 $y$ 帶進 Hebb learning 裡，我們會得到
\begin{gather*}
\tau_w\frac{dW_i}{dt}=\left(\sum_jx_ix_jW_j\right)+Nx_i\\
\tau_w\frac{d\bf{W}}{dt}=\hat{Q}\textbf{W} , \hat{Q}{ij}=\left<x_ix_j\right>
\end{gather*}
第二行的 $W_i$ 寫成向量 $\textbf{W}=(W_1,W_2,\cdots)$，還有因為權重變化需要的時間比較長，所以快速變化的輸入端的訊號及雜訊都要取平均。$\left<N\textbf{x}\right>$ 在第二行會消失是因為我們假設雜訊和輸入端的訊號是獨立的且雜訊的平均值為零 $\left<N\textbf{x}\right>=0$，從上面的 Hebb learning 我們就能夠解出輸入端的權重會如何成長，將上面的矩陣對角化我們可以得到
\begin{gather*}
\textbf{W}i=C_iexp\left(\frac{\lambda_i}{\tau_w}t\right)\textbf{p}_{i}
\end{gather*}
$\textbf{P}_{i}$ 是 $\hat{Q}$ 的第$i$個歸一化後的eigenvector，$\lambda_i$ 則是第 $i$ 個 eigenvalue。我們可以看到經過很長的時間後每一個特徵值大於零的特徵向量都會爆掉，但是實際上因為資源是有限的，所以神經元並不會無限制地長下去，Oja rule就是在描述資源有限的情況下神經元發展的上限，這邊我們簡單的假設 $1=\sum_i Wi^2$ ，這樣神經元的發展就能夠得到限制了。在有上限的情況下我們可以看出來到最後會剩下的 eigenvector 只剩對應到最大的 eigenvalue 的那一個特徵向量。

---

3. Dose Hebb learning leads to maximize mutual information?
現在我們大概知道了神經元生長的方式，是時候讓我們回到一開始的問題：Hebb learning 真的會將輸入端和輸出端的相互資訊最大化嗎？
輸入端和輸出端之間的相互資訊是$I(y;\textbf{x})=h(y)-h(y\mid \textbf{x})$
我們假設這些浮動的訊息都是常態分佈 $p(x)=\frac{1}{\sqrt{2\pi\sigma^2}}exp\left(\frac{-(x-\mu)^2}{2\sigma^2}\right)$
:::success
常態分佈的 Shannon entropy 為
\begin{gather*}
h(x)=\int\frac{1}{\sqrt{2\pi\sigma^2}}exp\left(\frac{-(x-\mu)^2}{2\sigma^2}\right)\left(\frac{1}{2}\log\left(2\pi\sigma^2\right)+\frac{(x-\mu)^2}{2\sigma^2}\right)dx\\
h(x)=\frac{1}{2}\log\left(2\pi\sigma^2\right)+\frac{1}{2}
\end{gather*}
所以 $h(y)=\frac{1}{2}\log(2\pi\sigma_y^2)+\frac{1}{2}$。 $h(y\mid\textbf{x})$ 是已經知道 $\textbf{x}$ 的情況下$y$ 的 Shannon entropy， 其實那就是雜訊的 Shannon entropy 若雜訊也是常態分佈的話 $h(y\mid \textbf{x})=\frac{1}{2}\log(2\pi\sigma_N^2)+\frac{1}{2}$。如此一來輸出端和輸入端的相互資訊
\begin{gather*}
I(y ;\textbf{x})=\log\left(\frac{\sigma_y}{\sigma_N}\right)
\end{gather*}
:::
因為雜訊$N$不是成長的過程可以控制的，所以若是神經生長的自然機制會最大化輸入端和輸出端之間的相互資訊，那麼唯一的可能就是最大化 $\sigma_y$，假設所有輸入端的訊號和雜訊都在0附近浮動，所以$\left<y\right>=0 , \sigma_y^2=\left<y^2\right>$

\begin{gather*}
\sigma_y^2=\sum_{ij}W_iW_j\left<x_ix_j\right>+2\sum_iW_i\left<Nx_i\right>\\
\sigma_y^2=\textbf{W}^T\hat{Q}\textbf{W} , Q_{ij}=\left<x_ix_j\right>
\end{gather*}
這邊除了丟掉平均值為零的部分還將$\left<Nx_i\right>$也丟掉了，因為輸入端的訊號應該要跟雜訊無關，所以$\left<Nx_i\right>＝\left<N\right>\left<x_i\right>=0$，那在什麼樣的突觸權重下$\sigma_y^2$會有最大值呢？我們可以用 lagrange multipliers 來找極值，對於突觸權重來說，因為不能無限生長，所以$W_i$的constrain是$1=\sum_iW_i^2$

\begin{gather*}
\frac{\partial(\textbf{W}^T\hat{Q}\textbf{W}-\alpha\textbf{W}^T\textbf{W})}{\partial W_i}=0\\
\frac{\partial \textbf{W}^T}{\partial W_i}\hat{Q}\textbf{W}+\textbf{W}^T\hat{Q}\frac{\partial\textbf{W}}{\partial W_i}-\alpha\left(\frac{\partial \textbf{W}^T}{\partial W_i}\textbf{W}+\textbf{W}^T\frac{\partial \bf{W}}{\partial W_i}\right)=0\\
\end{gather*}
:::info
經過一連串的整理後我們會得到
\begin{gather*}
\hat{Q}\bf{W}=\alpha\bf{W} , Q_{ij}=\left<x_ix_j\right>
\end{gather*}
之前Hebb learning 告訴我們最後穩定時$\bf{W}$會是$\hat{Q}$ 的其中一個特徵向量，這和相互資訊最大化所需要的條件是一樣的，所以Hebb learning 的確可以將輸出和輸入端的相互訊息最大化。
:::

---

## Numerical Simulation
我是用matlab來迭代，下面的程式碼是我主要的迭代內容，除此之外為了防止突觸權重爆掉我在每次迭代後都有將他歸一化。

y(i)=x(1)*w1(i)+x(2)*w2(i)+x(3);
w1(i+1)=w1(i)+y(i)*x1(i)*dt;
w2(i+1)=w2(i)+y(i)*x2(i)*dt;
wlong=(w1(i+1)^2+w2(i+1)^2)^(0.5);
w1(i+1)=w1(i+1)/wlong;
w2(i+1)=w2(i+1)/wlong;

除了迭代的內容之外，如何做出隨機變數也是重要的，因為$Q_{ij}=\left<x_ix_j\right>$會直接影響權重的生長方式，我是透過matlab內建的函數來當作我的隨機變數

    x(1)=2*randn;
    x(2)=randn;
    x(3)=randn; %x(3)是雜訊N

這種隨機變數的取法會導致$x_1,x_2$之間毫無相關，所以$Q_{ij}=\left<x_ix_j\right>$的特徵向量就會是(1,0) or(0,1)，從下面這張圖我們能夠看到結果的確是權重只往其中一個方向發展。這邊還有一個小細節值得注意，那就是我們要讓$x_1,x_2$的標準差有所不同，如果一樣的話$Q_{ij}$就會是單位矩陣，最後權重也會隨意的震盪。
在下面的程式裡我故意讓兩個輸入的資訊彼此之間有一定的關聯，最後權重長得就很平均了。

    s=randn;
    x(1)=0.8*s+0.2*randn;
    x(2)=0.8*s+0.2*randn;
    x(3)=randn; %x(3)是雜訊N