淺談獨立特徵（independent features）、潛在特徵（underlying features）提取、以及它們在網路安全中的應用
- 2019 年 10 月 3 日
- 筆記
1. 關於特徵提取
0x1：什麼是特徵提取
特徵提取研究的主要問題是，如何在數據集未明確表示結果的前提下，從中提取出重要的潛在特徵來。和無監督聚類一樣，特徵提取演算法的目的不是為了預測，而是要嘗試對數據進行特徵識別，以此得到隱藏在數據背後的深層次意義。
回想一下聚類演算法的基本概念，聚類演算法將數據集中的每一行數據分別分配給了某個組（group）或某個點（point），每一項數據都精確對應於一個組，這個組代表了組內成員的平均水平。
特徵提取是這種聚類思想更為一般的表現形式，它會嘗試從數據集中尋找新的數據行，將這些新找到的數據行加以組合，就可以重新構造出數據集。和原始數據集不一樣，位於新數據集中的每一行數據並不屬於某個聚類，而是由若干特徵的組合構造而成的。
從特徵提取的這個原理來看，特徵提取也可以理解為一種泛化的降維概念。
在這篇文章中，筆者會嘗試從底層的數學原理出發，闡述這些概念之間的聯繫和區別。其實不論是特徵提取、降維、聚類，都只是從不用的角度解決同一類問題，其背後的原理都是共通的。
0x2：獨立特徵提取的典型應用場景
1. 雞尾酒宴會問題
這是一個如何在多人談話時鑒別聲音。人類聽覺系統的一個顯著特徵就是：在一個人聲鼎沸屋子裡，儘管有許多不同的聲音混雜在一起湧入我們的耳朵，可我們還是能夠從中鑒別出某個聲音來，大腦非常擅長於從聽到的所有雜訊中分離出單獨的聲音。
同樣的目標，通過本文將要討論的獨立特徵提取技術，電腦就有可能完成同樣的任務。
2. 新聞主題分類
獨立特徵提取的一個重要應用是，對重複出現於一組文檔中的單詞使用組合進行模式識別（word-usage pattern recognition），這可以幫助我們有效地識別出，以不同組合形式出現於各個文檔中的主題。從文檔中提取出的這些主題，就可以作為一個泛化的一般化特徵，用於概括一整類文旦。同時，通過對文檔進行獨立特徵識別，我們會發現，一篇文章可以包含不止一個主題，同時，同一個主題也可以用於不止一篇文章。
例如，對於某個重大新聞，往往會有多家報社從不同角度都進行了報道，雖然各家報社的側重點各有不同，但是不可避免會有一些公共主題是交叉的。具體來說，例如美國大選川普當選總統，CNN和紐約時報都進行了報道，兩家報紙的側重點各有不同，但是都不約而同地談到了川普過去的從商經歷等。
需要特別注意的是，必須要保證要搜索的文檔中存在內容重疊的主題（acrossover topic），如果沒有任何文檔存在共同之處，那麼演算法將很難從中提取出重要的通用特徵，如果強行提取則最終得到的特徵將會變得毫無意義（這也體現了獨立特徵提取技術的降維和聚類本質）。
3. 股票市場數據分析
股票市場是一個集體智慧共同作用的典型場景，在股市經濟活動中，每一個自然人都遵循最大收益原則開展活動（凱恩斯的開不見的手宏觀經濟理論），通俗地說就是，每個投資者做每項投資的目的都是為了使自己的收益最大化。整個社會千千萬萬個投資者共同形成了一個股票市場。
基於這個前提下，我們假設股票市場數據背後潛藏著諸多原因，正是這些原因共同組成的結果，導致了證券市場的結果。我們可以將獨立特徵提取演算法應用於這些數據，尋找數據背後的原因，以及它們各自對結果所構成的影響。
獨立特徵提取技術可以用來對股市的成交量進行分析。所謂成交量，就是指在某一給定時間段內所買賣的股票份數，下圖是Yahoo!股票的走勢圖，
位於最上方的線代表了收盤價，下面的柱狀圖則給出了成交量。
我們發現，當股票價格有較大變化時，成交量在那幾天往往就會變得很高。這通常會發生在公司發表重要聲明或發布財務報告的時候。此外，當有涉及公司或業界新聞報道時，也會導致價格出現”突變“。在缺少外部事件影響的情況下，對於某隻股票而言，成交量通常（但不總是）保持不變的。
從成交量中提取出的獨立特徵，基本上就是明面上或者背後的某些”利好“或”不好“事件。
Relevant Link:
《集體智慧編程》Toby segaran - 第10章
2. 非負矩陣因式分解（non-negative matrix factorization，NMF）
0x1：以文章矩陣為例簡要說明NMF
NMF是一個非常數學化的概念，我們本章會詳細討論其底層數學原理，不過，筆者希望通過一個簡單易懂的例子來為讀者朋友先建立一個感性的直觀概念，為之後的原理討論作鋪墊。
假設我們現在通過詞頻語言模型，已經將一個文檔集（很多篇文章）轉換為了一個【M，N】矩陣，M代表有M篇文章，N代表有N個詞，元素的數值代表該單詞在對應文章中的詞頻計數。
我們的目標是要對著個矩陣進行因式分解，即，找到兩個更小的矩陣，使得二者相乘以得到原來的矩陣。這兩個矩陣分別是特徵矩陣和權重矩陣。
【特徵矩陣】
在該矩陣中，每個特徵對應一行，每個單詞對應一列。矩陣中的數字代表了某個單詞相對於某個特徵的重要程度。
由於每個特徵都應該對應於在一組文章中出現的某個主題，因此假如有一篇文章報道了一個新的電視秀節目，那麼也許我們會期望這篇文章相對於單詞”television“能夠有一個較高的權重值。
一個特徵矩陣示例
由於矩陣中的每一行都對應於一個由若干單詞聯合組成的特徵，因此很顯然，只要將不同數量的特徵矩陣按行組合起來，就有可能重新構造出文章矩陣來。
而如何組合不同的特徵矩陣，用多少特徵矩陣來組合，就是由權重矩陣來控制。
【權重矩陣】
該矩陣的作用是，將特徵映射到文章矩陣，其中每一行對應於一篇文章每一列對應於一個特徵。矩陣中的數字代表了每個特徵應用於每篇文章的程度。
一個權重矩陣示例
下圖給出了一個文章矩陣的重構過程，只要將權重矩陣與特徵矩陣相乘，就可以重構出文章矩陣，
遵照矩陣乘法法則，特徵矩陣的行數和權重矩陣的列數必須要相等。如果特徵數量與文章數量恰好相等，那麼最理想的結果就是能夠為每一篇文章都找到一個與之完美匹配的特徵。
在獨立特徵提取領域中，使用矩陣因式分解的面對，是為了縮減觀測數據（例如文章）的集合規模，並且保證縮減之後足以反映某些共性特徵。理想情況下，這個較小的特徵集能夠與不同的權重值相結合，從而完美地重新構造出原始的數據集。但實際情況中，這種可能性是非常小的，因此演算法的目標是要儘可能地重新構造出原始數據集來。
筆者思考：
筆者這裡帶領大家用形象化的思維來理解一下矩陣因式分解，我們將其想像為我們吃月餅，從中間將其掰開為兩半，一半是特徵矩陣，一半是權重矩陣。特徵矩陣和權重矩陣原本都不存在，因為我們一掰，憑空出現了2個小的月餅。那接下來的問題來了，我們能否隨意的掰這個月餅呢？答案是不行！這個月餅有自己的法則，只允許我們按照有限幾種方案來掰，因為該法則要求掰開後的兩半還必須能夠完美的拼回一個完整的月餅。
回到矩陣因式分解上來，矩陣的因式分解類似於大數的質因分解，一個矩陣只存在少量幾種因式分解方法。而要找到這幾種分解方案，就需要使用一些非常精巧的演算法，例如本文要介紹的乘法更新法則（multiplicative update rules）。
0x2：乘法更新法則（multiplicative update rules）- 尋找矩陣因式分解的一種演算法
我們的目標很明確，希望找到兩個矩陣（滿足M和N的約束即可），這兩個矩陣相乘要儘可能接近原始目標矩陣，我們也建立了損失函數difcost，通過計算最佳的特徵矩陣和權重矩陣，演算法嘗試盡最大可能來重新構造文章矩陣。我們通過difcost函數來度量最終的結果與理想結果的接近程度。
一種可行的優化方法是將其看作是一個優化問題，藉助模擬退火或者遺傳演算法搜索到一個滿意的題解，但是這麼做的搜索成本可能過於龐大（隨著原始矩陣尺寸的上升），一個更為有效的方法是，使用乘法更新法則（multiplicative update rules）。
這些法則產生了4個更新矩陣（update matrices），這裡我們將最初的文章矩陣稱為數據矩陣。
- hn：經轉置後的權重矩陣與數據矩陣相乘得到的矩陣
- hd：經轉置後的權重矩陣與原權重矩陣相乘，再與特徵矩陣相乘得到的矩陣
- wn：數據矩陣與經轉置後的特徵矩陣相乘得到的矩陣
- wd：權重矩陣與特徵矩陣相乘，再與經轉置後的特徵矩陣相乘得到的矩陣
為了更新特徵矩陣和權重矩陣，演算法需要做如下幾個操作，
- 首先將上述所有矩陣都轉換成數組
- 然後將特徵矩陣中的每一個值域hn中的對應值相乘，併除以hd中的對應值
- 再將權重矩陣中的每一個值域wn中的對應值相乘，併除以wd中的對應值
from numpy import * import numpy as np def difcost(a,b): dif=0 for i in range(shape(a)[0]): for j in range(shape(a)[1]): # Euclidean Distance dif+=pow(a[i,j]-b[i,j],2) return dif def factorize(v,pc=10,iter=50): ic=shape(v)[0] fc=shape(v)[1] # Initialize the weight and feature matrices with random values w=matrix([[random.random() for j in range(pc)] for i in range(ic)]) h=matrix([[random.random() for i in range(fc)] for i in range(pc)]) # Perform operation a maximum of iter times for i in range(iter): wh=w*h # Calculate the current difference cost=difcost(v,wh) if i%10==0: print cost # Terminate if the matrix has been fully factorized if cost==0: break # Update feature matrix hn=(transpose(w)*v) hd=(transpose(w)*w*h) h=matrix(array(h)*array(hn)/array(hd)) # Update weights matrix wn=(v*transpose(h)) wd=(w*h*transpose(h)) w=matrix(array(w)*array(wn)/array(wd)) return w,h if __name__ == '__main__': l1 = [[1,2,3], [4,5,6]] m1 = matrix(l1) m2 = matrix([[1,2], [3,4], [5,6]]) print "np.shape(m1*m2): ", np.shape(m1*m2) w, h = factorize(m1*m2, pc=3, iter=100) print "w: ", w print "h: ", h print "w*h: ", w*h print "m1*m2: ", m1*m2
可以看到，演算法成功地找到了權重矩陣和特徵矩陣，使得這兩個矩陣相乘的結果與原始矩陣幾乎完美匹配。
值得注意的是，MUR方法要求我們明確指定希望找到的特徵數。例如，在一段錄音中的兩個聲音，或者是當天的5大新聞主題。
0x3：理論概述
1. 演算法形式化描述
有了前兩節的鋪墊，現在我們來討論一些NMF的理論概念。
NMF（Non-negative matrix factorization），即對於任意給定的一個非負矩陣V，其能夠尋找到一個非負矩陣W和一個非負矩陣H，滿足條件V=W*H。從而將一個非負的矩陣分解為左右兩個非負矩陣的乘積。其中，
- V矩陣中每一列代表一個觀測(observation)，每一行代表一個特徵(feature)
- W矩陣稱為基矩陣，即特徵矩陣
- H矩陣稱為係數矩陣或權重矩陣，這時用係數矩陣H代替原始矩陣，就可以實現對原始矩陣進行降維，得到數據特徵的降維矩陣，從而減少存儲空間（從【M*N】降到【M*K】維）
分解過程如下圖所示，
2. 損失函數形式
對於如何衡量因式分解的效果，有3種損失函數形式，其中第一種我們上前面的例子中已經使用到了簡化版本。
【squared frobenius norm】
其中：
α為L1&L2正則化參數，?為L1正則化佔總正則化項的比例，||*||_1為L1範數。
我們前面章節中的difcost函數，就是當α和?都為0的的時候的一種損失函數特例。
【Kullback-Leibler (KL)】
【Itakura-Saito (IS)】
實際上，上面三個公式是beta-divergence family中的三個特殊情況（分別是當β=2,1,0），其原型是：
0x4：使用NMF對鳶尾花數據進行獨立特徵提取
這裡我們不使用手寫函數了，直接用sklearn封裝好的NMF庫。
# -*- coding: utf-8 -*- from sklearn.decomposition import NMF from sklearn.datasets import load_iris X, _ = load_iris(True) # can be used for example for dimensionality reduction, source separation or topic extraction nmf = NMF(n_components=4, # k value,默認會保留全部特徵 init=None, # W H 的初始化方法，包括'random' | 'nndsvd'(默認) | 'nndsvda' | 'nndsvdar' | 'custom'. solver='cd', # 'cd' | 'mu' beta_loss='frobenius', # {'frobenius', 'kullback-leibler', 'itakura-saito'}，一般默認就好 tol=1e-4, # 停止迭代的極限條件 max_iter=200, # 最大迭代次數 random_state=None, alpha=0., # 正則化參數 l1_ratio=0., # 正則化參數 verbose=0, # 冗長模式 shuffle=False # 針對"cd solver" ) # -----------------函數------------------------ print('params:', nmf.get_params()) # 獲取構造函數參數的值，也可以nmf.attr得到，所以下面我會省略這些屬性 nmf.fit(X) W = nmf.fit_transform(X) W = nmf.transform(X) nmf.inverse_transform(W) # -----------------屬性------------------------ H = nmf.components_ # H矩陣 print('reconstruction_err_', nmf.reconstruction_err_) # 損失函數值 print('n_iter_', nmf.n_iter_) # 實際迭代次數 # -----------------輸出------------------------ # H即為降維後的新矩陣（維度4是我們通過參數指定的） print "H: ", H print "W: ", W
再看一個更具體的影像處理的例子，使用NMF對影像進行降維，提取出影像中的關鍵元素，
已知Olivetti人臉數據共400個，每個數據是64*64大小，原始數據維度4096。
NMF分解得到的W矩陣相當於從原始矩陣中提取的特徵，那麼就可以使用NMF對400個人臉數據進行特徵提取。
通過設置k的大小，設置提取的特徵的數目。在本實驗中設置k=6，隨後將提取的特徵以影像的形式展示出來。
# -*- coding: utf-8 -*- from time import time from numpy.random import RandomState import matplotlib.pyplot as plt from sklearn.datasets import fetch_olivetti_faces from sklearn import decomposition n_row, n_col = 2, 3 n_components = n_row * n_col image_shape = (64, 64) rng = RandomState(0) # ############################################################################# # Load faces data dataset = fetch_olivetti_faces('./', True, random_state=rng) faces = dataset.data n_samples, n_features = faces.shape print("Dataset consists of %d faces, features is %s" % (n_samples, n_features)) def plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray): plt.figure(figsize=(2. * n_col, 2.26 * n_row)) plt.suptitle(title, size=16) for i, comp in enumerate(images): plt.subplot(n_row, n_col, i + 1) vmax = max(comp.max(), -comp.min()) plt.imshow(comp.reshape(image_shape), cmap=cmap, interpolation='nearest', vmin=-vmax, vmax=vmax) plt.xticks(()) plt.yticks(()) plt.subplots_adjust(0.01, 0.05, 0.99, 0.93, 0.04, 0.) # ############################################################################# estimators = [ ('Non-negative components - NMF', decomposition.NMF(n_components=n_components, init='nndsvda', tol=5e-3)) ] # ############################################################################# # Plot a sample of the input data plot_gallery("First centered Olivetti faces", faces[:n_components]) # ############################################################################# # Do the estimation and plot it for name, estimator in estimators: print("Extracting the top %d %s..." % (n_components, name)) t0 = time() data = faces estimator.fit(data) train_time = (time() - t0) print("done in %0.3fs" % train_time) components_ = estimator.components_ print('components_:', components_.shape, 'n**n', components_) plot_gallery('%s - Train time %.1fs' % (name, train_time), components_) plt.show()
原始400影像中的前6張（普通人臉）
NMF從原始400張影像中提取出的6個獨立通用特徵
這個例子非常形象直觀地說明了NMF提取出的特徵矩陣的含義。對於人臉影像來說，輸入NMF演算法的是代表了原始影像的像素矩陣，而矩陣因式分解得到的權重矩陣W則代表了一種像素權重重新分配分配方案。
具體地說，上述6張特徵向量，每個向量的維度都等於原始影像（4096）維，但區別在於，每個特徵中，各個像素的值是不同的，它們分別代表了不同的特徵。有的特徵側重於對鼻子部分的特徵提取，所以在鼻子區域的像素權重分配的就會多，其他的以此類推。
0x5：能否用模擬退火演算法或遺傳演算法這類隨機優化演算法得到和NMF因式分解近似的效果呢？
這個小節，我們來看一個有趣的想法，能夠用隨機優化技術直接得到矩陣的因式分解結果呢？我們來寫一個實驗程式碼看看結果：
# -*- coding: utf-8 -*- from numpy import * import numpy as np def difcost(a, b): dif = 0 for i in range(shape(a)[0]): for j in range(shape(a)[1]): # Euclidean Distance dif += pow(a[i, j] - b[i, j], 2) return dif def factorize(v, pc=10, iter=50): ic = shape(v)[0] fc = shape(v)[1] # Initialize the weight and feature matrices with random values w = matrix([[random.random() for j in range(pc)] for i in range(ic)]) h = matrix([[random.random() for i in range(fc)] for i in range(pc)]) # Perform operation a maximum of iter times for i in range(iter): wh = w * h # Calculate the current difference cost = difcost(v, wh) #if i % 10 == 0: # print cost # Terminate if the matrix has been fully factorized if cost == 0: break # Update feature matrix hn = (transpose(w) * v) hd = (transpose(w) * w * h) h = matrix(array(h) * array(hn) / array(hd)) # Update weights matrix wn = (v * transpose(h)) wd = (w * h * transpose(h)) w = matrix(array(w) * array(wn) / array(wd)) return w, h def difcost_2d(a, b): dif = 0 for i in range(shape(a)[0]): for j in range(shape(a)[1]): # Euclidean Distance dif += pow(a[i, j]-b[i, j], 2) return dif def hillclimb_2d(v, pc=2, costf=difcost_2d): ic = shape(v)[0] fc = shape(v)[1] # Initialize the weight and feature matrices with random values w = np.array([[random.random() for j in range(pc)] for i in range(ic)]) h = np.array([[random.random() for i in range(fc)] for i in range(pc)]) # Create a random solution sol_w = [] sol_h = [] for x_axi in v: sol_w.append([random.randint(0, y_axi) for y_axi in x_axi]) sol_h.append([random.randint(0, y_axi) for y_axi in x_axi]) # Main loop for w min_cost, current_cost, last_round_cost = 9999999999, 0, 0 patient_count = 0 while 1: best_tmp_w = [] best_tmp_h = [] for j in range(len(w)): # an tmp copy of x for w tmp_w = w # One away in each direction if sol_w[j][0] > w[j][0]: tmp_w[j][0] = sol_w[j][0] else: tmp_w[j][0] = w[j][0] if sol_w[j][1] > w[j][1]: tmp_w[j][1] = sol_w[j][1] else: tmp_w[j][1] = w[j][1] # record the candidate best_tmp_w.append(tmp_w) # now update h matrix for i in range(len(h)): # an tmp copy of x for h tmp_h = h # One away in each direction if sol_h[i][0] > h[i][0]: tmp_h[i][0] = sol_h[i][0] else: tmp_h[i][0] = h[i][0] if sol_h[i][1] > h[i][1]: tmp_h[i][1] = sol_h[i][1] else: tmp_h[i][1] = h[i][1] # record the candidate best_tmp_h.append(tmp_h) # hill climb one dimension one time # See what the best solution amongst the neighbors is i_w, j_h = 0, 0 for i in range(len(best_tmp_w)): for j in range(len(best_tmp_h)): current_cost = costf(v, best_tmp_w[i] * best_tmp_h[j]) if current_cost <= min_cost: i_w = i j_h = j min_cost = current_cost # update w, h w = best_tmp_w[i_w] h = best_tmp_h[j_h] # If there's no improvement, then we've reached the top if min_cost == last_round_cost and patient_count >= 99: print "min_cost == last_round_cost" break else: last_round_cost = min_cost patient_count += 1 print "patient_count: ", patient_count print "min_cost: ", min_cost print "last_round_cost: ", last_round_cost #print "w: ", w #print "h: ", h return w, h if __name__ == '__main__': m1 = np.array([[1, 2, 3], [4, 5, 6]]) m2 = np.array([[1, 2], [3, 4], [5, 6]]) v = np.dot(m1, m2) print "np.shape(v): ", np.shape(v) print "v: ", v print " " w_optimization, h_optimization = hillclimb_2d(v=v, pc=2, costf=difcost_2d) print "w_optimization: ", w_optimization print "h_optimization: ", h_optimization print " " print "w_optimization * h_optimization: ", w_optimization * h_optimization
從結果可以看到，隨機優化的結果並不好，之所以號稱”萬能優化演算法“的隨機優化演算法，在這個場景下不管用了，主要原因如下：
- 隨機優化技術（爬山法、模擬退火、遺傳演算法、SGD）本質上都是一種梯度導向（或近似梯度導向）的迭代式優化演算法，它們要求待優化的損失函數具有連續可微的特性。通俗地說就是，每次僅僅輕微修改某一個維度中的少量值，損失函數就能獲得一定程度的輕微優化。但是矩陣的因式分解這個任務對應的損失函數可能是完全不規則的損失曲線，輕微的修改帶來的損失值的變更是無法預料的，可能很小，也可能突然非常大，這非常不利於損失函數進行迭代優化
- 矩陣的因式分解對應的損失函數可能是一個非常不規則的損失曲線，因此隨機優化演算法很容易陷入局部最優值而無法跳出
反過來也看到，那什麼場景適合用隨機優化演算法呢？
- 影像檢測與識別：影像的像素之間是漸變的，模型的權重調整也是連續漸變的
- 文本檢測與分類：以惡意郵件檢測為例，要判斷一篇文檔是否是惡意的，和該文檔中出現的詞（word）以及詞頻（word frequency）有關，好的文檔和壞的文檔之間的差距往往就是某些詞是否出現，以及某些詞出現多少次這些因素決定，因此，損失函數也是連續可微的，非常適合用隨機優化演算法來得到近似全局最優解
Relevant Link:
https://en.wikipedia.org/wiki/Non-negative_matrix_factorization https://blog.csdn.net/jeffery0207/article/details/84348117 http://latex.91maths.com/ 《Learning the parts of objects by non-negative matrix factorization》D.D.Lee，H.S.Seung
3. PCA（Principal Component Analysis）
廣義上來說，PCA是一種基於線性降維的特徵提取方法。
我們還是用人臉影像來闡述這個概念，我們的目的是利用PCA演算法從一堆人臉影像數據中，找到一組共同特徵，即所謂的大眾臉特徵，
# -*- coding: utf-8 -*- import logging from time import time from numpy.random import RandomState import matplotlib.pyplot as plt from sklearn.datasets import fetch_olivetti_faces from sklearn.cluster import MiniBatchKMeans from sklearn import decomposition # Display progress logs on stdout logging.basicConfig(level=logging.INFO, format='%(asctime)s %(levelname)s %(message)s') n_row, n_col = 2, 3 n_components = n_row * n_col image_shape = (64, 64) rng = RandomState(0) # ############################################################################# # Load faces data dataset = fetch_olivetti_faces(shuffle=True, random_state=rng) faces = dataset.data n_samples, n_features = faces.shape # global centering faces_centered = faces - faces.mean(axis=0) # local centering faces_centered -= faces_centered.mean(axis=1).reshape(n_samples, -1) print("Dataset consists of %d faces" % n_samples) def plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray): plt.figure(figsize=(2. * n_col, 2.26 * n_row)) plt.suptitle(title, size=16) for i, comp in enumerate(images): plt.subplot(n_row, n_col, i + 1) vmax = max(comp.max(), -comp.min()) plt.imshow(comp.reshape(image_shape), cmap=cmap, interpolation='nearest', vmin=-vmax, vmax=vmax) plt.xticks(()) plt.yticks(()) plt.subplots_adjust(0.01, 0.05, 0.99, 0.93, 0.04, 0.) # ############################################################################# # PCA estimators = [ ('Eigenfaces - PCA using randomized SVD', decomposition.PCA(n_components=n_components, svd_solver='randomized', whiten=True), True) ] # ############################################################################# # Plot a sample of the input data plot_gallery("First centered Olivetti faces", faces_centered[:n_components]) # ############################################################################# # Do the estimation and plot it for name, estimator, center in estimators: print("Extracting the top %d %s..." % (n_components, name)) t0 = time() data = faces if center: data = faces_centered estimator.fit(data) train_time = (time() - t0) print("done in %0.3fs" % train_time) # print the components_ components_ = estimator.components_ print('components_:', components_.shape, 'n**n', components_) if hasattr(estimator, 'cluster_centers_'): components_ = estimator.cluster_centers_ else: components_ = estimator.components_ # Plot an image representing the pixelwise variance provided by the # estimator e.g its noise_variance_ attribute. The Eigenfaces estimator, # via the PCA decomposition, also provides a scalar noise_variance_ # (the mean of pixelwise variance) that cannot be displayed as an image # so we skip it. if (hasattr(estimator, 'noise_variance_') and estimator.noise_variance_.ndim > 0): # Skip the Eigenfaces case plot_gallery("Pixelwise variance", estimator.noise_variance_.reshape(1, -1), n_col=1, n_row=1) plot_gallery('%s - Train time %.1fs' % (name, train_time), components_[:n_components]) plt.show()
原始人臉影像中的前6張（普通人臉）
PCA從原始影像中提取出的6個獨立通用特徵（大眾臉特徵）
主成分矩陣
從運行結果中可以看到，PCA並沒有任何關於人臉的先驗資訊，但是PCA分解後的主成分矩陣，還是基本呈現了一個人臉的輪廓，也就是說，PCA降維後得到的主成分矩陣，代表了原始人臉影像數據集中的通用特徵。
主成分矩陣的每一行的維數都和原始影像的維數相同，區別在於不同像素點的權重分配不同，這點和NMF是一樣的。
筆者提醒：
嚴格來說，PCA是一種降維演算法，並不是一種獨立特徵提取演算法，但是在滿足特定前提假設條件下，PCA可以作為一種獨立特徵提取演算法，
- 數據集中不同特徵之間是彼此線性無關，或線性相關度不大的
- 數據集中包含的噪音數據不能太大，否則PCA很可能將噪音數據當做特徵提取出來
Relevant Link:
https://www.cnblogs.com/LittleHann/p/6558575.html#_label3_1_4_1 https://scikit-learn.org/stable/auto_examples/decomposition/plot_faces_decomposition.html#sphx-glr-auto-examples-decomposition-plot-faces-decomposition-py
4. SVD與NMF之間的異同分析
前面我們已經各自討論了SVD和NMF演算法的原理和應用，從表面形式上看，它們二者似乎都是在做矩陣分解，這章我們來梳理總結一下它們二者各自的異同。
0x1：特徵矩陣可解釋性方面的差異
一般把SVD和NMF分解得到的U矩陣中的列稱為特徵向量，從形式上看，二者的特徵向量有明顯的區別：
- NMF的特徵向量由於有非負的特點，因此不同向量之間的內積比大於零，不可能完全正交，說明NMF分解的特徵向量存在一定資訊冗餘。同時非負性也帶來了很好的物理解釋性，比如
- 影像領域，一個非負特徵向量可以解釋為一副特徵圖，每個元素代表一個像素
- 文本領域，一個非負特徵向量可以解釋為一個”主題“，每個元素代表某個單詞在主題中的重要程度
- SVD分解的特徵向量彼此正交，但失去了非負的特點，可解釋性變差
0x2：特徵矩陣提供對原數據的近似逼近能力的差異
- SVD分解的特徵向量在不同的維度都能夠對樣本向量做最好的近似，使近似後的結果與樣本的誤差在歐氏距離定義下達到全局最小
- NMF在分解過程中基數目需要預定確定，每次分解只是在一個確定維度下對樣本的局部近似，其與樣本的近似程度只能在歐式距離上達到局部最小，並且受初始值影響明顯
0x3：樣本特徵提取效率方面的差異
- 在體現樣本特徵方面，SVD所抽取的ui對應最大奇異值σi，則樣本在ui上的投影能量綜合為Ei = (XTui)T(XTui) = σi2。很明顯，對應奇異值越大的向量具有更大的樣本投影能量，不同的特徵ui在體現樣本特徵重要性方面有逐次之分
- NMF所抽取的特徵向量彼此間重要程度差別不大，無主次之分
總體來說，在Forebenius範數意義下，SVD分解具有全局最優意義上的處理能力，但數據的可解釋性不強；而NMF需要迭代計算，運算速度較慢，但分解的可解釋性強，帶有輔助約束的演算法可使分解的矩陣有更大的稀疏性。另外，NMF分解針對不同類型的數據進行處理性能差異明顯。
Relevant Link:
http://xueshu.baidu.com/usercenter/paper/show?paperid=2e99407850fe26f810cd1ddfa556caff&site=xueshu_se
5. 獨立成分分析（Independent component analysis，ICA）
0x1：從雞尾酒會問題（cocktail party problem）說起
雞尾酒問題是一個經典的盲源訊號分離（blind signal separation）問題。
假設在一個開party的房間里有兩個人同時說話，房間里兩個不同的位置上各放一個麥克風，各自記錄了一段聲音（時間訊號），x1(t)和x2(t)，這兩段聲音訊號都包含了兩個人說的話的混合，我們怎樣得到他們分別說了什麼話呢？
記源訊號是s1(t)和s2(t)，這個問題可以用以下等式描述：
如果s1，s2是統計獨立的隨機訊號，那麼就可以用x1，x2去還原。ICA正是利用統計上的獨立性從混合訊號中分解出源訊號的。
在討論混合矩陣、權重特徵這些具體概念，我們先將aij形象地理解是說話者跟麥克風的距離，如下圖，
現在假設聲源訊號分別為”sinusoidal signal“和”square signal“，將其混合為一個混合訊號，並加入一定的噪音干擾，來看一下PCA和ICA訊號分解的結果，
import numpy as np import matplotlib.pyplot as plt from scipy import signal from sklearn.decomposition import FastICA, PCA # ############################################################################# # Generate sample data np.random.seed(0) n_samples = 2000 time = np.linspace(0, 8, n_samples) s1 = np.sin(2 * time) # Signal 1 : sinusoidal signal s2 = np.sign(np.sin(3 * time)) # Signal 2 : square signal S = np.c_[s1, s2] S += 0.2 * np.random.normal(size=S.shape) # Add noise S /= S.std(axis=0) # Standardize data # Mix data A = np.array([[1, 1], [0.5, 2]]) # Mixing matrix X = np.dot(S, A.T) # Generate observations # print X print "X: ", X # Compute ICA ica = FastICA(n_components=2) S_ = ica.fit_transform(X) # Reconstruct signals A_ = ica.mixing_ # Get estimated mixing matrix # We can `prove` that the ICA model applies by reverting the unmixing. assert np.allclose(X, np.dot(S_, A_.T) + ica.mean_) # For comparison, compute PCA pca = PCA(n_components=2) H = pca.fit_transform(X) # Reconstruct signals based on orthogonal components # ############################################################################# # Plot results plt.figure() models = [X, S, S_, H] names = ['Observations (mixed signal)', 'True Sources', 'ICA recovered signals', 'PCA recovered signals'] colors = ['red', 'steelblue', 'orange'] for ii, (model, name) in enumerate(zip(models, names), 1): plt.subplot(4, 1, ii) plt.title(name) for sig, color in zip(model.T, colors): plt.plot(sig, color=color) plt.subplots_adjust(0.09, 0.04, 0.94, 0.94, 0.26, 0.46) plt.show()
可以看到，ICA和PCA都成功地對原始混合訊號進行了訊號分解，得到了原始的信源分訊號。
0x2：ICA演算法形式化定義
類似於EM演算法的思想，我們用隱變數模型（latent variables model）來描述ICA：
- x：由n個觀測訊號x1，x2，…，xn組成，x是可觀測的訊號序列
- s：由n個獨立訊號成分s1，s2，…，sn組成，s是不可觀測的隱變數序列
- A：n*n的混合矩陣，用於對隱變數序列進行混合疊加，即
ICA模型成立的基本假設如下：
- si之間是統計獨立的，它們是隱含變數，不能被直接觀測到，A和s都是未知的
- si服從非高斯分布
- 混合矩陣A可逆
ICA演算法的核心目標就是，基於觀測得到的x序列，估計A矩陣和s隱序列。
0x3：ICA演算法的基本前提
ICA可分解的前提是，隱變數序列si彼此獨立，而高斯的不相關和獨立是等價的，所以我們這裡討論高斯不相關性。
假設混合矩陣是方陣，如果s1和s2都是服從高斯分布，那麼x1，x2也是高斯的，且不相關，方差為1 。聯合概率密度如下：
概率分布圖如下，
從圖中可以看出，分布是對稱的，這樣就失去了方向資訊，也就是A的列向量的資訊。A也就預測不出來了。所謂的矩陣分解，本質就是在提取矩陣中的特徵向量的方向資訊。
可以證明，一個高斯分布經過正交變換以後仍然是那個高斯分布，而且變數之間仍是獨立的，所以如果變數都是高斯分布的，我們只能知道A是一個正交變換，不能完全確定A。
所以，對數據集進行ICA分解的大前提是，數據集獨立訊號之間需要彼此獨立。
0x4：ICA估計的理論推導
我們假設每個si有概率密度ps，那麼給定時刻原訊號（未知）的聯合分布就是，
這個公式代表一個前提：每個發聲源都是各自獨立的。有了ps(si），我們就可以求得觀測序列的聯合概率分布p(x)，
左邊是每個取樣訊號x（n維向量）的概率，右邊是每個原訊號概率的乘積|W|倍。
現在，我們必須對si假設一個概率分布函數F(x)，才能繼續求解W和s（EM演算法中的E步驟）。但是F(x)不能是隨機任意形式的，要滿足兩個性質：
- 單調遞增和在[0,1]區間內
- 不能是高斯分布的概率密度，ICA可分解的大前提就是，si序列彼此之間互相獨立
研究表明，sigmoid函數很合適，它的定義域為負無窮到正無窮，值域為0到1，同時還是一個對稱函數（即期望/均值為0）。
其導數為：
現在
我們就求W了（EM演算法中的M步驟），在給定訓練樣本{x(i)(x(i)1,x(i)2,x(i)3,…,,x(i)n);i = 1,…,m}後，對W的對數似然估計如下：
最後求導之後的結果公式如下，
其中α是梯度上升速率，當迭代求出W後，便可通過：s(i)=Wx(i)求出原始訊號，原始訊號的列向量即為獨立分解的分訊號。
0x5：演算法運行流程
- 數據前置預處理
- 中心化：將數據零均值化，使數據處於坐標中心
- 白化/球化：白化是一種重要的預處理過程，其目的是為了降低輸入數據的冗餘性，使得經過白化處理的輸入數據（i）特徵之間相關性較低 （ii）所有特徵具有相同的方差
- Expectation（E步驟）
- 假定si的概率密度ps，給定某時刻原訊號的聯合概率分布
- 求解p(x)
- 假定si的累計分布函數形式，例如sigmoid
- maximization（M步驟）
- 基於對si的函數形式假設，對W進行最大似然估計
- 似然函數對W求偏導，得到梯度
0x6：FastICA/固定點（Fixed-Point）演算法 – ICA的衍生優化演算法
我們知道，ICA演算法的核心目的是從觀測序列中得到一組彼此相對獨立的si序列。因此度量si序列提取效果的好壞是衡量其獨立性。
根據中心極限定理，若一隨機變數X由許多相互獨立的隨機變數Si（i=1,2,…,N）之和混合組成，只要Si具有有限的均值和方差，則不論其為何種分布，隨機變數X較Si更接近高斯分布。因此，在分離過程匯總，可通過對分離結果的非高斯性度量來表示分離結果間的互相獨立性，當非高斯性度量達到最大時，則表明已完成對各個獨立分量的分離。
度量費高斯性的方法有如下幾種：
- 峭度
- 似然最大
- 負熵最大
Relevant Link:
https://blog.csdn.net/u014485485/article/details/78452820 https://danieljyc.github.io/2014/06/13/%E6%9C%BA%E5%99%A8%E5%AD%A6%E4%B9%A015-3--%E7%8B%AC%E7%AB%8B%E6%88%90%E5%88%86%E5%88%86%E6%9E%90ica%EF%BC%88independent-component-analysis%EF%BC%89/ https://www.dataivy.cn/blog/%E7%8B%AC%E7%AB%8B%E6%88%90%E5%88%86%E5%88%86%E6%9E%90independent-component-analysis_ica/ https://www.cnblogs.com/Determined22/p/6357291.html https://scikit-learn.org/stable/auto_examples/decomposition/plot_ica_blind_source_separation.html#sphx-glr-auto-examples-decomposition-plot-ica-blind-source-separation-py https://www.jianshu.com/p/de396e8cce15
6. 自動編碼器（Auto-Encoder）
一言以蔽之，自動編碼器是一種無監督的神經網路模型，它可以用於像
- 學習輸入數據隱含特徵（編碼(coding)）
- 獨立特徵提取：自動編碼器學習到的新特徵可以送入有監督學習模型中
- 特徵降維：類似主成分分析PCA
- 生成與訓練樣本不同的新數據：這樣自動編碼器（變分自動編碼器，Variational Autoencoders）即為生成式模型
- 數據降噪：通過保留核心有效特徵，將原始數據中的噪點去除
之類的目的，同時AE用學習到的新特徵可以重構出原始輸入數據，稱之為解碼(decoding)。
一個最基本的AE架構如下圖，包括編碼和解碼兩個過程，
自動編碼器的編碼與解碼
自動編碼器是將輸入 進行編碼，得到新的特徵
，並且希望原始的輸入
能夠從新的特徵
重構出來。編碼過程如下：
可以看到，和神經網路結構一樣，其編碼過程本質上線性組合之後加上非線性的激活函數。
同時，利用新的特徵 ，可以對輸入
重構，即解碼過程：
我們希望重構出的 和儘可能一致，可以採用最小化負對數似然的損失函數來訓練這個模型：
一般情況下，我們會對自動編碼器加上一些限制，常用的是使 ，這稱為綁定權重（tied weights），本文所有的自動編碼器都加上這個限制。有時候，我們還會給自動編碼器加上更多的約束條件，去噪自動編碼器以及稀疏自動編碼器就屬於這種情況，因為大部分時候單純地重構原始輸入並沒有什麼意義，我們希望自動編碼器在近似重構原始輸入的情況下能夠捕捉到原始輸入更有價值的資訊，在這方面，embedding詞向量嵌入有類似的思想。
這章我們介紹幾種不同類型的Auto-Encoder，從不同角度理解AE的思想內涵。
0x1：全連接神經網路自編碼器
# -*- coding: utf-8 -*- from keras.layers import Input, Dense from keras.models import Model from keras.datasets import mnist import numpy as np import matplotlib.pyplot as plt ######################## 先建立一個全連接的編碼器和解碼器 ######################## # this is the size of our encoded representations encoding_dim = 32 # 32 floats -> compression of factor 24.5, assuming the input is 784 floats # this is our input placeholder input_img = Input(shape=(784,)) # "encoded" is the encoded representation of the input encoded = Dense(encoding_dim, activation='relu')(input_img) # "decoded" is the lossy reconstruction of the input decoded = Dense(784, activation='sigmoid')(encoded) # this model maps an input to its reconstruction autoencoder = Model(input=input_img, output=decoded) ######################## 可以單獨的使用編碼器和解碼器 # this model maps an input to its encoded representation encoder = Model(input=input_img, output=encoded) # create a placeholder for an encoded (32-dimensional) input encoded_input = Input(shape=(encoding_dim,)) # retrieve the last layer of the autoencoder model decoder_layer = autoencoder.layers[-1] # create the decoder model decoder = Model(input=encoded_input, output=decoder_layer(encoded_input)) ######################## 訓練自編碼器，來重構MNIST中的數字，這裡使用逐像素的交叉熵作為損失函數，優化器為adam autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy') ######################## 準備MNIST數據，將其歸一化和向量化，然後訓練 (x_train, _), (x_test, _) = mnist.load_data() x_train = x_train.astype('float32') / 255. x_test = x_test.astype('float32') / 255. x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:]))) x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:]))) print x_train.shape print x_test.shape autoencoder.fit(x_train, x_train, nb_epoch=50, batch_size=256, shuffle=True, validation_data=(x_test, x_test)) ######################## 可視化一下重構出來的輸出 # encode and decode some digits # note that we take them from the *test* set encoded_imgs = encoder.predict(x_test) decoded_imgs = decoder.predict(encoded_imgs) n = 10 # how many digits we will display plt.figure(figsize=(20, 4)) for i in range(1, n): # display original ax = plt.subplot(2, n, i) plt.imshow(x_test[i].reshape(28, 28)) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) # display reconstruction ax = plt.subplot(2, n, i + n) plt.imshow(decoded_imgs[i].reshape(28, 28)) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) plt.show()
AE解碼後的結果如下，
可以看到，該全連接神經網路成功地從784維的影像矩陣中提取出了有效的32維低維特徵資訊，並且基於這32維特徵，幾乎無損地重建出了原始的輸入影像。從這個例子中，我們可以感受到一些深度學習可解釋性方面的內涵。
0x2：稀疏自編碼器：為碼字加上稀疏性約束
上一節中我們的隱層有32個神經元，這種情況下，一般而言自編碼器學到的是PCA的一個近似。但是如果我們對隱層單元施加稀疏性約束（正則化約束）的話，會得到更為緊湊的表達，只有一小部分神經元會被激活。在Keras中，我們可以通過添加一個activity_regularizer達到對某層激活值進行約束的目的：
# -*- coding: utf-8 -*- from keras.layers import Input, Dense from keras.models import Model from keras.datasets import mnist from keras import regularizers import numpy as np import matplotlib.pyplot as plt ######################## 先建立一個全連接的編碼器和解碼器 ######################## # this is the size of our encoded representations encoding_dim = 64 # 64 floats -> compression of factor 24.5, assuming the input is 784 floats input_img = Input(shape=(784,)) # add a Dense layer with a L1 activity regularizer encoded = Dense(encoding_dim, activation='relu', kernel_regularizer=regularizers.l1(10e-5))(input_img) decoded = Dense(784, activation='sigmoid')(encoded) autoencoder = Model(input=input_img, output=decoded) ######################## 可以單獨的使用編碼器和解碼器 # this model maps an input to its encoded representation encoder = Model(input=input_img, output=encoded) # create a placeholder for an encoded (32-dimensional) input encoded_input = Input(shape=(encoding_dim,)) # retrieve the last layer of the autoencoder model decoder_layer = autoencoder.layers[-1] # create the decoder model decoder = Model(input=encoded_input, output=decoder_layer(encoded_input)) ######################## 訓練自編碼器，來重構MNIST中的數字，這裡使用逐像素的交叉熵作為損失函數，優化器為adam autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy') ######################## 準備MNIST數據，將其歸一化和向量化，然後訓練 (x_train, _), (x_test, _) = mnist.load_data() x_train = x_train.astype('float32') / 255. x_test = x_test.astype('float32') / 255. x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:]))) x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:]))) print x_train.shape print x_test.shape autoencoder.fit(x_train, x_train, epochs=100, batch_size=256, shuffle=True, validation_data=(x_test, x_test)) ######################## 可視化一下重構出來的輸出 # encode and decode some digits # note that we take them from the *test* set encoded_imgs = encoder.predict(x_test) decoded_imgs = decoder.predict(encoded_imgs) n = 10 # how many digits we will display plt.figure(figsize=(20, 4)) for i in range(1, n): # display original ax = plt.subplot(2, n, i) plt.imshow(x_test[i].reshape(28, 28)) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) # display reconstruction ax = plt.subplot(2, n, i + n) plt.imshow(decoded_imgs[i].reshape(28, 28)) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) plt.show()
隱層設置為32維
隱層設置為64維
可以看到，隱層的維數越高，AE解碼還原後的影像辨識度就相對越高。但總的來說，增加了稀疏約束後，編碼出來的碼字（隱層特徵向量）都會更加稀疏。
在保持32維隱層的情況下，稀疏自編碼器的在10000個測試圖片上的碼字均值為3.33，而之前的為7.30。
0x3：深度自編碼器（DAE）：堆疊自動編碼器
我們可以將多個自編碼器疊起來，本質上是提高深度神經網路的複雜度，提升特徵提取效率，
# -*- coding: utf-8 -*- from keras.layers import Input, Dense from keras.models import Model from keras.datasets import mnist from keras import regularizers import numpy as np import matplotlib.pyplot as plt ######################## 先建立一個全連接的編碼器和解碼器 ######################## # this is the size of our encoded representations encoding_dim = 32 # 64 floats -> compression of factor 24.5, assuming the input is 784 floats input_img = Input(shape=(784,)) encoded = Dense(128, activation='relu')(input_img) encoded = Dense(64, activation='relu')(encoded) encoded = Dense(32, activation='relu')(encoded) decoded = Dense(64, activation='relu')(encoded) decoded = Dense(128, activation='relu')(decoded) decoded = Dense(784, activation='sigmoid')(decoded) autoencoder = Model(input=input_img, output=decoded) autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy') ######################## 訓練自編碼器，來重構MNIST中的數字，這裡使用逐像素的交叉熵作為損失函數，優化器為adam autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy') ######################## 準備MNIST數據，將其歸一化和向量化，然後訓練 (x_train, _), (x_test, _) = mnist.load_data() x_train = x_train.astype('float32') / 255. x_test = x_test.astype('float32') / 255. x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:]))) x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:]))) print x_train.shape print x_test.shape autoencoder.fit(x_train, x_train, epochs=100, batch_size=128, shuffle=True, validation_data=(x_test, x_test)) ######################## 可視化一下重構出來的輸出 # encode and decode some digits # note that we take them from the *test* set decoded_imgs = autoencoder.predict(x_test) n = 10 # how many digits we will display plt.figure(figsize=(20, 4)) for i in range(1, n): # display original ax = plt.subplot(2, n, i) plt.imshow(x_test[i].reshape(28, 28)) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) # display reconstruction ax = plt.subplot(2, n, i + n) plt.imshow(decoded_imgs[i].reshape(28, 28)) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) plt.show()
效果比AE要好一些，當然，訓練時間也增加了。
0x4：卷積自編碼器：用卷積層搭建自編碼器
當輸入是影像時，使用卷積神經網路基本是一個較好的策略，在工程項目中，用於處理影像的自動編碼器幾乎都是卷積自動編碼器。
- 卷積自編碼器的編碼器部分由卷積層和MaxPooling層構成
- MaxPooling負責空域下取樣
- 卷積層負責提取細節特徵資訊
- 而解碼器由卷積層和上取樣層構成
- 卷積層負責資訊傳遞
- 上取樣層負責將卷積資訊反卷回原始訊號
# -*- coding: utf-8 -*- from keras.layers import Input, Dense from keras.models import Model from keras.datasets import mnist import numpy as np from keras.layers import Input, Dense, Convolution2D, MaxPooling2D, UpSampling2D from keras.models import Model import matplotlib.pyplot as plt ######################## 先建立一個編碼器和解碼器 ######################## input_img = Input(shape=(1, 28, 28)) x = Convolution2D(16, 3, 3, activation='relu', border_mode='same')(input_img) x = MaxPooling2D((2, 2), border_mode='same')(x) x = Convolution2D(8, 3, 3, activation='relu', border_mode='same')(x) x = MaxPooling2D((2, 2), border_mode='same')(x) x = Convolution2D(8, 3, 3, activation='relu', border_mode='same')(x) encoded = MaxPooling2D((2, 2), border_mode='same')(x) # at this point the representation is (8, 4, 4) i.e. 128-dimensional x = Convolution2D(8, 3, 3, activation='relu', border_mode='same')(encoded) x = UpSampling2D((2, 2))(x) x = Convolution2D(8, 3, 3, activation='relu', border_mode='same')(x) x = UpSampling2D((2, 2))(x) x = Convolution2D(16, 3, 3, activation='relu')(x) x = UpSampling2D((2, 2))(x) decoded = Convolution2D(1, 3, 3, activation='sigmoid', border_mode='same')(x) autoencoder = Model(input_img, decoded) autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy') ######################## 訓練自編碼器，來重構MNIST中的數字，這裡使用逐像素的交叉熵作為損失函數，優化器為adam autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy') ######################## 準備MNIST數據，將其歸一化和向量化，然後訓練 (x_train, _), (x_test, _) = mnist.load_data() x_train = x_train.astype('float32') / 255. x_test = x_test.astype('float32') / 255. x_train = np.reshape(x_train, (len(x_train), 1, 28, 28)) x_test = np.reshape(x_test, (len(x_test), 1, 28, 28)) print x_train.shape print x_test.shape #print "x_test: ", x_test autoencoder.fit(x_train, x_train, nb_epoch=100, batch_size=128, shuffle=True, validation_data=(x_test, x_test) ) ######################## 可視化一下重構出來的輸出 # encode and decode some digits # note that we take them from the *test* set decoded_imgs = autoencoder.predict(x_test) n = 10 # how many digits we will display plt.figure(figsize=(20, 4)) for i in range(1, n): # display original ax = plt.subplot(2, n, i) plt.imshow(x_test[i].reshape(28, 28)) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) # display reconstruction ax = plt.subplot(2, n, i + n) plt.imshow(decoded_imgs[i].reshape(28, 28)) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) plt.show()
0x5：AE隱層的特徵向量可視化
這個小節，不介紹新的AE，我們來看看中間的碼字長什麼樣，我們將中間碼字reshape成4*32，
# -*- coding: utf-8 -*- from keras.layers import Input, Dense from keras.models import Model from keras.datasets import mnist from keras import regularizers import numpy as np import matplotlib.pyplot as plt ######################## 先建立一個編碼器和解碼器 ######################## # this is the size of our encoded representations encoding_dim = 128 # 64 floats -> compression of factor 24.5, assuming the input is 784 floats input_img = Input(shape=(784,)) # add a Dense layer with a L1 activity regularizer encoded = Dense(encoding_dim, activation='relu', kernel_regularizer=regularizers.l1(10e-5))(input_img) decoded = Dense(784, activation='sigmoid')(encoded) autoencoder = Model(input=input_img, output=decoded) ######################## 可以單獨的使用編碼器和解碼器 # this model maps an input to its encoded representation encoder = Model(input=input_img, output=encoded) # create a placeholder for an encoded (32-dimensional) input encoded_input = Input(shape=(encoding_dim,)) # retrieve the last layer of the autoencoder model decoder_layer = autoencoder.layers[-1] # create the decoder model decoder = Model(input=encoded_input, output=decoder_layer(encoded_input)) ######################## 訓練自編碼器，來重構MNIST中的數字，這裡使用逐像素的交叉熵作為損失函數，優化器為adam autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy') ######################## 準備MNIST數據，將其歸一化和向量化，然後訓練 (x_train, _), (x_test, _) = mnist.load_data() x_train = x_train.astype('float32') / 255. x_test = x_test.astype('float32') / 255. x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:]))) x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:]))) print x_train.shape print x_test.shape autoencoder.fit(x_train, x_train, epochs=100, batch_size=256, shuffle=True, validation_data=(x_test, x_test)) ######################## 可視化一下中間碼字 # encode and decode some digits # note that we take them from the *test* set encoded_imgs = encoder.predict(x_test) decoded_imgs = decoder.predict(encoded_imgs) n = 10 # how many digits we will display plt.figure(figsize=(20, 8)) for i in range(1,n): ax = plt.subplot(1, n, i) plt.imshow(encoded_imgs[i].reshape(4, 4 * 8).T) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) plt.show()
0x6：使用自動編碼器進行影像去噪
從神經網路結構的角度看，AE本質上是一個”matrix-to-matrix”的”結構的濾波器，我們在很多語言翻譯神經網路中也會看到類似結構。基於這種濾波器結構的AE網路，我們可以實現影像降噪方面的目的。
我們把訓練樣本用雜訊污染，然後使解碼器解碼出乾淨的照片，以獲得去噪自動編碼器。首先我們把原圖片加入高斯雜訊，然後把像素值clip到0~1。
可以看到，降噪效果很不錯，AE起到了濾波器的作用。
筆者思考：
在實際工程項目中，數據集中的噪點是一個很常見的問題，例如
- 影像中的白噪點
- 文本文檔中出現的一些不常見詞
- 由於取樣過程引入的系統性噪點數據
當然，演算法本身也有一些緩解手段來應對噪點問題（例如SVM中的soft-boundary），但是依舊無法完全規避數據噪點帶來的過擬合和欠擬合問題。一個可行的做法是，對原始數據先通過AE進行降噪處理後，再輸出機器學習模型（例如決策樹），這樣會得到更好的模型性能。
0x7：變分自編碼器（Variational autoencoder，VAE）：編碼數據集的概率分布
我們已經知道，AE是一個獨立特徵提取的神經網路，在前面的章節中，我們介紹了通過AE直接在隱層得到碼字特徵向量，這是一種一步到位的做法。這小節，我們對前面的AE進行一些小修改，使得編碼器學習到輸入數據的隱變數模型（概率分布模型）。隱變數模型是連接顯變數集和隱變數集的統計模型，隱變數模型的假設是顯變數是由隱變數的狀態控制的，各個顯變數之間條件獨立。
簡單來說，變分編碼器學習數據概率分布的一組參數（編碼階段）。並通過在這個概率分布中取樣，你可以生成新的輸入數據（解碼階段），即變分編碼器是一個生成模型。
變分編碼器的簡要工作過程如下：
- 首先，編碼器網路將輸入樣本x轉換為隱空間的兩個參數，記作z_mean和z_log_sigma，得到了這2個參數，就等於得到了隱藏的正態分布，z = z_mean + exp(z_log_sigma)*epsilon，epsilon是一個服從正態分布的張量。
- 然後，我們隨機從隱藏的正態分布中取樣得到數據點z。
- 最後，使用解碼器網路將隱空間映射到顯空間，即將z轉換回原來的輸入數據空間。
參數藉由兩個損失函數來訓練，
- 重構損失函數，該函數要求解碼出來的樣本與輸入的樣本相似（與之前的自編碼器相同）
- 學習到的隱分布與先驗分布的KL距離，作為一個正則化選項是可選的，它對學習符合要求的隱空間和防止過擬合有幫助
# -*- coding: utf-8 -*- from __future__ import absolute_import from __future__ import division from __future__ import print_function from keras.layers import Lambda, Input, Dense from keras.models import Model from keras.datasets import mnist from keras.losses import mse, binary_crossentropy from keras.utils import plot_model from keras import backend as K import numpy as np import matplotlib.pyplot as plt import argparse import os # reparameterization trick # instead of sampling from Q(z|X), sample epsilon = N(0,I) # z = z_mean + sqrt(var) * epsilon def sampling(args): """Reparameterization trick by sampling from an isotropic unit Gaussian. # Arguments args (tensor): mean and log of variance of Q(z|X) # Returns z (tensor): sampled latent vector """ z_mean, z_log_var = args batch = K.shape(z_mean)[0] dim = K.int_shape(z_mean)[1] # by default, random_normal has mean = 0 and std = 1.0 epsilon = K.random_normal(shape=(batch, dim)) return z_mean + K.exp(0.5 * z_log_var) * epsilon def plot_results(models, data, batch_size=128, model_name="vae_mnist"): """Plots labels and MNIST digits as a function of the 2D latent vector # Arguments models (tuple): encoder and decoder models data (tuple): test data and label batch_size (int): prediction batch size model_name (string): which model is using this function """ encoder, decoder = models x_test, y_test = data #os.makedirs(model_name, exist_ok=True) filename = os.path.join(model_name, "vae_mean.png") # display a 2D plot of the digit classes in the latent space z_mean, _, _ = encoder.predict(x_test, batch_size=batch_size) plt.figure(figsize=(12, 10)) plt.scatter(z_mean[:, 0], z_mean[:, 1], c=y_test) plt.colorbar() plt.xlabel("z[0]") plt.ylabel("z[1]") plt.savefig(filename) plt.show() filename = os.path.join(model_name, "digits_over_latent.png") # display a 30x30 2D manifold of digits n = 30 digit_size = 28 figure = np.zeros((digit_size * n, digit_size * n)) # linearly spaced coordinates corresponding to the 2D plot # of digit classes in the latent space grid_x = np.linspace(-4, 4, n) grid_y = np.linspace(-4, 4, n)[::-1] for i, yi in enumerate(grid_y): for j, xi in enumerate(grid_x): z_sample = np.array([[xi, yi]]) x_decoded = decoder.predict(z_sample) digit = x_decoded[0].reshape(digit_size, digit_size) figure[i * digit_size: (i + 1) * digit_size, j * digit_size: (j + 1) * digit_size] = digit plt.figure(figsize=(10, 10)) start_range = digit_size // 2 end_range = (n - 1) * digit_size + start_range + 1 pixel_range = np.arange(start_range, end_range, digit_size) sample_range_x = np.round(grid_x, 1) sample_range_y = np.round(grid_y, 1) plt.xticks(pixel_range, sample_range_x) plt.yticks(pixel_range, sample_range_y) plt.xlabel("z[0]") plt.ylabel("z[1]") plt.imshow(figure, cmap='Greys_r') plt.savefig(filename) plt.show() # MNIST dataset # 使用MNIST庫來訓練變分編碼器 (x_train, y_train), (x_test, y_test) = mnist.load_data() image_size = x_train.shape[1] original_dim = image_size * image_size x_train = np.reshape(x_train, [-1, original_dim]) x_test = np.reshape(x_test, [-1, original_dim]) x_train = x_train.astype('float32') / 255 x_test = x_test.astype('float32') / 255 # network parameters input_shape = (original_dim, ) intermediate_dim = 512 batch_size = 128 latent_dim = 2 epochs = 50 # VAE model = encoder + decoder # build encoder model # 建立編碼網路，將輸入影射為隱分布的參數： inputs = Input(shape=input_shape, name='encoder_input') x = Dense(intermediate_dim, activation='relu')(inputs) z_mean = Dense(latent_dim, name='z_mean')(x) z_log_var = Dense(latent_dim, name='z_log_var')(x) # use reparameterization trick to push the sampling out as input # note that "output_shape" isn't necessary with the TensorFlow backend # 從這些參數確定的分布中取樣，這個樣本相當於之前的隱層值 z = Lambda(sampling, output_shape=(latent_dim,), name='z')([z_mean, z_log_var]) # instantiate encoder model encoder = Model(inputs, [z_mean, z_log_var, z], name='encoder') encoder.summary() plot_model(encoder, to_file='vae_mlp_encoder.png', show_shapes=True) # build decoder model # 將取樣得到的點映射回去重構原輸入 latent_inputs = Input(shape=(latent_dim,), name='z_sampling') x = Dense(intermediate_dim, activation='relu')(latent_inputs) outputs = Dense(original_dim, activation='sigmoid')(x) # instantiate decoder model decoder = Model(latent_inputs, outputs, name='decoder') decoder.summary() plot_model(decoder, to_file='vae_mlp_decoder.png', show_shapes=True) # instantiate VAE model # 到目前為止我們做的工作需要實例化三個模型： # 1. 一個端到端的自動編碼器，用於完成輸入訊號的重構 # 2. 一個用於將輸入空間映射為隱空間的編碼器 # 3. 一個利用隱空間的分布產生的樣本點生成對應的重構樣本的生成器 outputs = decoder(encoder(inputs)[2]) vae = Model(inputs, outputs, name='vae_mlp') if __name__ == '__main__': parser = argparse.ArgumentParser() help_ = "Load h5 model trained weights" parser.add_argument("-w", "--weights", help=help_) help_ = "Use mse loss instead of binary cross entropy (default)" parser.add_argument("-m", "--mse", help=help_, action='store_true') args = parser.parse_args() models = (encoder, decoder) data = (x_test, y_test) # VAE loss = mse_loss or xent_loss + kl_loss if args.mse: reconstruction_loss = mse(inputs, outputs) else: reconstruction_loss = binary_crossentropy(inputs, outputs) reconstruction_loss *= original_dim kl_loss = 1 + z_log_var - K.square(z_mean) - K.exp(z_log_var) kl_loss = K.sum(kl_loss, axis=-1) kl_loss *= -0.5 vae_loss = K.mean(reconstruction_loss + kl_loss) vae.add_loss(vae_loss) vae.compile(optimizer='adam') vae.summary() plot_model(vae, to_file='vae_mlp.png', show_shapes=True) if args.weights: vae.load_weights(args.weights) else: # train the autoencoder vae.fit(x_train, epochs=epochs, batch_size=batch_size, validation_data=(x_test, None)) vae.save_weights('vae_mlp_mnist.h5') plot_results(models, data, batch_size=batch_size, model_name="vae_mlp") x_test_encoded = encoder.predict(x_test, batch_size=batch_size) plt.figure(figsize=(6, 6)) plt.scatter(x_test_encoded[:, 0], x_test_encoded[:, 1], c=y_test) plt.colorbar() plt.show() # display a 2D manifold of the digits n = 15 # figure with 15x15 digits digit_size = 28 figure = np.zeros((digit_size * n, digit_size * n)) # we will sample n points within [-15, 15] standard deviations grid_x = np.linspace(-15, 15, n) grid_y = np.linspace(-15, 15, n) for i, yi in enumerate(grid_x): for j, xi in enumerate(grid_y): z_sample = np.array([[xi, yi]]) * [1.0, 1.0] x_decoded = decoder.predict(z_sample) digit = x_decoded[0].reshape(digit_size, digit_size) figure[i * digit_size: (i + 1) * digit_size, j * digit_size: (j + 1) * digit_size] = digit plt.figure(figsize=(10, 10)) plt.imshow(figure) plt.show()
編碼器網路結構
解碼器網路結構
總體VAE網路結構
因為我們的隱空間只有兩維，所以我們可以可視化一下。我們來看看2D平面中不同類的近鄰分布：
上圖每種顏色代表一個數字，相近聚類的數字代表他們在結構上相似。
因為變分編碼器是一個生成模型，我們可以用它來生成新數字。我們可以從隱平面上取樣一些點，然後生成對應的顯變數，即MNIST的數字：
可以看到，和上面的2D概率是對應的，不同的數字佔據了不同的概率分布區域。
筆者思考：
從這個實驗中，我們可以更加深刻的理解深度神經網路的可解釋性的內涵，基本上來說，深度神經網路強大的擬合能力來源於它的高維非線性，同時因為正則化約束的存在，隱層中間態碼字的權重分配會傾向於集中。這樣，神經網路就可以同時具備對複雜問題和簡單問題的擬合能力，總體上上說，深度神經網路的策略是：”盡量精確匹配，但是夠用就行“。所以在MNIST這個問題上，AE只需要很少的隱層碼字就可以提取出MNIST手寫數字中的主要特徵，但是如果面對的是一個更複雜的問題，AE就相對地需要更多的碼字來進行有效編碼。
相關的思想，田淵棟博士在它的論文《Luck Matters: Understanding Training Dynamics of Deep ReLU Networks》中進行了詳細討論，論文非常精彩，推薦讀者朋友們。
Relevant Link:
https://blog.csdn.net/a819825294/article/details/53516980 https://zhuanlan.zhihu.com/p/31742653 https://zhuanlan.zhihu.com/p/67782029 https://github.com/keras-team/keras/blob/master/examples/variational_autoencoder.py https://keras-cn.readthedocs.io/en/latest/legacy/blog/autoencoder/
7. AE（Auto-Encoder）和PCA（principle component analysis）在降維效果方面的對比
這個章節，我們來做一個有趣的實驗，通過影像處理這個具體問題，來探究AE和PCA在降維和獨立特徵提取方面性能的優劣。
0x1：AE和PCA各自的理論區別
- AE的理論基礎是低維複合非線性函數對高維數據流形的近似逼近
- PCA的理論基礎是最大熵定理，通過尋找熵值最大的top K個特徵向量基，從而得到一個新的K維向量空間，將原始數據集投影到該新向量空間中，從而實現降維
從線性代數角度來看，這兩種技術似乎有些異曲同工之妙，我們來通過實驗驗證一下這個猜想。
0x2：PCA和AE對MNIST手寫數字的降維的效果對比
我們分別調用PCA和AE的API，將MNIST手寫數字數據測試集降維到128維，並觀察前10個降維後的影像，
# -*- coding: utf-8 -*- import numpy as np import pandas as pd from pandas import Series,DataFrame import matplotlib.pyplot as plt from sklearn.decomposition import PCA from keras.datasets import mnist from keras.layers import Input, Dense from keras.models import Model from keras import regularizers if __name__ == '__main__': # load mnist data (x_train, _), (x_test, _) = mnist.load_data() x_train = x_train.astype('float32') / 255. x_test = x_test.astype('float32') / 255. x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:]))) x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:]))) print x_train.shape print x_test.shape # 1. PCA: dimensionality reduction into 128 dim pca = PCA(n_components=128, whiten=True) pca.fit(x_train, x_train) X_train_pca = pca.transform(x_train) x_test_pca = pca.transform(x_test) print "x_test_pca: ", x_test_pca print "np.shape(x_test_pca): ", np.shape(x_test_pca) # show top100 sample which is been pca reduction samples = x_test_pca[:8] plt.figure(figsize=(12, 18)) for i in range(len(samples)): plt.subplot(10, 10, i + 1) plt.imshow(samples[i].reshape(32, 4), cmap='gray') plt.axis('off') plt.show() # 2. AE: dimensionality reduction into 128 dim encoding_dim = 128 input_img = Input(shape=(784,)) encoded = Dense(encoding_dim, activation='relu')(input_img) decoded = Dense(784, activation='sigmoid')(encoded) autoencoder = Model(input=input_img, output=decoded) encoder = Model(input=input_img, output=encoded) encoded_input = Input(shape=(encoding_dim,)) autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy') autoencoder.fit(x_train, x_train, epochs=100, batch_size=256, shuffle=True, validation_data=(x_test, x_test)) encoded_imgs = encoder.predict(x_test) n = 10 # how many digits we will display plt.figure(figsize=(20, 8)) for i in range(1, n): ax = plt.subplot(1, n, i) plt.imshow(encoded_imgs[i].reshape(4, 32).T) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) plt.show()
PCA降維結果
AE降維結果
從肉眼上看，兩種技術的降維效果差不多。那怎麼定量評估這兩種方案的實際效果（保真度和可還原度）呢？我們可以將它們降維後的數據集輸入有監督學習演算法（例如SVC），並評估有監督模型的預測性能，以此來評價前置降維演算法的性能優劣。
# -*- coding: utf-8 -*- import numpy as np import pandas as pd from pandas import Series,DataFrame import matplotlib.pyplot as plt from sklearn.decomposition import PCA from keras.datasets import mnist from keras.layers import Input, Dense from keras.models import Model from keras import regularizers from sklearn.svm import SVC from sklearn.metrics import accuracy_score if __name__ == '__main__': # load mnist data (x_train, y_train), (x_test, y_test) = mnist.load_data() x_train = x_train.astype('float32') / 255. x_test = x_test.astype('float32') / 255. x_train = x_train.reshape((len(x_train), np.prod(x_train.shape[1:]))) x_test = x_test.reshape((len(x_test), np.prod(x_test.shape[1:]))) print x_train.shape print "y_train: ", y_train print x_test.shape print "y_test: ", y_test # 1. PCA: dimensionality reduction into 128 dim pca = PCA(n_components=128, whiten=True) pca.fit(x_train, y_train) X_train_pca = pca.transform(x_train) x_test_pca = pca.transform(x_test) print "x_test_pca: ", x_test_pca print "np.shape(x_test_pca): ", np.shape(x_test_pca) # show top100 sample which is been pca reduction samples = x_test_pca[:8] plt.figure(figsize=(12, 18)) for i in range(len(samples)): plt.subplot(10, 10, i + 1) plt.imshow(samples[i].reshape(32, 4), cmap='gray') plt.axis('off') plt.show() # crate supervised learning model svc = SVC(kernel='rbf') svc.fit(X_train_pca, y_train) print 'accuracy for PCA: {0}'.format(accuracy_score(y_test, svc.predict(x_test_pca))) # 2. AE: dimensionality reduction into 128 dim encoding_dim = 128 input_img = Input(shape=(784,)) encoded = Dense(encoding_dim, activation='relu')(input_img) decoded = Dense(784, activation='sigmoid')(encoded) autoencoder = Model(input=input_img, output=decoded) encoder = Model(input=input_img, output=encoded) encoded_input = Input(shape=(encoding_dim,)) autoencoder.compile(optimizer='adadelta', loss='binary_crossentropy') print "x_train: ", x_train autoencoder.fit(x_train, x_train, epochs=100, batch_size=256, shuffle=True, validation_data=(x_test, x_test)) X_train_ae = encoder.predict(x_train) x_test_ae = encoder.predict(x_test) n = 10 # how many digits we will display plt.figure(figsize=(20, 8)) for i in range(1, n): ax = plt.subplot(1, n, i) plt.imshow(X_train_ae[i].reshape(4, 32).T) plt.gray() ax.get_xaxis().set_visible(False) ax.get_yaxis().set_visible(False) plt.show() # crate supervised learning model svc = SVC(kernel='rbf') svc.fit(X_train_ae, y_train) print 'accuracy for AE: {0}'.format(accuracy_score(y_test, svc.predict(x_test_ae)))
分別用PCA和AE降維後的訓練集訓練得到的SVC有監督模型，對同樣降維後測試集進行預測（降維演算法這裡充當一個前置樣本預處理器），結果如下，
- PCA降維後的預測精確度：0.9782
- AE降維後的預測精確度：0.9779
從實驗結果來看，AE和PCA的降維效果不相上下，不存在誰比誰有質的飛躍。讀者朋友還可以用DAE代替PCA做同樣的實驗，結果也差不多，DAE的效果有輕微提升，但沒有質的飛躍。
從這個實驗中可以看出，我們不需要盲目崇拜深度學習，雖然現在深度學習是最熱門的一個前沿論文領域，但是有著堅實理論基礎的傳統機器學習演算法，依然有著強大的生命力以及不輸深度學習的效果。我們學習要注重學習底層原理，不要盲目追新。
8. 獨立特徵提取在細分領域的衍生演算法
本文介紹的幾種獨立特徵提取演算法基本上是和具體細分領域無關的通用演算法。其實在文本降維、文本語義理解、話題提取、影像處理、音頻訊號處理等領域還有很多衍生的獨立特徵提取演算法，例如：
- 影像處理領域
- 方向梯度直方圖（Histogram of Oriented Gradient, HOG）特徵：種在電腦視覺和影像處理中用來進行物體檢測的特徵描述子。它通過計算和統計影像局部區域的梯度方向直方圖來構成特徵
- SIFT（尺度不變特徵變換）：在不同的尺度空間上查找關鍵點(特徵點)，並計算出關鍵點的方向。SIFT所查找到的關鍵點是一些十分突出、不會因光照、仿射變換和噪音等因素而變化的點，如角點、邊緣點、暗區的亮點及亮區的暗點等
- SURF：SURF主要是把SIFT中的某些運算作了簡化。SURF把SIFT中的高斯二階微分的模板進行了簡化，使得卷積平滑操作僅需要轉換成加減運算，這樣使得SURF演算法的魯棒性好且時間複雜度低
- ORB：ORB特徵基於FAST角點的特徵點檢測與描述技術，具有尺度與旋轉不變性，同時對雜訊及透視仿射也具有不變性，良好的性能使得用ORB在進行特徵描述時的應用場景十分廣泛
- LBP(Local Binary Pattern)：局部二值模式是一種描述影像局部紋理的特徵運算元，具有旋轉不變性與灰度不變性等顯著優點。LBP特徵描述的是一種灰度範圍內的影像處理操作技術，針對的是輸入源為8位或16位的灰度影像
- 語義提取
- PLSA
- LDA
- HDP
- 向量空間模型（VSM）：通過為語義單元賦予不同的權重以反映它們在語義表達能力上的差異，將文本看作有一組正交詞條構成的矢量空間，將文本的語義單元看作是高維空間的維度
- wordembedding（詞向量）。VSM背後的數學基礎還是矩陣因式分解理論
- 概率模型：概率模型基於特徵的概率分布表示文本數據，同時也考慮特徵之間的概率依賴關係。例如二元獨立概率模型、二元一階相關概率模型、雙柏松分布概率模型以及概率網路資訊模型等
- 音頻訊號處理
- 過零率：過零率體現的是訊號過零點的次數，體現的是頻率特性。因為需要過零點，所以訊號處理之前需要中心化處理
- 短時平均幅度差：音頻具有周期特性，平穩雜訊情況下利用短時平均幅度差可以更好地觀察周期特性
9. 矩陣因式分解（獨立特徵提取）在安全領域的應用
本章內容由於涉及到一些工作上的內容和實驗結果，筆者這裡只放置一些簡要思路，詳細的實驗過程和實驗結果，有興趣的朋友可以給我發郵件或者微信，我們可以離線交流。
0x1：基於Webshell文檔獨立特徵提取的模型可解釋性研究
- 通過對詞向量矩陣進行因式分解，得到webshell黑樣本的整體核心特徵（特徵矩陣的列向量），以此獲得webshell檢測模型可解釋性方面的描述，即哪些詞（word）佔據了主要的權重分配
- 通過對詞向量矩陣進行因式分解，得到每篇文檔的特徵主題（權重矩陣），將每篇文檔（黑樣本/白樣本）都轉換為一個主題列表，用一個更高維抽象的視角來看webshell黑白樣本
程式碼里的blacksamples是筆者自己準備的一份webshell黑樣本，讀者朋友請自行準備。
# -*- coding: utf-8 -*- from time import time from sklearn.feature_extraction.text import TfidfVectorizer, CountVectorizer from sklearn.decomposition import NMF, LatentDirichletAllocation from sklearn.datasets import fetch_20newsgroups from sklearn.externals import joblib import os import pickle import pandas import numpy as np n_samples = 2000 n_features = 1000 n_components = 10 n_top_words = 20 def print_top_words(model, feature_names, n_top_words): for topic_idx, topic in enumerate(model.components_): message = "Topic #%d: " % topic_idx message += " ".join([feature_names[i] for i in topic.argsort()[:-n_top_words - 1:-1]]) print(message) print() def load_webshell_apidata(): op_dict = { 'apicalls': [], 'path': [] } rootDir = "./blacksamples" for lists in os.listdir(rootDir): if lists == '.DS_Store': continue pickle_path = os.path.join(rootDir, lists, "vector.pickle") if os.path.exists(pickle_path): op = pickle.load(open(pickle_path, 'rb')) op_dict['apicalls'].append(' '.join(op['apicalls'])) # add as a apicall string "api api" op_dict['path'].append(op['path']) df = pandas.DataFrame(data=op_dict) return df print("Loading dataset...") t0 = time() dataset_apiifo = load_webshell_apidata() dataset = dataset_apiifo['apicalls'].values.tolist() dataset_path = dataset_apiifo['path'].values.tolist() #print "dataset: ", dataset data_samples = dataset[:n_samples] print "data_samples[:10]: ", data_samples[:10] print("done in %0.3fs." % (time() - t0)) print(" ") # Use tf-idf features for NMF. print("Extracting tf-idf features for NMF...") tfidf_vectorizer = TfidfVectorizer(max_df=0.95, min_df=2, max_features=n_features, stop_words='english') t0 = time() tfidf = tfidf_vectorizer.fit_transform(data_samples) print("done in %0.3fs." % (time() - t0)) print(" ") # Use tf (raw term count) features for LDA. print("Extracting tf features for LDA...") tf_vectorizer = CountVectorizer(max_df=0.95, min_df=2, max_features=n_features, stop_words='english') t0 = time() tf = tf_vectorizer.fit_transform(data_samples) print("done in %0.3fs." % (time() - t0)) print(" ") # Fit the NMF model print("Fitting the NMF model (Frobenius norm) with tf-idf features, " "n_samples=%d and n_features=%d..." % (n_samples, n_features)) t0 = time() nmf = NMF(n_components=n_components, random_state=1, alpha=.1, l1_ratio=.5).fit(tfidf) print("done in %0.3fs." % (time() - t0)) print(" ") print("nTopics in NMF model (Frobenius norm):") tfidf_feature_names = tfidf_vectorizer.get_feature_names() print_top_words(nmf, tfidf_feature_names, n_top_words) # Fit the NMF model print("Fitting the NMF model (generalized Kullback-Leibler divergence) with " "tf-idf features, n_samples=%d and n_features=%d..." % (n_samples, n_features)) t0 = time() nmf = NMF(n_components=n_components, random_state=1, beta_loss='kullback-leibler', solver='mu', max_iter=1000, alpha=.1, l1_ratio=.5).fit(tfidf) print("done in %0.3fs." % (time() - t0)) print(" ") print("nTopics in NMF model (generalized Kullback-Leibler divergence):") tfidf_feature_names = tfidf_vectorizer.get_feature_names() print_top_words(nmf, tfidf_feature_names, n_top_words) print("Fitting LDA models with tf features, " "n_samples=%d and n_features=%d..." % (n_samples, n_features)) lda = LatentDirichletAllocation(n_components=n_components, max_iter=5, learning_method='online', learning_offset=50., random_state=0) t0 = time() lda.fit(tf) print("done in %0.3fs." % (time() - t0)) print(" ") print("nTopics in LDA model:") tf_feature_names = tf_vectorizer.get_feature_names() print_top_words(lda, tf_feature_names, n_top_words)
實驗結果如下：
Topics in NMF model (Frobenius norm): Topic #0: stristr eval preg_match headers_sent error_reporting base64_decode gzinflate mail isref explode include defined require trim getcwd getenv getdir getfakedownloadpath get_info getdoccomment Topic #1: stripos eval gzuncompress setcookie time base64_decode header require_once dirname defined gzinflate require mail explode trim str_ireplace file_exists ini_get rtrim ltrim Topic #2: function_exists eval create_function preg_match gzdecode ob_start header defined error_reporting base64_decode gzinflate mail __lambda_func explode require_once block_bot chr set_time_limit strrev strtr Topic #3: base64_encode eval ini_set error_reporting stream_context_create implode strpos md5 base64_decode file_get_contents gzinflate mail preg_match explode fwrite getfilename fileperms filesize getfakedownloadpath getextension Topic #4: strlen __lambda_func str_repeat ceil str_replace create_function base64_decode xor_data__mut rtrim ltrim str_split trim mt_rand assert getcwd getextension getenv getdir getfakedownloadpath getfilename Topic #5: file_get_contents base64_decode eval mail gzinflate explode date_default_timezone_set strrev set_time_limit trim assert strtr getcwd getdir ymsumrlfso get_magic_quotes_gpc getenv getextension getfakedownloadpath getfilename Topic #6: ord chr strlen array_map eval preg_replace base64_decode range krzbnlqhvt getcwd getextension getenv getdoccomment getfakedownloadpath getfilename gethostbyname getdir get_current_user get_magic_quotes_gpc get_info Topic #7: define eval gzinflate __lambda_func set_time_limit get_magic_quotes_gpc islogin md5 header base64_decode mb_ereg_replace loginpage ob_start error_reporting dirname file_exists isref gethostbyname set_magic_quotes_runtime version_compare Topic #8: substr ltrim trim rtrim sizeof qfyuecrkyr ymsumrlfso str_replace str_split assert md5 stripos strlen wzgygjawpn __lambda_func error_reporting array_flip preg_replace strrpos get_info Topic #9: ini_set strlen set_time_limit unserialize xor_data__mut set_magic_quotes_runtime session_start printlogin md5 base64_decode wsologin implode get_magic_quotes_gpc assert clearstatcache defined ignore_user_abort error_reporting eval getdir () Fitting the NMF model (generalized Kullback-Leibler divergence) with tf-idf features, n_samples=2000 and n_features=1000... done in 0.111s. Topics in NMF model (generalized Kullback-Leibler divergence): Topic #0: eval error_reporting stristr preg_match defined gzuncompress base64_decode define require require_once headers_sent ini_get readdir opendir basename function_exists include header time create_function Topic #1: stripos gzuncompress setcookie time require_once set_error_handler dirname file_exists str_ireplace array_map version_compare unlink chdir gmdate ymsumrlfso get_magic_quotes_gpc getcwd getdir getdoccomment getenv Topic #2: create_function function_exists header gzdecode defined ob_start preg_match base64_decode __lambda_func chdir error_reporting pack xor_data__mut getfilename getfakedownloadpath gethostbyname getextension getenv getdoccomment ymsumrlfso Topic #3: stream_context_create base64_encode file_get_contents ini_set strpos md5 implode error_reporting function_exists eval gethostbyname getfilename getcwd getfakedownloadpath getinfo getextension getenv getdoccomment getdir ymsumrlfso Topic #4: str_replace strlen ceil str_repeat base64_decode __lambda_func substr string_n_random stream_context_create strrpos explode file_get_contents reads file_read mt_rand session_start get_magic_quotes_gpc version_compare md5 extension_loaded Topic #5: base64_decode file_get_contents strrev gzinflate mail strtr explode str_rot13 header set_time_limit date_default_timezone_set php_uname eval get_current_user __construct in_array htmlspecialchars html_n file_exists extension_loaded Topic #6: strlen unserialize ord filesize chr fclose is_file gzinflate base64_decode fopen ini_set range set_time_limit xor_data__mut fread mt_rand strstr function_exists explode readdir Topic #7: define strtolower mb_ereg_replace set_time_limit require urldecode get_magic_quotes_gpc md5 dirname strpos file_exists explode require_once ob_start islogin include strdir rand header ini_get Topic #8: substr trim assert str_split rtrim preg_replace sizeof array_flip wzgygjawpn ltrim md5 getwritingscriptpath gzcompress removeoldlogs isbot ymsumrlfso gmdate strrpos var_dump getfakedownloadpath Topic #9: ini_set set_time_limit strpos implode set_magic_quotes_runtime block_bot is_array class_exists unserialize ignore_user_abort wsologin xor_data__mut dirname microtime preg_match str_ireplace printlogin file_put_contents session_start clearstatcache () Fitting LDA models with tf features, n_samples=2000 and n_features=1000... done in 0.553s. Topics in LDA model: Topic #0: ymsumrlfso substr sizeof mt_rand string_n_random preg_replace wzgygjawpn chr eval strstr function_exists str_replace explode strtolower dirname reads base64_encode require_once file_read file_get_contents Topic #1: strlen ord stripos base64_decode eval gzuncompress xor_data__mut preg_replace __lambda_func str_replace chr create_function ltrim trim rtrim str_split ceil str_repeat substr ini_set Topic #2: mb_ereg_replace urldecode assert explode array_map loginpage array_flip print_r chr ini_get fopen ord eval is_callable in_array str_replace strrpos set_magic_quotes_runtime var_dump substr Topic #3: time is_dir filemtime fileowner fileperms filegroup posix_getpwuid is_link posix_getgrgid filesize basename str_replace urlencode realpath function_exists round define htmlspecialchars readdir is_file Topic #4: eval base64_encode ini_set error_reporting file_get_contents md5 base64_decode implode preg_match strpos stream_context_create strdir set_time_limit chop define str_replace getinfo set_magic_quotes_runtime session_start unserialize Topic #5: array_map ord chr krzbnlqhvt substr sizeof strlen xor_data__mut preg_replace base64_decode eval function_exists error_reporting str_replace str_split __lambda_func explode ini_set create_function unserialize Topic #6: makename trydir getfilename exists file_exists getdir call_user_func removeoldlogs remove md5 substr file_put_contents write gmdate strtotime read dds_debug round str_replace __construct Topic #7: is_file readdir filesize round preg_match strrev str_rot13 strtr basename is_dir posix_getgrgid posix_getpwuid is_link fileperms filegroup fileowner filemtime base64_decode count str_replace Topic #8: qfyuecrkyr strpos substr sizeof header ord define isref checkrequest strlen trim _runmagicquotes preg_replace chr base64_decode run clientip is_array function_exists isallowdip Topic #9: eval wzgygjawpn base64_decode stristr gzinflate function_exists preg_match create_function stripos error_reporting define header defined substr ob_start strtolower headers_sent gzdecode __lambda_func gzuncompress ()
我們列印了NMF提取出的10個獨立主題其對應的top權重詞，可以看到，例如eval、base64_decode這種詞是經常出現的詞。這也符合我們的先驗預期，畢竟如果一個文件中經常出現這些詞，則有很大可能說明該文件是一個黑文件。
0x2：基於ICA/NMF濾波器的異常登錄流量檢測
我們知道，一個時間區間內伺服器的登錄流水（ssh/rdp）有可能同時包含有惡意攻擊者和合法管理員的登錄行為，為了降噪的目的，我們可以在登錄流水日誌進入檢測模組之前，進行ICA/NMF因式分解，從原始流量中先分解出攻擊流量和正常流量（類似文章前面提到的MNIST手寫數字降噪），之後對分別對攻擊流量和正常流量的檢測就會變得更容易。
Relevant Link:
https://scikit-learn.org/stable/modules/decomposition.html#nmf https://scikit-learn.org/stable/auto_examples/applications/plot_topics_extraction_with_nmf_lda.html#sphx-glr-auto-examples-applications-plot-topics-extraction-with-nmf-lda-py
