文献学习(part101)--CONVEX BICLUSTERING

文献学习(part101)--CONVEX BICLUSTERING

文章目录

​​CONVEX BICLUSTERING​​

​​摘要​​​​Introduction​​​​双聚类的凸形式​​

​​Splitting Methods for Convex Clustering​​

​​凸公式和COBRA解的性质(略)​​​​COBRA法双簇聚类估计​​

​​Splitting Methods for Convex Clustering​​​​cvxclustr(Validation​​​​解决Hold-Out问题​​

​​仿真实验​​

CONVEX BICLUSTERING

摘要

在双聚类问题中，我们寻求同时对观察和特征进行分组. 虽然聚类在从文本挖掘到协同过滤等领域都有应用，但高维基因组数据中识别结构的问题激发了这项工作.

在这种情况下，双聚类使我们能够识别仅在实验条件子集内共同表达的基因子集. 我们给出了双聚类问题的凸公式(目标函数为凸)，它具有唯一的全局最小值和一个保证收敛的迭代算法​​COBRA​​.

我们的方法在单个参数发生变化时生成可能的双聚类的完整解决方案路径. 我们还展示了如何将选择这个调整参数的问题简化为解决凸双聚类问题中一个微不足道的优化.

我们工作的关键贡献是：

它的简单性、可解释性和算法收敛保证——这些特性可以说是目前替代算法所缺乏的.我们展示了我们的方法的优势，包括在模拟和真实的微阵列数据上稳定和可重复地识别双簇.

Introduction

在双聚类问题中，我们试图同时对数据矩阵中的观察值(列)和特征(行)进行分组. 这种数据有时被描述为双向的，或可转换的，以将行和列放在平等的基础上，并强调在行和列变量中揭示结构的愿望.

双聚类用于广泛领域的可视化和探索性分析。例如，在文本挖掘中，双聚类可以识别相对于词的子组具有相似属性的文档的子组(Dhillon，2001)。

在协同过滤中，它可以用于识别对产品子集有相似偏好的客户子群(霍夫曼和普济查，1999)。

关于双聚簇方法的综述可在(Madeira and Oliveira, 2004;Tanay, Sharan and Shamir, 2005;Busygin, Prokopyev和Pardalos, 2008)中找到.

在这项工作中，我们专注于双聚类来识别高维癌症基因组数据中的模式. 虽然癌症，如乳腺癌，临床上可能表现为同质性疾病，但它通常在分子水平上由几个不同的亚型组成. 癌症研究的一个基本目标是确定具有类似分子原形的癌性肿瘤亚型，以及确定每个亚型的特征基因. 识别这些模式是针对患者特定癌症亚型制定个性化治疗策略的第一步.

亚型发现可以作为一个双聚类问题提出，其中基因表达数据被划分成一个棋盘状的模式， 突出患者组和区分这些患者的基因组之间的联系. 双聚类在subtype discovery方面取得了显著的成功. 例如,例如，双聚类乳腺癌数据已经确定了一些基因组，这些基因组的表达水平将患者分为五个亚型，具有不同的生存结果(Sørlie et al , 2001). 这些亚型已在许多研究中重现(Sø rlie等人，2003年). 受这些结果的鼓舞，科学家已经在其他癌症中寻找分子亚型，如卵巢癌(Tothill et al.， 2008).

不幸的是，这些额外研究中的许多发现并不像瑟利等人(2001年)确定的那样具有可复制性. 未能复制这些其他结果可能反映了真正缺乏有生物学意义的分组. 但另一种可能性可能与目前用于识别双簇的计算方法固有的问题有关.

虽然许多双聚簇基因组数据的方法已被提出(Madeira和Oliveira, 2004;Busygin, Prokopyev和Pardalos, 2008)，最流行的癌症基因组数据聚簇的方法是聚簇树状图.

这种方法在病人(列)和基因(行)上执行层次聚类(Hastie, Tibshirani和Friedman, 2009)。

然后根据这些单独的行和列树状图重新排序矩阵，并将其可视化为沿着行和列轴绘制**树状图的热图。**图2显示了肺癌研究中表达数据的聚类树状图示例：

该数据包括56个个体中500个基因的表达水平，这是Lee等人(2010年)研究的数据子集. 受试者分为四个亚组:正常、类癌、结肠或小细胞. 这一简单的策略似乎可以合理地恢复临床诊断，并识别出以这些亚型为特征的失调基因.

然而，作为一种算法，聚类树图有两个特点，使其不太适合产生可再现的结果:

树状图是通过贪婪地融合观测值(特征)来降低某些标准而构建的. 因此，该算法可以返回相对于该标准仅局部最优的双聚类. 由于解决方案可能会因算法的初始化方式而异，因此这类过程往往会在多次初始化后运行，但即使如此，也不能保证会达到全局最优.该算法也不稳定，因为数据中的小扰动会导致簇分配的大变化.

更复杂的聚类方法已经被提出，一些基于数据矩阵的奇异值分解(SVD) ，以及其他基于graph cuts的方法.

有些方法在本质上类似于聚类树状图，直接对行和列进行聚类.

虽然这些方法可能在经验性能方面提供有价值的改进，但它们都没有解决困扰聚类树图的两个基本问题.

此外，科学家可能不愿使用这些方法，因为与简单的聚类谱系图相比，它们的输出通常不太容易可视化. 从可再现的研究角度来看，双聚类方法应该给出相同的，理想情况下唯一的答案，而不管算法是如何初始化的，并且对于数据中的扰动保持稳定.

本文将双聚类作为一个凸优化问题，提出了一种新的凸双聚类算法(​​COBRA​​)来迭代求解. COBRA输出的结果保留了聚类树图的简单可解释性和可视化，并且与现有技术相比还具有几个关键优势：

(1)稳定性和唯一性:COBRA为凸规划产生唯一的全局最小化，该最小化在数据中是连续的. 这意味着COBRA总是将数据映射到单个双聚类分配，并且这个解决方案是稳定的.(2)简单性:COBRA采用单个调谐参数(a single tuning paramete)来控制双簇的数量.(3) 数据自适应性:COBRA admits a simple and principled dataadaptiveprocedure for choosing the tuning parameter that involves solving a convex matrix completion problem.

回到我们激发肺癌的例子，图2说明了COBRA的结果，其调整参数是根据我们的数据自适应过程选择的.

在通过系列化包(Hahsler，Hornik和Buchta，2008)对列和行进行重新排序后，通过“TSP”选项，出现了更清晰的双聚类模式。

双聚类的凸形式

我们的目标是识别数据矩阵中相互关联的行组和列组. 如图2所示，当行和列根据它们的分组被重新排序时，棋盘模式出现，即由行和列组定义的矩阵分区的元素趋向于显示均匀的强度.

这个双聚簇模型相当于一个棋盘平均模型(Madeira和Oliveira, 2004)。这个模型与Tan和Witten(2014)的假设最为相似，他们提出了一些方法来估计一些双群集平均条目是稀疏的棋盘状结构.

棋盘模型是穷尽的，因为每个矩阵元素被分配给一个给定双簇. 这与其他双聚类模型不同，后者识别出可能重叠的行和列子集，但不是穷尽的，这些通常使用svd类方法进行估计.

相比之下，估计行和列分区是一个组合难问题，也是双聚类的主要目标.

这个任务类似于回归问题中的最佳子集选择(Hastie, Tibshirani和Friedman, 2009). 然而，正如最佳子集选择问题已经成功地通过解决一个凸代理问题，即Lasso来解决，我们将设计一个凸松弛的组合困难的问题，选择行和列划分.

当只有行或列被聚类时，(1)中的目标函数进行最小化，就相当于在L2范数融合惩罚下解决一个凸聚类问题。

反过来，凸聚类问题可以看作是融合Lasso的推广。

The convex clustering problem in turn can be seen as a generalization of the Fused Lasso.

在本文中，我们将自己限制在L2范数，因为它是旋转不变的。一般来说，当数据的坐标表示发生微小变化时，我们不希望程序的双聚类输出发生微小变化。

如图2所示，受惩罚的估计显示出期望的棋盘结构。

请注意，这种收缩过程从根本上不同于聚类树图这样的方法，后者独立地对行和列进行聚类。

通过耦合行和列聚类，我们的公式明确地寻求具有棋盘平均结构的解决方案。

我们现在解决选择权重的问题。对权重的明智选择使我们能够：

考虑到这些目标，我们建议权重具有以下属性：

我们现在讨论关于权重建议背后的基本原理：

参考文献

Splitting Methods for Convex Clustering

如前所述，将正权重限制到最近的邻居可以提高计算效率和聚类质量。虽然定义权重的两个因素作用相似，但是它们的组合增加了聚类路径对数据的局部密度的敏感性。

凸公式和COBRA解的性质(略)

COBRA法双簇聚类估计

另一种常用的使非光滑凸函数最小化的迭代方法是交替方向乘子法(ADMM)。

虽然ADMM算法是可行的，但我们采用了Bauschke和Combettes提出的更直接的Dykstra-like proximal算法(DLPA)近似算法计算凸双聚类路径，因为它产生了一个简单的元算法，可以利用快速求解凸聚类问题。

DLPA算法

算法1中的每个邻近映射对应于解决凸聚类问题。

命题1

本文利用Chi和Lange (2015)提出的交替最小化算法(AMA)来解决凸聚类问题。该算法对拉格朗日对偶问题进行投影梯度上升。

它的主要优点是，当我们使用所描述的稀疏高斯核权重时，它需要计算量和与数据矩阵X的大小成线性关系。

第二个优点是，硬聚类分配很容易从分裂方法中使用的变量中获得。

尽管如此，DLPA框架使得交换更有效的解算器变得轻而易举。

参考文献

Splitting Methods for Convex Clustering

GitHub

cvxclustr(Proximal mapping: h(x) = tau * ||x||_2temp = function(x, tau) { lv = norm(as.matrix(x),type="f") if (lv == 0) { px = x } else { px = max(0, 1-tau/lv)*x } return(px)}prox_L2 = cmpfun(temp)## Proximal mapping: h(x) = tau * ||x||_1temp = function(x,tau) { return(sign(x)*sapply(abs(x)-tau,FUN=function(x) {max(x,0)}))}prox_L1 = cmpfun(temp)## Projection onto L2 ball of radius tautemp = function(x,tau) { lv = norm(as.matrix(x),'f') u = x if (lv > tau) { u = (tau/lv)*x } return(u)}proj_L2 = cmpfun(temp)## Projection onto L-infinty ball of radius tautemp = function(x,tau) { u = as.matrix(abs(x)) cbind(u,tau) u = apply(cbind(u,tau),1,min) u = u*sign(x) return(u)}proj_Linf = cmpfun(temp)## Projection onto a simplex of size ztemp = function(v,z) { mu = sort(v, decreasing=T) partial_sum <- mu[1] rho = 1 for (j in 2:length(v)) { partial_sum <- partial_sum + mu[j] if (j*mu[j] - partial_sum + z <= 0) { rho = j - 1 break } rho = j } theta <- mean(mu[1:rho]) - z/rho w <- v - theta w[w <= 0] = 0 return(w)}project_to_simplex = cmpfun(temp)## Projection onto L1 ball of radius tautemp = function(x,tau) { w = project_to_simplex(abs(x),tau) return(w*sign(x))}proj_L1 = cmpfun(temp)temp = function(x,tau,type) { if (type == 1) { return(proj_L2(x,tau)) } else if (type==2) { return(prox_L1(x,tau)) } else if (type==3) { return(prox_L2(x,tau)) } else if (type==4) { return(proj_Linf(x,tau)) } else if (type==5) { return(proj_L1(x,tau)) }}prox_R = cmpfun(temp)

cluster_path_preprocess.R

## Clusterpath preprocessinglibrary(compiler)temp = function(i,j,p) { return(p*(i-1) - i*(i-1)/2 + j -i)}tri2vec = cmpfun(temp)temp = function(k,p) { i = ceiling(0.5*(2*p-1 - sqrt((2*p-1)^2 - 8*k))) j = k - p*(i-1) + i*(i-1)/2 + i return(as.matrix(cbind(i,j)))}vec2tri = cmpfun(temp)temp = function(w,p) { sizes1 = double(p) sizes2 = double(p) P = vec2tri(which(w > 0),p) M1 = matrix(0,nrow(P),p) M2 = matrix(0,nrow(P),p) for (i in 1:p) { group1 = which(P[,1] == i) sizes1[i] = length(group1) if (sizes1[i] > 0) { M1[1:sizes1[i],i] = group1 } group2 = which(P[,2] == i) sizes2[i] = length(group2) if (sizes2[i] > 0) { M2[1:sizes2[i],i] = group2 } } M1 = M1[1:max(sizes1),,drop=FALSE] M2 = M2[1:max(sizes2),,drop=FALSE] return(list(ix=P,M1=M1,M2=M2,s1=sizes1,s2=sizes2))}compactify_edges = cmpfun(temp)temp = function(X,mu) { p = ncol(X) k = 1 w = matrix(0,p*(p-1)/2,1) for (i in 1:(p-1)) { for (j in (i+1):p) { w[k] = exp(-mu*norm(as.matrix(X[,i,drop=FALSE]-X[,j,drop=FALSE]),'f')^2) k = k+1 } } return(weights=w)}kernel_weights = cmpfun(temp)temp = function(w,k,p) { i = 1 neighbors = tri2vec(i,(i+1):p,p) keep = neighbors[sort(w[neighbors],decreasing=TRUE,index.return=TRUE)$ix[1:k]] for (i in 2:(p-1)) { group_A = tri2vec(i,(i+1):p,p) group_B = tri2vec(1:(i-1),i,p) neighbors = c(group_A,group_B) knn = neighbors[sort(w[neighbors],decreasing=TRUE,index.return=TRUE)$ix[1:k]] keep = union(knn,keep) } i = p neighbors = tri2vec(1:(i-1),i,p) knn = neighbors[sort(w[neighbors],decreasing=TRUE,index.return=TRUE)$ix[1:k]] keep = union(knn,keep) w[-keep] = 0 return(w)}knn_weights = cmpfun(temp)rm(temp)

ama_loss.R

temp = function(X,U,gamma,w,ix) { data_fidelity = 0.5*(norm(as.matrix(X-U),'f')^2) center_distances = U[,ix[,1],drop=FALSE] - U[,ix[,2],drop=FALSE] penalty = sum(w*apply(center_distances,2,FUN=function(x) {norm(as.matrix(x),'2')})) return(data_fidelity + gamma * penalty)}loss_primalR = cmpfun(temp)temp = function(X,Lambda,ix,M1,M2,s1,s2) { q = nrow(X) p = ncol(X) nK = ncol(Lambda) L = matrix(0,q,p) for (i in 1:p) { if (s1[i] > 0) { jx = M1[1:s1[i],i] L[,i] = apply(Lambda[,jx,drop=FALSE],1,sum) } if (s2[i] > 0) { jx = M2[1:s2[i],i] L[,i] = L[,i] - apply(Lambda[,jx,drop=FALSE],1,sum) } } first_term = 0.5*norm(L,'f')^2 second_term = sum((X[,ix[,1]] - X[,ix[,2]])*Lambda) return(-first_term-second_term)}loss_dualR = cmpfun(temp)rm(temp)

ama_updates.R

temp = function(X,Lambda,M1,M2,s1,s2) { p = ncol(X) U = X for (i in 1:p) { if (s1[i] > 0) { ix = M1[1:s1[i],i] U[,i] = U[,i] + apply(Lambda[,ix,drop=FALSE],1,sum) } if (s2[i] > 0) { ix = M2[1:s2[i],i] U[,i] = U[,i] - apply(Lambda[,ix,drop=FALSE],1,sum) } } return(U)}update_UR = cmpfun(temp) temp = function(U,Lambda,w,gamma,nu,ix) { q = nrow(U) nK = ncol(Lambda) V = matrix(0,q,nK) for (kk in 1:nK) { i = ix[kk,1] j = ix[kk,2] z = U[,i] - U[,j] - (1/nu)*Lambda[,kk] V[,kk] = prox_L2(z, w[kk]*gamma/nu) } return(V)}update_VR = cmpfun(temp)temp = function(Lambda,U,V,nu,ix) { Z = U[,ix[,1],drop=FALSE] - U[,ix[,2],drop=FALSE] - V return(Lambda - nu*Z)}update_LambdaR = cmpfun(temp)temp = function(Lambda,U,w,gamma,nu,ix) { V = update_VR(U,Lambda,w,gamma,nu,ix) return(update_LambdaR(Lambda,U,V,nu,ix))}update_Lambda_combinedR = cmpfun(temp)rm(temp)

算法1 凸双聚类算法(COBRA)

设U0 = XP0 = 0矩阵Q0 = 0矩阵while (收敛条件): Ym = prox(Um.T + pm.T) ## 行凸聚类 Pm_p1 = Um + Pm -Ym.T Um_p1 = prox(Ym.T + Qm.T) ## 列凸聚类 Qm_p1 = Ym + Qm - Um_p1.T

模型的选择

估计数据集中的簇或双簇的数量是一项重大挑战，对于许多现有的双聚类方法：

必须选择的许多调整参数；以及双聚类分配并不总是随着双聚类的数量而平滑变化。

例如，稀疏奇异值分解方法(Lee等人，2010)需要三个调谐参数:两个参数控制稀疏奇异值分解的左右奇异向量的稀疏性，一个参数控制其秩。

此外，对于大型数据集，选择双聚类的数量和其他调整参数可能是一个主要的计算负担。例如，稀疏双聚类方法(Tan和Witten，2014)使用交叉验证来选择行和列分区的数量。

Hold-Out Validation

这个想法最初是由Wold (1978)提出的，用于主成分分析中的模型选择，最近被用于选择矩阵完成问题的调整参数(Mazumder，Hastie和Tibshirani，2010)。

解决Hold-Out问题

在(4)中定义的问题可以看作是一个凸矩阵补全问题。

算法2 缺少数据的COBRA

初始化U0while (收敛条件): M = P_theta_c(X) + P_theta(Um) Um_p1 = COBRA(M)

这种相似性并不是巧合，因为这两个过程都是优化-最小化(MM)算法(Lange，Hunter and Yang，2000)的实例，它们在光滑二次项上应用相同的优化。

算法2有以下收敛保证.

命题2

仿真实验

我们比较了COBRA和另外两种同样假定棋盘平均结构的双聚类方法。

第一个是聚类树状图;使用在dynamicTreeCut包中实现的广泛使用的动态树切割算法(Langfelder, Zhang和Horvath, 2008)，对聚类树形图进行了复杂的双聚类分配。

第二种是稀疏双聚簇方法(Tan and Witten, 2014)，在sparseBC包中实现。所有参数均根据R包中的方法选取。

我们使用以下三种衡量标准:

兰德指数(Rand index, RI)调整后的兰德指数(ARI)信息变化(VI)

实验结果：

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。