資料科學領域中，尤其是深度學習，很常用到多維度的資料結構。對於初入這個領域的人，理解這個就像是在做瑞文氏標準推理測驗一樣。
不過早在資料科學或是計量領域在處理多維度資料的時候，物理學跟數學界早就已經見怪不怪。這種「多維度」的概念被稱作張量 (tensor)。
話說這個詞被誤用得很徹底，不過對資訊界來說就當作是 N 維的資料。實際上張量要搭配線性轉換時的性質才能被好好定義，但反正這在資料科學中似乎不太重要。
在物理，尤其是廣義相對論中，因為多維度元素之間的相乘、相加實在太過頻繁，如果每一個都要寫出一個總和符號 Σ 還有上下標，數學式就會變得其醜無比。因此愛因斯坦就發展出了一個人類史上最重要的數學語法 — 愛因斯坦求和約定，而這也將與今天最後介紹的工具有關。
這一篇文章中，我會先舉一個多重維度的例子，以及傳統上用 numpy 的 broadcasting 來思考的話，會是多麼痛苦的事情。接著我會介紹同樣的問題，如何利用數學語言定義，與此同時透過 numpy/cupy/torch 裡面都有的 einsum 語法，來解決這個計算問題。
Youbike 流量與附近其他站點的關係
此例子來自我自己一個 project 中，特徵工程的一部分。
Youbike 是台灣的共享單車，光台北市就有將近 1300以上個站點，而台北市政府開放資料則有提供過去幾年的騎乘紀錄，包含起點、終點、與時間（去標籤，且解析度降為以小時為單位）。以此則可以研究影響 Youbike 流量的因子。
直觀上，Youbike 的流量不只與附近地理關係有關。有時候附近其他 Youbike 如果空了，也會來「求助於」遠一點的 Youbike。身為台大人這點再清楚不過了 — — 社科院晚上如果一車難求，真的會願意往 118 （台大附近的美食街）方向走一小段騎腳踏車。以我自己為例，我最多願意走 400 公尺找 youbike，再遠我寧可等有人騎過來（BTW，社科院到公館捷運站大概要走 1.2公里）。
因此，距離社科院 400 公尺內的任何其他 Youbike 站點的流量，無論是借出還是還車，都會受社科院這邊的站點的流量影響。在建立變數時，就會需要這兩個變數
- 附近 400 公尺內的總合借出量（排除自己）
- 附近 400 公尺內的總和還車輛（排除自己）
這兩個變數需要為每一個站點，每一個時間（小時為單位）建立。
直覺做法
我們先專心考慮一個地點、一個時段的情況。我要先找到在所有台北市的 Youbike 中，距離社科院 400 公尺內的站點有哪些。接著，我要看看在這個時段內，這些附近站點借出以及還車的數量有多少，全部加起來，就會是我需要的資料。聽起來十分簡單，但重點在於有很多車柱與時間。
為了之後方便解釋，假設現在考慮的社科院叫做 i，其他Youbike 站點叫做 j，時間是 t，而流量種類是 k（k=0 是借出、k=1 是還車）。
用數學的語言講，如果 W(i,j) 代表「 i 與 j 之間是否在 400 公尺以內」，y(i,t,k) 代表「i 在 t 時點的第 k 種流量」，那我要做的，其實就是算出
如果 i,j 在附近，則 W(i,j) = 1，那麼 j 站點當時的流量 y(j,t,k) 就會被加上去。但如果 i,j 很遠，則 W(i,j) = 0，那麼 y(j,t,k) 就不會被納入加總。
如果要利用 for 來算出每一個站點每一個時點的 s(i,t,k)的話，大概會這這個步驟：
- 對於其中一個時間[t]
- 對於其中一個Youbike 站點 [i]，比方說「社科院」
- 對於除了社科院以外的所有其他站點[j]，查看他的距離。如果小於 400，再下一步
- 加總鄰近站點的借出量[k=0] 與 還車量[k=1]
- 重複 3. -> 2. -> 1.，跑遍所有 j, i, t
不過，如果以小時計算，光一年就有8,760小時，搭配 1300 個Youbike 站點總共 C(1300,2) * 8760 = 7396506000 次計算。不要忘記借出還車都要計算一次，所以光一年其實有14,793,012,000 個計算或是資料擷取要處理。
這時候就算做法直覺，最好還是依賴一下有優化過的做法，因此我們選擇用既有的數值計算套件 Numpy 來處理比較好～
Numpy
用 Numpy 的話，最好就不要再分維度進行 for 迴圈，而要「一整坨一起處理」。這個操作叫 Vectorization。 由於 Numpy 會針對一整坨東西的操作進行最佳化處理，所以效率上會快很多。
比方說，如果一個向量每一個元素都要 +1，課本上的做法就是
old_array = np.array([1,2,3,4,5])
new_array = old_array + 1 # 雖然一個是向量一個是純量，但 Numpy 會自動處理這種把一整個向量加上一個純量的概念叫 Broadcasting，這使得Numpy 可以自動化「一整坨」進行操作，也就是 Vectorization。
不過這個最大的問題是，上面這個 Youbike 的問題到底要怎麼一整坨處理呢？
處理的做法就是手動「擴增維度」。從圖像來看的話，或許最容易理解。但如果太抽象，這部分可以跳過。
先來看看總流量的資料。我們一共有 1300 個Youbike 站點，姑且先簡化為 N，也就是 i 與 j 都是從 0 到 N-1。資料從 2020一路到 2023一共有 T 個小時。而每一個站點在每一個時間點，我們都有進/出的數量 （k=0 與 k=1）。換句話說，這是一個 N x T x 2 的資料，我們叫他 y：
為了後續方便起見，我將一個格子裡面塞了兩個數字，分別代表借出與還車數量。
至於站點跟站點之間，暫時不管他會不會移動位置，或是有刪減站點，一律以 2023 年的經緯度資料為主的話，我們可以算出「是否在 400公尺內」的一個矩陣。假如社科院[i] 跟法學院[j] 之間距離在 400公尺內，那麼在(i, j) 這個「座標」上面，就標記一個 1，反之則為 0。以圖論的角度來看，這相當於把距離 400內的兩個站點畫上一條線，而這個 N x N 矩陣就是這個圖的鄰接矩陣(Adjacency matrix)，我們叫他 W：
當然，因為要看鄰居，所以自己跟自己 (i,i) 之間就不能考慮。
我們接著要來讓兩個資料可以很好的「融合」在一起。對於進出流量的資料(i,t,k) 來說，我們想考慮當下每一個站點跟每一個站點之間的組合(i,j)。對於站點跟站點之間的資料(i,j)，我們則想考慮每一個時間點其中一方站點的進出資料(j,t,k)。也就是說，對於兩筆資料，我們都想要「往外延伸」去考慮另一個維度的事情，也就是我們需要先某種程度上「延伸一個維度」：
再把資料，根據這個維度重複過去：
左邊的盒子，是將「每個站點每個時刻的進與出」，沿著第二個維度複製出來 N 個重複的。
右邊的盒子，是將「每個站點與站點之間靠近與否」，沿著第三個（與第四個）維度複製出來。
這兩個盒子之間在延伸了之後，維度都是 N x N x T x 2，是個四階的張量。如果把左邊一筆資料，乘上右邊相對應位置的資料，就代表著：第 i 個站點，在 t 這個時間點下，對另外一個站點[j]在某流量[k]上的貢獻。
舉例來說，如果在右邊這個盒子中， 第 i 個站點跟第 j 個站點之間的關係不存在 (等於0），則相乘後，第 i 個站點，在 t 這個時間點下，來自第 j 個站點任何種類的流量（進或是出），完全不會對第 i 個站點產生貢獻。
最後再沿著第二維相加就可以，因為第二維代表的是與 i 不同的其他站點，所以加總之後就可以得到 i 在這時間點的附近流量總和。
如果熟悉矩陣的話，這有一點類似矩陣乘法，只不過因為想要批次處理，所以維度變多了。實際上在做的是固定時間與流量後，針對中間 2,3 維的矩陣相乘 。
程式碼
首先，我們將一個三維資料 y，沿著第二維（程式的語言中是第 1 = 2–1 維，因為是從 0 開始數）增加一個維度，其他依序。我們又將一個二維的距離資料 W，沿著三四維，增加兩個維度：
extended_y = y[:, np.newaxis, ...]
extended_W = W[:, :, np.newaxis, np.newaxis]利用 np.newaxis第一個變數 extended_y 現在形狀從 (N, T, 2)變成 (N, 𝟙, T, 2)，而 extended_W 則是 (N, N)變成(N, N, 𝟙, 𝟙)
雖然形狀不一樣，但是如果兩兩配對下，數字一樣或是其中一個是 1，在 numpy 中都是可以自動「延伸」的（詳見 Broadcasting）。因此相乘時
extended_y ： (N, 𝟙, T, 2) -> (N, N, T, 2)
extended_W ：(N, N, 𝟙, 𝟙) ->(N, N, T, 2)
換句話說，相乘時
extended_y * extended_W就會是上面，元素跟元素之間相乘的結果。如果要看對單一站點，在各時點的周維流量加總，則需要對第二個維度（還記得我們在第二個維度上，延伸出了其他站點，所以現在需要加總起來）加總
s = np.sum(extended_y * extended_W , axis = 1)這樣，就完成了我們所想要的
各站點、各時點、附近 400 公尺內的總合借出/還車量
這時他還是一個(N x 1 x T x 2)的資料，不過利用 reshape適當重整形狀之後，他會是一個 N x T x 2 的特徵，可以放回 pandas 的資料中。
而每一個元素則會是
有沒有更簡單的做法？
上面的作法，是非常進階的 broadcasting 技巧，不只要對於Numpy的語法非常熟悉，也要很清楚的判斷各階段操作時的資料形狀。幸運的時候，出錯會出現錯誤訊息，而不幸的話，可能因為不清楚 broadcasting 到底背後自動重複了哪些維度，而出來一個似是而非卻沒錯誤訊息的結果，且這輩子別想 debug。
另外，延伸資料維度的做法非常耗記憶體，雖然如果成功的話速度極快，但以上面的例子，會需要用到 101GB 的記憶體空間，根本跑不動。
這時候，一個鮮少被使用的函數 np.einsum() 就可以派出用場了。
愛因斯坦求和記號
還記得我們要求的：
用愛因斯坦求和記號，可以直接簡化成
沒錯，就只是把加總符號Σ 省略掉。不過觀察等號左邊，可以發現沒有 j ，但右邊兩項都有 j，表示一定是被加總加掉了！
因此，這一項張量間的操作，可以完全由座標代號之間的關係來定義：
ij, jtk -> itk
而在 numpy/torch/tensorflow 中，都有一個可以直接利用上面這個神秘字串的函式，那就是 einsum 。我們完全不用自己延伸維度或是加總某一方向，只需要：
s = np.einsum("ij, jtk -> itk", W, y)就可以了！而如果資料很大，可以再加入一個 optimize = True 進行最佳化。
他不只可以將平常需要 broadcast 到不同維度的過程省略，還可以自動計算節省記憶體的做法，同時也兼顧效率。
更多 einsum 的用法，大力推薦這篇：https://zhuanlan.zhihu.com/p/27739282
另外，Numpy 是純 CPU 計算，如果想要利用 GPU 來操作上面的話，也可以用 cupy 或是專門進行深度網路的 torch / tensorflow ，都有提供 einsum 這個函式。經過實驗，利用 Numpy 執行上面的操作，平均每次需要 47.5秒，而如果透過 GPU（利用 cupy），則只需要29.9毫秒，也就是1600倍的加速！而相比原本用 for 迴圈跑的 4分鐘以上的話，則保底 8000倍以上的加速！
結論
當資料量過大時，用太多層 for 會非常緩慢，甚至佔記憶體。而傳統上進階的教學則是建議搭配 broadcasting 變成高維度之後，透過 vectorization 來實現「單指令流多資料流 (SIMD)」，藉此加速運算。然而碰到資料量太大的時候，記憶體隨便都會不足。
einsum 這個函數則不但可以跳過複雜的腦內激盪，直接利用數學表達式，更有直接寫好的最佳化演算法可以計算，也確保不會記憶體不足（應該）。而如果能透過 cupy 或是 torch 等等套件將陣列移轉至 GPU上進行平行運算，則這原本要將近14,793,012,000次操作的工作，則只需要不到 30 毫秒。
