作者：Mathieu Carrière
翻譯：孫韜淳
校對：和中華
本文約4500字，建議閱讀10分鐘
本文簡要介紹了機器學習中拓撲數據分析的力量並展示如何配合三個Python庫：Gudhi，Scikit-Learn和Tensorflow進行實踐。
標籤：數據可視化
Hi大家好。今天，我想強調下在機器學習中拓撲數據分析（TDA，Topological Data Analysis）的力量，並展示如何配合三個Python庫：Gudhi，Scikit-Learn和Tensorflow進行實踐。
拓撲數據分析？
首先，讓我們談談TDA。它是數據科學中相對小眾的一個領域，尤其是當與機器學習和深度學習對比的時候。但是它正迅速成長，並引起了數據科學家的注意。很多初創企業和公司正積極把這些技術整合進它們的工具箱中（比如IBM，Fujitsu，Ayasdi），原因則是近年來它在多種應用領域的成功，包括生物學、時間序列、金融、科學可視化、計算機圖形學等。未來我可能會寫一個關於TDA一般用途和最佳實踐的帖子，所以請大家等待下。
TDA：
https://en.wikipedia.org/wiki/Topological_data_analysis
IBM：
https://researcher.watson.ibm.com/researcher/view_group.php?id=6585
Fujitsu：
https://www.fujitsu.com/global/about/resources/news/press-releases/2016/0216-01.html
Ayasdi：
https://www.ayasdi.com/platform/technology/
生物學：
https://www.ncbi.nlm.nih.gov/pubmed/28459448
時間序列：
https://www.ams.org/journals/notices/201905/rnoti-p686.pdf
金融：
https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2931836
科學可視化：
https://topology-tool-kit.github.io/
計算機圖形學：
http://www.lix.polytechnique.fr/~maks/papers/top_opt_SGP18.pdf
TDA的目標是對你數據的拓撲性質進行計算和編碼，這意味著記錄數據集中多樣的連接成分，環，腔和高維結構。這非常有用，主要是因為其他描述符不可能計算這類信息。所以TDA真的儲存了一組你不可能在其他地方找到的數據特徵。現實情況是這類特徵已被證明對提升機器學習預測能力很有用，所以如果你以前還沒見過或聽過這類特徵，我來帶你快速瞭解一下。
我已經寫過很多這個主題的文章，你可以在Medium找到關於TDA的很多其他帖子，所以我不打算浪費時間在數學定義上面，而是通過解釋TDA文獻中的典型例子，來展示如何在你的數據集上應用TDA。
文章：
https://towardsdatascience.com/mixing-topology-and-deep-learning-with-perslay-2e60af69c321
帖子：
https://towardsdatascience.com/applied-topological-data-analysis-to-deep-learning-hands-on-arrhythmia-classification-48993d78f9e6
TDA的參考示例：點雲分類
這個數據集在一篇開創性的TDA文章上介紹過。它由通過下述動力系統生成的軌跡來得到的點雲集組成：
開創性的TDA文章
http://jmlr.org/papers/v18/16-337.html
一個動力系統的方程
這意味著我們將從一個單位正方形內隨機抽取一個初始點，並通過上面的方程生成一個點的序列。這將給我們一個點雲。現在我們可以根據意願重複這個操作，得到一堆點雲。這些點雲的一個有趣的屬性在於，根據你用來生成點序列的r參數的值，點雲會有非常不一樣且有意思的結構。比如，如果r=3.5，得到的點雲似乎覆蓋了整個單位正方形，但如果r=4.1，單位正方形的一些區域就是空的：換句話說，在你的點雲裡有好多洞。這對我們是個好消息：TDA可以直接計算這些結構是否能出現。
r=3.5（左）和r=4.1（右）計算出的點雲。相當明顯的是後者有個洞，但前者沒有
TDA跟蹤這些洞的方式實際上相當簡單。想象給定半徑為R的每個球的圓心都在你點雲的每個點上。如果R=0，這些球的並集就是點雲本身。如果R為無窮，那麼球的並集是整個單位正方形。但如果R被很精心的選擇，球的並集可能存在很多拓撲結構，比如，洞。
球並集的例子。對於中間圖的並集，它清晰的組成了一個洞。整張圖片被我“不要臉”地借用自我之前的一個帖子
帖子
https://towardsdatascience.com/a-concrete-application-of-topological-data-analysis-86b89aa27586
那麼，為了避免人工選擇R的“好值”，TDA將針對每一個可能的R值（從0到無窮）計算球的並集，並記錄每個洞出現或者消失時的半徑，並對一些點使用這些半徑值作為二維座標。TDA的輸出則是另一個點雲，其中每個點代表一個洞：這叫做Rips持續圖（Rips persistence diagram）。假設點雲在一個numpy數組X中儲存（shape為N*2），通過Gudhi，這個圖可以用兩行代碼計算出來：
<code>import gudhirips = gudhi.RipsComplex(points=X).create_simplex_tree()dgm = rips.persistence()/<code>
這個漂亮的持續圖由r=4.1對應的點雲計算出。紅色的點代表相連的成分，藍色的點代表洞
接下來我們將解決的任務則是給定點雲預測r的值。
通過Gudhi+Scikit-Learn進行拓撲機器學習
持續圖很簡潔，是不是？它們存在的問題則是，從不同點雲計算出的持續圖可能有不同數量的點（因為點雲可能有不同數量的洞）。所以如果你想用Scikit-Learn從持續圖中預測r，不幸的是，沒有直接的方法，因為這些庫預期輸入是一個結構化的向量。這也是為什麼目前大量的工作是關於將這些持續圖轉化為固定長度的歐幾里得向量，或者是開發對應的核。這很棒，但是你應該使用哪種呢？
不要擔心！Gudhi再一次給你解決辦法。通過它的表達（representation）模塊，你不僅可以計算所有的向量和核，甚至也可以使用Scikit-Learn來交叉驗證並且（或）選擇最佳的一種。就像下面這麼簡單:
表達
https://gudhi.inria.fr/python/latest/representations.html
<code>import gudhi.representations as tdafrom sklearn.pipeline import Pipelinefrom sklearn.svm import SVCfrom sklearn.ensemble import RandomForestClassifier as RFfrom sklearn.neighbors import KNeighborsClassifier as kNNfrom sklearn.model_selection import GridSearchCV
pipe = Pipeline([("TDA", tda.PersistenceImage()), ("Estimator", SVC())])param = [{"TDA": [tda.SlicedWassersteinKernel()], "TDA__bandwidth": [0.1, 1.0], "TDA__num_directions": [20], "Estimator": [SVC(kernel="precomputed")]}, {"TDA": [tda.PersistenceWeightedGaussianKernel()], "TDA__bandwidth": [0.1, 0.01], "TDA__weight": [lambda x: np.arctan(x[1]-x[0])], "Estimator": [SVC(kernel="precomputed")]}, {"TDA": [tda.PersistenceImage()], "TDA__resolution": [ [5,5], [6,6] ], "TDA__bandwidth": [0.01, 0.1, 1.0, 10.0], "Estimator": [SVC()]}, {"TDA": [tda.Landscape()], "TDA__resolution": [100], "Estimator": [RF()]}, {"TDA": [tda.BottleneckDistance()], "TDA__epsilon": [0.1], "Estimator: [kNN(metric="precomputed")]} ]model = GridSearchCV(pipe, param, cv=3)model = model.fit(diagrams, labels)/<code>
在前面的代碼中，我嘗試了帶切片Wasserstein核和持續權重Gaussian核的核SVM、帶有Persistence Images的C-SVM，帶有Persistence Landscapes的隨機森林，和一個帶有所謂的持久圖之間瓶頸距離（bottleneck distance）的簡單KNN。在Gudhi中還有許多其他的可能，所以你一定要試試！如果想了解更多細節你也可以看看Gudhi的Tutorial。
帶切片Wasserstein核：
http://proceedings.mlr.press/v70/carriere17a/carriere17a.pdf
持續權重Gaussian核：
http://proceedings.mlr.press/v48/kusano16.html
Persistence Images：
http://jmlr.org/papers/v18/16-337.html
Persistence Landscapes：
http://www.jmlr.org/papers/volume16/bubenik15a/bubenik15a.pdf
Gudhi的Tutorial：
https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-representations.ipynb
用Gudhi和Tensorflow/Pytorch進行拓撲優化
我很確信你目前已經成為了TDA的愛好者。如果你仍不相信，我還有其他的東西給你，這是受這篇論文啟發。想象你現在想解決一個更難的問題：我想讓你給我一個點雲，這個點雲的持續圖有儘可能多的點。換句話說，你需要生成一個有好多洞的點雲。
https://arxiv.org/abs/1905.12200
我可以看見你額頭上出汗了。但我是很仁慈的，轉眼間就能讓你知道Gudhi（1）可以做這個。想一想：當你生成一個持續圖時，這個圖中不同點的座標並不受全部的初始點雲影響，是不是？對於這個持續圖的一個給定點p，p的座標僅依賴於在初始點雲中組成p對應洞的點的位置，以一種簡單的方式：這些座標僅是球的並集使得這個洞出現或者消失時候的半徑；或者，等價表達是，這些點中的最大的成對距離。而Gudhi（2）可以通過它的persistence_pairs()函數找出這些關係。梯度則可以簡單的定義成歐幾里得距離函數的導數（正式定義見這篇論文）。
Gudhi（1）：
http://gudhi.gforge.inria.fr/python/latest/
Gudhi（2）：
https://gudhi.inria.fr/python/latest/
https://sites.google.com/view/hiraoka-lab-en/research/mathematical-research/continuation-of-point-cloud-data-via-persistence-diagram
接下來讓我們寫兩個函數，第一個從點雲中計算Rips持續圖，第二個計算持續圖點集的導數。為了可讀性我簡化了一點點代碼，實際的代碼可以從這裡找到。
<code>https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-optimization.ipynb
def compute_rips(x): rc = gd.RipsComplex(points=x) st = rc.create_simplex_tree() dgm = st.persistence() pairs = st.persistence_pairs() return [dgm, pairs]
def compute_rips_grad(grad_dgm, pairs, x): grad_x = np.zeros(x.shape, dtype=np.float32) for i in range(len(dgm)): [v0a, v0b] = pairs[i][0] [v1a, v1b] = pairs[i][1] grad_x[v0a,:]+=grad_dgm[i,0]*(x[v0a,:]-x[v0b,:])/val0 grad_x[v0b,:]+=grad_dgm[i,0]*(x[v0b,:]-x[v0a,:])/val0 grad_x[v1a,:]+=grad_dgm[i,1]*(x[v1a,:]-x[v1b,:])/val1 grad_x[v1b,:]+=grad_dgm[i,1]*(x[v1b,:]-x[v1a,:])/val1 return grad_x/<code>
現在讓我們把函數封裝進Tensorflow函數中（對Pytorch同樣簡單），並定義一個損失loss，這個損失是持續圖點到其對角線的距離的相反數。這將迫使圖有很多點，它們的縱座標比橫座標大得多。這樣的話，一個點雲會有很多大尺寸的洞。
<code>import tensorflow as tffrom tensorflow.python.framework import opsdef py_func(func, inp, Tout, stateful=True, name=None, grad=None): rnd_name = "PyFuncGrad" + str(np.random.randint(0, 1e+8)) tf.RegisterGradient(rnd_name)(grad) g = tf.get_default_graph() with g.gradient_override_map({"PyFunc": rnd_name}): return tf.py_func(func, inp, Tout, stateful=stateful, name=name)def Rips(card, hom_dim, x, Dx, max_length, name=None): with ops.op_scope([x], name, "Rips")as name: return py_func(compute_rips, [x], [tf.float32], name=name, grad=_RipsGrad)def _RipsGrad(op, grad_dgm): pairs = op.outputs[1] x = op.inputs[0] grad_x = tf.py_func(compute_rips_grad, [grad_dgm,pairs,x], [tf.float32])[0] return [None, None, grad_x, None, None]
tf.reset_default_graph()x = tf.get_variable("X", shape=[n_pts,2], initializer=tf.random_uniform_initializer(0.,1.), trainable=True)dgm, pairs = Rips(x)loss = -tf.reduce_sum(tf.square(dgm[:,1]-dgm[:,0]))opt = tf.train.GradientDescentOptimizer(learning_rate=0.1)train = opt.minimize(loss)/<code>現在我們開始優化！這是epochs 0,20,90的結果：
好多洞，好漂亮。。我們是不是在夢裡。如果你想往前看看，使用其它的損失，查閱這個Gudhi的tutorial。
https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-optimization.ipynb
最後的話
這個帖子僅是一瞥由Gudhi，Scikit-Learn和Tensorflow提供的眾多可能性。我希望我可以使你相信，在你的流程中整合TDA已經成為很簡單的事情。即使許多TDA應用已經在文獻中出現，肯定還有更多的應用需要去發現！
原文標題：
The Holy Trinity of Topological Machine Learning: Gudhi, Scikit-Learn and Tensorflow
原文鏈接：
https://towardsdatascience.com/the-holy-trinity-of-topological-machine-learning-gudhi-scikit-learn-and-tensorflow-pytorch-3cda2aa249b5
校對：林亦霖
譯者簡介
孫韜淳，首都師範大學大四在讀，主修遙感科學與技術。目前專注於基本知識的掌握和提升，期望在未來有機會探索數據科學在地學應用的眾多可能性。愛好之一為翻譯創作，在業餘時間加入到THU數據派平臺的翻譯志願者小組，希望能和大家一起交流分享，共同進步。
—完—
關注清華-青島數據科學研究院官方微信公眾平臺“ THU數據派 ”及姊妹號“ 數據派THU ”獲取更多講座福利及優質內容。
閱讀更多 數據派THU 的文章
