獨家 | 拓撲機器學習的神聖三件套：Gudhi，Scikit-Learn和Tensorflow（附鏈接&代碼）

• 2020 年 3 月 26 日
• 筆記

作者：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，這個圖可以用兩行代碼計算出來：

import gudhi  rips = gudhi.RipsComplex(points=X).create_simplex_tree()  dgm = rips.persistence()

這個漂亮的持續圖由r=4.1對應的點雲計算出。紅色的點代表相連的成分，藍色的點代表洞

接下來我們將解決的任務則是給定點雲預測r的值。

通過Gudhi+Scikit-Learn進行拓撲機器學習

持續圖很簡潔，是不是？它們存在的問題則是，從不同點雲計算出的持續圖可能有不同數量的點（因為點雲可能有不同數量的洞）。所以如果你想用Scikit-Learn從持續圖中預測r，不幸的是，沒有直接的方法，因為這些庫預期輸入是一個結構化的向量。這也是為什麼目前大量的工作是關於將這些持續圖轉化為固定長度的歐幾里得向量，或者是開發對應的核。這很棒，但是你應該使用哪種呢？

不要擔心！Gudhi再一次給你解決辦法。通過它的表達（representation）模塊，你不僅可以計算所有的向量和核，甚至也可以使用Scikit-Learn來交叉驗證並且（或）選擇最佳的一種。就像下面這麼簡單:

表達

https://gudhi.inria.fr/python/latest/representations.html

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)

在前面的代碼中，我嘗試了帶切片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持續圖，第二個計算持續圖點集的導數。為了可讀性我簡化了一點點代碼，實際的代碼可以從這裡找到。

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

現在讓我們把函數封裝進Tensorflow函數中（對Pytorch同樣簡單），並定義一個損失loss，這個損失是持續圖點到其對角線的距離的相反數。這將迫使圖有很多點，它們的縱坐標比橫坐標大得多。這樣的話，一個點雲會有很多大尺寸的洞。

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)

現在我們開始優化！這是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