來源 | 數據派THU
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 gudhirips = gudhi.RipsComplex(points=X).create_simplex_treedgm = 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 importPipelinefrom sklearn.svm importSVCfrom sklearn.ensemble importRandomForestClassifier as RFfrom sklearn.neighbors importKNeighborsClassifier as kNNfrom sklearn.model_selection importGridSearchCV
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
defcompute_rips(x): rc = gd.RipsComplex(points=x) st = rc.create_simplex_tree dgm = st.persistence pairs = st.persistence_pairsreturn[dgm, pairs]
defcompute_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,:])/val1returngrad_x
現在讓我們把函數封裝進Tensorflow函數中（對Pytorch同樣簡單），並定義一個損失loss，這個損失是持續圖點到其對角線的距離的相反數。這將迫使圖有很多點，它們的縱坐標比橫坐標大得多。這樣的話，一個點雲會有很多大尺寸的洞。
import tensorflow as tffrom tensorflow.python.framework importopsdefpy_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_graphwithg.gradient_override_map({"PyFunc": rnd_name}):returntf.py_func(func, inp, Tout, stateful=stateful, name=name)defRips(card, hom_dim, x, Dx, max_length, name=None):withops.op_scope([x], name, "Rips")asname:returnpy_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_graphx = 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
【end】
曠視提Circle Loss，統一優化視角，革新深度特徵學習範式 | CVPR 2020清華學霸組團的工業AIoT創企再獲數千萬融資：玩家應推動在邊緣 AI 晶片上跑算法騰訊內測全新 Tim 3.0，支持微信登錄；滴滴順風車上線夜間服務；Angular 9.1發布為何你的 SaaS 想法總是失敗？沒想清楚這 4 個原因可能會繼續失敗！GitHub 疑遭中間人攻擊，無法訪問，最大暗網託管商再被黑！萬字好文：智能合約編寫之Solidity的編程攻略，建議收藏！你點的每個「在看」，我都認真當成了AI