其實是個常見的 vectorization trick，今天看一些 ML 的 code 又看見了一次，只是這次竟然沒有反應過來，想說就寫一篇 blog 當作加深自己印象也好。
其實這個 trick 的動機很單純: 給定一個 N x M的二維陣列，代表 N 個 M 維資料點，如何快速有效率地計算一個 N x N 二維陣列 D ? 其中 D[i, j] 是第 i 個點跟 j 個點的 Euclidean 距離。
一個直覺但比較沒效率的寫法就是直接用個 for-loop 去實作:
但如果我們利用 numpy 的 broadcasting，我們可以省去那個無效率的 for-loop:
一個簡單的測試:
有向量化的版本大概是 loop 的 100 倍快。
但眼尖的讀者可能會發現我並沒有去平方根。如果試著取平方根的話，你會發現 numpy 跳出了 RuntimeWarning:
把 dist_vec 印出來發現裡面竟然有負的值:
這也指出了一個向量化後的版本造成的數值上不穩定的問題，主要原因是因為我們現在是把很多個浮點數相加之後，再把這些和去做加減運算，基於浮點數本身精度的問題，可能就會產生這些無效的值，這是在實作時需要注意的點。那如果運用在實際 ML 的模型實作上，通常我就會建議把資料做一次 normalize。
不過最後我使用的小改動是 np.sqrt(np.maximum(sq + sq[:, None] + prod, 0.0)) 。效能上沒差太多，但就結果上跑個簡單測試:
大概可以有小數點下 6 位左右的精度，在一些狀況下或許夠用，但就差強人意，目前我還沒想到好的解法就是了。
感謝 Julia 社群的 Peter Cheng 提供了一個更好的寫法，雖然說他寫的是 Julia 但我硬轉成 numpy 了 XD
我補在我 gist 裡，另外一個叫 distance_matrix_trick_v2.py 的檔案。
補個跑的結果:
