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