IA2 - $n$維球上的點

馬可夫過程 (Markov process) 生成三維球面上的點，此時點的位置不是均勻分佈，而是猶如醉漢走路 (drunkard's walk) 的情況。馬可夫過程在針對自然界運作或是社會結構演變的模擬中扮演著核心的角色！

一、二維圓盤的case

這次的專題非常簡單，我們先從二維的情況開始，二維球面上的點即為是分散在圓周上的點。我們知道如何產生散布在$(x,y)$上的點 (圖1)，但如何將之收縮到單位圓的圓周上呢？

圖1. 均勻散布在邊長為1的正方形內的藍點以及單位圓 (紅) 的範圍

　在圖1我們產生了1000個點，以及對應單位圓圈起來的面積。這次我們並不是要估計$\pi$，而是要找出哪些點分佈在紅色的圓周上，亦即產生出散布在單位圓上的點。你可能想說可以用以下的方式來實現 (以二維的case為例)：

1. 產生$N$的點的座標分別在$(-1,1)$之間
   $x={\tt rand}(-1,1)$
   $y={\tt rand}(-1,1)$
2. ${\tt if}~(x^2+y^2=1):~{\tt return}~(x,y)$
   ${\tt else:~return}~0$

上述的算法在第2步就是去篩出隨機生成的$N$個點中，哪些是符合半徑等於$1$的狀況 (即這些點在紅色的圓上)，如果不等於$1$則將之刪去。我不做計算，但我可以打包票，上面的算法接受率 (acceptance rate) 是超級低，生成數以萬計的點可能都不見得有一個是符合半徑等於$1$的條件。
　但換個角度來看生成的點，點座標$(x,y)$其實也代表著向量，既然單位圓有半徑為$1$，我們只要將每個點的向量都做歸一化的動作，我們就可以得到單位向量$(x^\prime,y^\prime)$ (即向量距原點的長度為$1$)，所以可以保證歸一化後的單位向量一定會在圖1的紅色圓上，而未歸一化前的點座標則用來決定在極座標上的角度$\phi$，即 $$\phi=\tan^{-1}\left(\frac{y}{x}\right).$$ 示意圖如圖2。這樣不僅保證所有的點都在單位圓上，又由於點是隨機生成的，在圓周上的位置$\phi$也是在$0$和$2\pi$之間均勻分佈。

圖2. 歸一化前與歸一化後的點，以及在極座標上的角度$\phi$

所以我們改寫上述的方式得到了算法1：

1. 產生$N$個點的座標分別在$(-1,1)$之間
   $x={\tt rand}(-1,1)$
   $y={\tt rand}(-1,1)$
2. 計算歸一化常數$r=\sqrt{x^2+y^2}$
3. ${\tt return}~(x/r,y/r)$

上述的算法可以用以下python的腳本實現：

def unitcircle(N):
    """
    Generating points on the unit circle.
    N: How many points we want to generate?
    """
    
    # Step 1
    X,Y = np.random.uniform(-1,1,size=(2,N))
    # Step 2: Normalization constant
    r = np.sqrt(X**2+Y**2)    
    return X/r,Y/r

而生成的點會像圖3一樣。透過算法1，我們隨機生成了100散佈在單位圓上的點，不像一開始用哪些點的半徑等於1的low acceptance rate的方式來取得，在算法1內每個點都可以保證在單位圓上！接下來我們進入三維的情形。

圖3. 透過算法1生成100在單位圓上的點

二、$n$維球面的case

原則上由算法1的方式，我們可以很簡變得拓展到三維球面的狀況，只需要在算法1的第1步多加上一個$z$軸的位置，並且讓第2步歸一化常數變成計算$r=\sqrt{x^2+y^2+z^2}$即可，結果如圖4。

圖4. 單位球球面上均勻分佈的紫色點點

　所以與其再寫一遍三維的算法，我們直接推廣到任意維的case，我們會發現一維 (點)、二維 (線) 和三維 (面) 都可以直接從這個算法導出。我們用$x_i$其中$i=1,2,3,\cdots,d$表示$d$維空間的第一個到第$d$個維的卡氏座標，所以算法2如下：

1. 生成$N$個$d$維的點座標
   $x_1={\tt rand}(-1,1)$
   $x_2={\tt rand}(-1,1)$
   $\qquad\quad\vdots$
   $x_d={\tt rand}(-1,1)$
2. 計算歸一化常數$r=\sqrt{\sum_{i=1}^d x_i^2}$
3. ${\tt return}~(x_1/r,x_2/r,\cdots,x_d/r)$

算法2我們可以用vector-notation的方式來實現：

def nsphere(N,dim=2):
    """
    Generating points on the unit sphere in arbitrary dimensions.
    N: How many points we want to generate?
    dim: Dimension of the sphere.
    """
    
    # Step 1
    X = np.random.uniform(-1,1,size=(dim,N))
    # Step 2: Normalization constant
    r = np.linalg.norm(X,axis=0)
    
    # The following returns:
    # 1st row: all N points' positions on x-axis
    # 2nd row: all N points' positions on y-axis
    # ...
    # d-th row: all N points' positions on d-axis
    # It can be projected to traditional (x,y,z,...,d) coordinate by transpose T
    return X/r

由上述python的腳本發現藉由numpy的便利性以及向量化的特性，使得生成任意維度的$N$個點都可以非常簡單的達成，算法1本身便包含在算法2內，但算法2的呈現比算法1更簡潔。其中我們預設維度參數dim是$2$，這就會生成圖3的情況，取$3$的話就是圖4，若dim$>3$其結果無法可視化 (原則上第4個維度可以作為時間，所以4維的情形並非完全不可視，而是需要以動畫的方式呈現，至於維度大於$4$我們無法直接看到高維球的全貌，只能看到它在$4$維以下的切面)，但電腦可以以邏輯的方式理解高維度的情形，人類也是如此，只不過看數字比看圖像來得令人頭痛😱。

三、小結

在這裡我們提供一個簡便生成任意維度單位球面上點的算法，而且也藉由python便利的性質，雖然我們用到許多駭人的術語，像是高維球之類的😩，但實際上操作起來並沒有很困難，此份章節我們也幾乎不用到數學，完全體現了第一原理模擬，從最簡單的概念出發，我們可以逼近自然界運作的方式，當然生成球面上的點和自然運作沒有什麼關聯就是😅。

　除了均勻分佈外，如本專題一開始透過馬可夫過程生成球面上的點，馬可夫過程表示每一個步驟所得到的結果會與它的上一步有關，但僅僅如此而已，與更早的過程完全無關，是無記憶性的 (memory-less)。但由於必須給定初始條件 (可以透過隨機生成初始條件)，我們看到點並不是均勻隨機分佈在球面上，而是會呈現一幅醉漢走路的狀態。馬可夫過程在許多模擬自然或社會運作的計算中扮演著核心的角色，有機會我們再回到這個議題上。


Source on github: I&A2_nSphere.ipynb