Majorization Minimization by Coordinate Descent for Concave Penalized Generalized Linear Models

Dingfeng Jiang; Jian Huang

doi:10.1007/s11222-013-9407-3

统计计算。作者手稿；PMC 2014年10月9日提供。

以最终编辑形式发布为：

统计计算。2014年9月；24(5): 871–883.

在线发布2013年6月6日。数字对象标识：2007年10月17日/11222-013-9407-3

预防性维修识别码：项目经理4191872

NIHMSID公司：NIHMS532755号

PMID：25309048

凹惩罚广义线性模型的坐标下降优化最小化

丁凤江和黄健

作者信息版权和许可信息 PMC免责声明

摘要

最近的研究已经证明了一类凹罚函数在变量选择中的理论吸引力，包括平滑剪裁的绝对偏差和极小极大凹罚函数。然而，计算高维模型中的凹惩罚解是一项困难的任务。我们提出了一种坐标下降优化最小化（MMCD）算法来计算广义线性模型中的凹惩罚解。与使用惩罚函数的局部二次或局部线性近似的现有算法相比，市场管理与分销部寻求通过二次损失来优化负对数似然，但不使用惩罚的任何近似。该策略可以避免每次更新解时计算比例因子，从而提高坐标下降的效率。在一定的正则性条件下，我们建立了MMCD的理论收敛性。我们使用SCAD和MCP惩罚对惩罚logistic回归模型实现了该算法。模拟研究和数据示例表明，在协变量数量远大于样本量的高维环境中，MMCD对于惩罚逻辑回归的工作速度足够快。

关键词：逻辑回归，第页≫n个模型，平滑剪裁的绝对偏差惩罚，minimax凹惩罚，变量选择

1引言

变量选择是统计学中的一个基本问题。惩罚方法在变量选择方面具有吸引人的理论性质第页≫n个型号，带n个样本量和第页变量的数量。已经提出了几项重要处罚措施。示例包括ℓ₁惩罚或最小绝对收缩和选择运算符（Lasso）(Tibshirani（1996）)，平滑剪裁绝对偏差（SCAD）惩罚(范和李（2001）)和最小最大凹罚（MCP）(张（2010）). SCAD和MCP是凹惩罚，具有预言性质，这意味着它们可以正确地选择重要变量，并以高概率估计变量系数，就好像在某些稀疏性条件和其他适当的正则性条件下提前知道模型一样。

惩罚回归模型的计算算法已经取得了相当大的进展。埃夫隆等(2004)结果表明，改进后的LARS算法可以在线性模型中有效计算整个Lasso解路径。此修改的LARS算法与之前由奥斯本、普雷斯内尔和图拉赫（2000）对于SCAD罚款，范和李（2001）提出了一种局部二次逼近（LQA）算法。LQA算法的一个缺点是，一旦系数在任何迭代步骤设置为零，它将永久保持为零，并且相应的变量将从最终模型中删除。亨特和李（2005）使用优化最小化（MM）方法，通过将分母从零上界来优化LQA的扰动版本。如何选择扰动的大小以及扰动如何影响稀疏性需要在特定模型中确定。邹和李（2008）提出了一种局部线性逼近（LLA）算法。LLA算法通过重复使用拉索惩罚算法来近似凹惩罚解。Schifano、Strawderman和Wells（2010年）通过MM方法将LLA算法的思想推广到多重惩罚，并通过迭代软阈值（MIST）算法证明了其最小化的收敛性。张（2010）开发了PLUS算法，用于计算线性回归模型中的凹惩罚最小二乘解，包括MCP解。

在过去的几年中，人们已经认识到坐标下降算法（CDA）可以有效地计算拉索解第页≫n个模型。该算法在应用数学中有着悠久的历史，其根源在于求解线性系统的高斯-西德尔方法(Warge（1963年）;Ortega和Rheinbold（1970）;曾（2001）). CDA通过一次处理一个坐标（或一个坐标块）来优化目标函数，在所有坐标中迭代循环，直到达到收敛。它特别适用于每个坐标都有简单解决方案，但在更高维度上缺乏解决方案的问题。拉索惩罚线性模型的CDA与LARS具有很强的竞争力，特别是在高维情况下(Friedman、Hastie、Höfling和Tibshirani（2007）;Wu和Lange（2008）;Friedman、Hastie和Tibshirani（2010）).

坐标下降法也用于计算凹面惩罚解路径。Breheny和Huang（2011）比较了不同组合的CDA和LLA(n、第页)以及协变量矩阵的各种设计。他们的结果表明，在他们考虑的设置下，CDA比LLA-LARS算法收敛得更快。Mazumder、Friedman和Hastie（2011）证明了CDA比LLA具有更好的收敛性。Breheny和Huang（2011）为了克服MCP惩罚广义线性模型（GLM）解的计算中尺度因子不断变化的困难，提出了一种自适应尺度调整技术。然而，自适应缩放方法不能应用于SCAD惩罚。此外，尚不清楚其解是否达到原始目标函数的局部最优点。

我们提出了一种基于坐标下降的优化最小化（MMCD）算法来计算凹惩罚GLM模型的解，重点是逻辑回归。MMCD为每个坐标寻求一个闭合形式的解，并通过优化损失函数来避免计算缩放因子。在合理的正则性条件下，我们建立了MMCD的收敛性。该算法特别适用于逻辑回归模型，因为可以找到简单有效的优化。本文的结构如下。第2节定义了GLM中的凹惩罚解。第3节描述了所提出的MMCD算法，解释了优化的优点并研究了其收敛性。第四节在凹惩罚logistic模型中实现了MMCD算法。通过仿真研究，比较了MMCD算法及其竞争对手的计算效率和选择性能。第5节给出了结束语。

2 GLM的凹面惩罚解决方案

让 ${{(年_{我}, {x个}_{我})}_{我 = 1}^{n个}}$ 为观测数据，其中年_我是响应变量x个_我是一个(第页+1）-预测器的维向量，第一个元素为1，对应于截距。我们考虑具有标准链接函数的GLM，其中年_我与…有关x个_我通过线性组合 $η_{我} = {x个}_{我}^{T型} β$ ，β=（β₀, β₁, …, β_第页)^T型∈ ℝ^第页+1.此处β₀是截距。的条件密度函数年_我鉴于x个_我是（f）_我(年_我)=经验{(年_我θ_我− ψ(θ_我))/ϕ_我+c（c）(年_我, ϕ)}. 此处_我>0是色散参数。ψ（θ）的形式取决于指定的模型。例如，逻辑模型中的ψ（θ）=log（1+exp（θ））。（标度）负对数似然函数为

ℓ (β) \propto \frac{1}{n个} \sum_{我 = 1}^{n个} {ψ ({x个}_{我}^{T型} β) - 年_{我} {x个}_{我}^{T型} β} .

（1）

在这里x个_我0= 1, 1 ≤我≤n个.对于另一个第页变量，我们假设它们是标准化的，即， ${‖ {x个}^{j个} ‖}_{2}^{2} / n个 = 1$ 具有x个^j个= (x个_1j个, …,x个_新泽西州)^T型, 1 ≤j个≤第页.在这里v（v）‖₂是ℓ₂a的范数n个-量纲向量v（v）标准化允许对每个变量均匀应用惩罚。

我们考虑凹惩罚GLM准则

问 (β; λ, γ) = ℓ (β) + \sum_{j个 = 1}^{第页} ρ (| β_{j个} |; λ, γ),

(2)

其中，ρ是惩罚参数λ≥0和控制ρ形状的正则化参数γ的惩罚函数。请注意，在(2)截距β₀不受处罚。我们重点关注两个凹面点球，SCAD和MCP。SCAD惩罚(范和李（2001）)定义为

ρ (t吨; λ, γ) = λ {¦Β}_{0}^{| t吨 |} 1_{{x个 \leq λ}} + \frac{{(γ λ - x个)}_{+}}{(γ - 1) λ} 1_{{x个 > λ}} d日 x个,

(3)

λ≥0且γ>2。在这里1_x个∈A类指示灯功能和x个₊=x个1_{x个≥0｝表示的非负部分x个MCP罚款(张（2010）)定义为

ρ (t吨; λ, γ) = λ {¦Β}_{0}^{| t吨 |} {(1 - \frac{x个}{γ λ})}_{+} d日 x个

（4）

λ≥0且γ>1。

对于SCAD和MCP惩罚，正则化参数γ控制凹度，较小的γ对应更凹的惩罚。这两种惩罚的开始是采用与拉索相同的惩罚程度，然后逐渐将惩罚减至零|t吨|增加。当γ→ ∞, SCAD和MCP惩罚集中于ℓ₁处罚。

为了对这些惩罚有一个基本的理解，考虑一个定义为惩罚一元线性回归解的阈值算子，

θ̂ (λ, γ) = \underset{θ}{argmin（最小值）} {\frac{1}{2 n个} \sum_{我 = 1}^{n个} {(年_{我} - {x个}_{我} θ)}^{2} + ρ (θ; λ, γ)} .

让 ${θ̂}_{L（左） S公司} = \sum_{我 = 1}^{n个} {x个}_{我} 年_{我} / \sum_{我 = 1}^{n个} {x个}_{我}^{2}$ 是最小二乘解。通过表示软阈值运算符S公司(t吨，λ）=符号(t吨) (|t吨| − λ)₊对于λ≥0(多诺霍和约翰斯通（1994）)，何处标志(t吨)=1,0，-1，如果t吨>分别为0、=0、<0。然后，SCAD和MCP处罚针对 $\hat{θ}$ （λ，γ）如下所示，

{θ̂}_{S公司 C类 A类 D类} (λ, γ) = S公司 (θ_{L（左） S公司}, λ) 1_{{| {θ̂}_{L（左） S公司} | \leq 2 λ}} + \frac{γ - 1}{γ - 2} S公司 ({θ̂}_{L（左） S公司}, λ γ / (γ - 1)) 1_{{2 λ < | {θ̂}_{L（左） S公司} | \leq γ λ}} + {θ̂}_{L（左） S公司} 1_{{| {θ̂}_{L（左） S公司} | > λ γ}},

(5)

{θ̂}_{M（M） C类 P（P）} (λ, γ) = \frac{γ}{γ - 1} S公司 ({θ̂}_{L（左） S公司}, λ) 1_{{| {θ̂}_{L（左） S公司} | \leq λ γ}} + {θ̂}_{L（左） S公司} 1_{{| {θ̂}_{L（左） S公司} | > λ γ}} .

(6)

观察SCAD和MCP是否使用LS溶液，如果| $\hat{θ}$ _{LS（负载感应）}|>λγ；MCP仅在以下情况下应用缩放的软阈值操作| $\hat{θ}$ _{LS（负载感应）}|≤λγ，而SCAD应用软阈值操作，如果| $\hat{θ}$ _{LS（负载感应）}|<2λ，如果2λ<| $\hat{θ}$ _{LS（负载感应）}| ≤ λγ. 这些阈值化算子是下文所述MMCD算法的基本构建块。

3坐标下降的最大化最小化

对于GLM，二次近似值为ℓ（β）在给定估计的附近 $\tilde{β}$ 导致加权最小二乘损失，

ℓ^{秒} (β | β̃) = \frac{1}{2 n个} \sum_{我 = 1}^{n个} {w个}_{我} {({z（z）}_{我} - {x个}_{我}^{T型} β)}^{2},

(7)

具有 ${w个}_{我} (β̃) = ψ̈ ({x个}_{我}^{T型} β)$ 和 ${z（z）}_{我} (β̃) = ψ̈ {({x个}_{我}^{T型} β̃)}^{- 1} {年_{我} - ψ̇ ({x个}_{我}^{T型} β̃)} + {x个}_{我}^{T型} β̃$ 式中ψ̇（θ）和 $\ddot{ψ}$ （θ）是ψ（θ）相对于（w.r.t.）θ的一阶和二阶导数。使用ℓ^秒(β| $\tilde{β}$ )在criteria函数中，CDA更新j个通过固定剩余的坐标k个(k个≠j个)坐标。让 $β_{j个}^{米} = {({β̂}_{0}^{米 + 1}, \dots, {β̂}_{j个}^{米 + 1}, {β̂}_{j个 + 1}^{米}, \dots, {β̂}_{第页}^{米})}^{T型}$ ，CDA更新 ${β̂}_{j个 - 1}^{米}$ 到 ${β̂}_{j个}^{米}$ 通过最小化标准

{β̂}_{j个}^{米 + 1} = \underset{β_{j个}}{argmin（最小值）} 问^{秒} (β_{j个} | {β̂}_{j个 - 1}^{米}) = \underset{β_{j个}}{argmin（最小值）} \frac{1}{2 n个} \sum_{我 = 1}^{n个} {w个}_{我} {({z（z）}_{我} - \sum_{秒 < j个} {x个}_{我 j个} {β̂}_{秒}^{米 + 1} - {x个}_{我 j个} β_{j个} - \sum_{秒 > j个} {x个}_{我 j个} β_{秒}^{米})}^{2} + ρ (| β_{j个} |; λ, γ),

(8)

哪里w个_我和z（z）_我依靠 $({β̂}_{j个 - 1}^{米}, {x个}_{我}, 年_{我})$ . Thej个然后通过求解该方程得到坐标向极小值，

\frac{1}{n个} \sum_{我 = 1}^{n个} {w个}_{我} {x个}_{我 j个}^{2} β_{j个} + ρ' (| β_{j个} |) sgn公司 (β_{j个}) - \frac{1}{n个} \sum_{我 = 1}^{n个} {w个}_{我} {x个}_{我 j个} ({z（z）}_{我} - {x个}_{我}^{T型} {β̂}_{j个 - 1}^{米}) - \frac{1}{n个} \sum_{我 = 1}^{n个} {w个}_{我} {x个}_{我 j个}^{2} {β̂}_{j个}^{米} = 0,

(9)

带有ρ′(|t吨|)ρ的一阶导数(|t吨|)水电阻温度|t吨|的t吨≠ 0. 请注意，对于t吨=0，则ρ′（0+）=λ。符号sgn(t吨)是的次微分|t吨|也就是说，

sgn公司 (t吨) = {\begin{matrix} 1, & 如果 t吨 > 0; \\ - 1, & 如果 t吨 < 0; \\ \in [- 1, 1], & 如果 t吨 = 0 . \end{matrix}

我们指的是兰格（2004）有关子差速器的详细说明。

通过定义比例因子 $δ_{j个} ≜ {n个}^{- 1} \sum_{我 = 1}^{n个} {w个}_{我} {x个}_{我 j个}^{2}$ .直接求解(9)由于MCP处罚j个协调解决方案为，

{β̂}_{j个}^{米 + 1} = \frac{S公司 (τ_{j个}, λ)}{δ_{j个} - 1 / γ} 1_{{| τ_{j个} | \leq δ_{j个} γ λ}} + \frac{τ_{j个}}{δ_{j个}} 1_{{| τ_{j个} | > δ_{j个} γ λ}},

(10)

具有 $τ_{j个} = {n个}^{- 1} \sum_{我 = 1}^{n个} {w个}_{我} {x个}_{我 j个} ({z（z）}_{我} - {x个}_{我}^{T型} {β̂}_{j个 - 1}^{米}) + δ_{j个} {β̂}_{j个}^{米}$ .观察当|τ_j个|≤λ，我们有 ${β̂}_{j个}^{米 + 1} = 0$ 在线性模型中，w个_我= 1,我= 1, …,n个因此，比例因子δ_j个=1（对于标准化预测值）。然而，在GLM中w个_我在 $({β̂}_{j个 - 1}^{米}, {x个}_{我}, 年_{我})$ 导致比例因子δ_j个从一个迭代到另一个迭代。这是有问题的，因为δ_j个−1/γ可能非常小，并且不能保证是正的。因此，直接应用CDA可能不具有数值稳定性，并可能导致不合理的解决方案。

为了克服这个困难，Breheny和Huang（2011）提出了一种自适应缩放方法，该方法使用

{β̂}_{j个}^{米 + 1} = \frac{S公司 (τ_{j个}, λ)}{δ_{j个} (1 - 1 / γ)} 1_{{| τ_{j个} | \leq γ λ}} + \frac{τ_{j个}}{δ_{j个}} 1_{{| τ_{j个} | > γ λ}},

(11)

对于j个协调解决方案。这相当于应用一个新的正则化参数γ*=γ/δ_j个在每个坐标方向的迭代中。因此，对于惩罚变量，有效的正则化参数是不相同的，直到算法收敛时才知道。数值上，比例因子δ_j个需要额外的计算，这对于大型计算机来说是不可取的第页此外，小δ_j个可能会导致趋同问题。自适应重缩放方法不能应用于SCAD惩罚，因为缩放的软阈值操作仅适用于SCAD阈值操作符的中间子句，如所示(5).

观察到在GLM中，比例因子δ_j个等于损失函数的二阶偏导数，即， $\nabla_{j个}^{2} ℓ (β) = \sum ψ̈ ({x个}_{我}^{T型} β) {x个}_{我 j个}^{2} / n个 = \sum {w个}_{我} {x个}_{我 j个}^{2} / n个$ 。MMCD算法寻求的上限为 $\nabla_{j个}^{2} ℓ (β)$ 。因此，我们需要找到一个M（M）这样δ_j个≤M（M）对于每个坐标。在二阶导数上找到一致界的想法也是由Böhning和Lindsay（1988）.

对于MM方法，使用δ的上限_j个等同于查找代理函数 $ℓ^{M（M） M（M）} (β_{j个} | {β̂}_{j个 - 1}^{米})$ 主要是 $ℓ^{秒} (β_{j个} | {β̂}_{j个 - 1}^{米})$ ，其中

ℓ^{M（M） M（M）} (β_{j个} | {β̂}_{j个 - 1}^{米}) = ℓ ({β̂}_{j个 - 1}^{米}) + \nabla_{j个} ℓ ({β̂}_{j个 - 1}^{米}) (β_{j个} - {β̂}_{j个}^{米}) + \frac{1}{2} M（M） {(β_{j个} - {β̂}_{j个}^{米})}^{2},

(12)

和

ℓ^{秒} (β_{j个} | {β̂}_{j个 - 1}^{米}) = ℓ (β_{j个 - 1}^{米}) + \nabla_{j个} ℓ ({β̂}_{j个 - 1}^{米}) (β_{j个} - {β̂}_{j个}^{米}) + \frac{1}{2} \nabla_{j个}^{2} ℓ ({β̂}_{j个 - 1}^{米}) {(β_{j个} - {β̂}_{j个}^{米})}^{2},

(13)

带二阶偏导数 $\nabla_{j个}^{2} ℓ ({β̂}_{j个 - 1}^{米})$ 用泰勒展开式的上界代替M（M）请注意，对于给定的 ${β̂}_{j个 - 1}^{米}, Ş_{j个}^{2} ℓ ({β̂}_{j个 - 1}^{米})$ 是一个常量，但是 $\nabla_{j个}^{2} ℓ ({β̂}_{j个 - 1}^{米})$ 从迭代到迭代的更改。多数化使用上限来避免这些变化。无论当前估计的 ${β̂}_{j个 - 1}^{米}$ 这与MM方法的常见实现不同，其中优化是在当前估计点构造的。还要注意的是，为了更好地适应CDA，协调应用了多数化。MM方法的下降特性确保了 $ℓ^{M（M） M（M）} (β_{j个} | {β̂}_{j个 - 1}^{米})$ 导致原始目标函数的下降序列。有关MM算法的更多详细信息，请参阅兰格、亨特和杨（2000）;亨特和兰格（2004）.

考虑到δ的优化_j个，一些代数表明j个SCAD和MCP的协调解决方案如下

S公司 C类 A类 D类 : {β̂}_{j个}^{米 + 1} = \frac{1}{M（M）} S公司 (τ_{j个}, λ) 1_{{| τ_{j个} | \leq (1 + M（M）) λ}} + \frac{S公司 (τ_{j个}, γ λ / (γ - 1))}{M（M） - 1 / (γ - 1)} 1_{{(1 + M（M）) λ < | τ_{j个} | \leq M（M） γ λ}} + \frac{1}{M（M）} τ_{j个} 1_{{| τ_{j个} | > M（M） γ λ}},

(14)

M（M） C类 P（P） : {β̂}_{j个}^{米 + 1} = \frac{S公司 (τ_{j个}, λ)}{M（M） - 1 / γ} 1_{{| τ_{j个} | \leq M（M） γ λ}} + \frac{1}{M（M）} τ_{j个} 1_{{| τ_{j个} | > M（M） γ λ}},

(15)

具有 $τ_{j个} = M（M） {β̂}_{j个}^{米} + {n个}^{- 1} \sum_{我 = 1}^{n个} {x个}_{我 j个} (年_{我} - ψ̇ ({x个}_{我}^{T型} {β̂}_{j个 - 1}^{米})), j个 = 1, \dots, 第页$ 。截距的解决方案是

{β̂}_{0}^{米 + 1} = τ_{0} / M（M）,

(16)

具有 $τ_{0} = M（M） {β̂}_{0}^{米} + {n个}^{- 1} \sum_{我 = 1}^{n个} {x个}_{我 0} (年_{我} - ψ̇ ({x个}_{我}^{T型} {β̂}^{米}))$ ，其中 ${β̂}^{米} = {({β̂}_{0}^{米}, {β̂}_{1}^{米}, \dots, {β̂}_{第页}^{米})}^{T型}$ .英寸(14)和(15)，我们希望确保分母为正，即。M（M）−1/（γ−1）>0和M（M）− 1/γ > 0. 这自然导致对惩罚的限制，inf_t吨ρ〃(|t吨|;λ, γ) > −M（M），式中ρ〃(|t吨|;λ、 γ）是ρ的二阶导数(|t吨|;λ、 γ）水电阻|t吨|。对于SCAD和MCP，通过选择合适的γ来满足条件。自inf以来_t吨ρ″(|t吨|;λ、对于SCAD惩罚，γ）=−1/（γ−1），以及inf_t吨ρ〃(|t吨|;λ、 γ）=−1/对于MCP，我们需要γ>1+1/M（M）对于SCAD和γ>1/M（M）用于MCP。

通过采用以下技巧，MMCD算法可以进一步提高效率。设η=（η₁, …, η_n个)^T型和 $X（X） = {({x个}_{1}^{T型}, \dots, {x个}_{n个}^{T型})}^{T型}$ 、和 ${η̂}_{j个}^{米} = X（X） {β̂}_{j个}^{米}$ 是对应于 ${β̂}_{j个}^{米}$ 可以通过使用以下等式来实现进一步的效率

{η̂}_{j个 + 1}^{米} = {η̂}_{j个}^{米} + {x个}^{j个 + 1} ({β̂}_{j个 + 1}^{米 + 1} - {β̂}_{j个 + 1}^{米}) = {η̂}_{j个}^{米} + ({β̂}_{j个 + 1}^{米} - {β̂}_{j个}^{米}) {x个}^{j个 + 1} .

(17)

该等式使O（运行）(净现值)操作到O（运行）(n个)一个。由于这一步涉及每个坐标的每次迭代，因此它对减少计算成本非常重要。

我们还可以使用LQA优化每次迭代中的惩罚(Fan和Li（2001）)，扰动LQA(亨特和李（2005）)或LLA(邹和李（2008）). 然而，这并没有充分利用CDA。实际上，罚函数的近似值需要额外的迭代才能收敛，这是不必要的，因为存在精确的协调解。因此，市场管理部、商业险核保部和分销部使用罚款的确切形式，只对损失进行控制。

现在，我们总结了给定（λ，γ）的MMCD。假设以下条件成立：

的二阶偏导数ℓ（β） w.r.t.β_j个为标准化一致有界X（X），即存在一个实数M（M）>0，这样 $\nabla_{j个}^{2} ℓ (β) \leq M（M）$ 对于所有β和j个= 0, …,第页.
inf公司_t吨ρ″(|t吨|;λ, γ) > −M（M），带ρ〃(|t吨|;λ、 γ）是ρ的二阶导数(|t吨|;λ、 γ）相对湿度|t吨|.

然后，市场管理部、商业险核保部和分销部按以下方式进行，

给定的初始值为 $\hat{β}$ ⁰,MMCD计算相应的线性分量 $\hat{η}$ ⁰.
对于m= 0,1, …,MMCD通过以下方式更新截距 (16),和使用 (14) 或 (15) 要更新 ${β̂}_{j个}^{米}$ 到 ${β̂}_{j个 + 1}^{米}$ 惩罚变量。每次迭代后，它还计算相应的线性分量 ${η̂}_{j个 + 1}^{米}$ 使用 (17).然后MMCD循环通过所有坐标，使得 $\hat{β}$ ^米已更新为 $\hat{β}$ ^米+1.
MMCD检查收敛标准。如果算法收敛，则停止迭代，否则重复步骤2直到算法收敛.

3.1 MMCD的收敛性分析

定理1建立了在一定条件下，MMCD解收敛到目标函数的极小值。

定理1 考虑目标函数 (2),其中给定数据(y、 X（X）)位于一个紧集合中，并且X的任何两列都不相同。对于给定的(λ, γ),假设惩罚ρ(|t吨|;λ, γ) ≡ ρ(t吨)满足ρ(t吨) = ρ(−t吨)和ρ′(|t吨|)是一个非负的一致有界函数ρ′(|t吨|)是的一阶导数（假设存在）ρ(|t吨|)水反应堆。|t吨|.还假设MMCD算法中规定的两个条件成立.

然后由MMCD算法生成的序列{β^米}收敛到函数Q的最小值(β;λ, γ).

上的条件(y、 X（X）)是温和的。标准化X（X）只要没有列为零，就可以执行。定理1的证明见附录.它扩展了Mazumder、Friedman和Hastie（2011）其SparseNet算法具有最小二乘损失。定理1涵盖了更一般的损失函数。

4惩罚Logistic回归的MMCD

如引言中所述，MMCD算法特别适用于逻辑回归，这是生物统计应用中使用最广泛的模型之一。对于逻辑回归模型，响应年是0或1的向量，其中1表示感兴趣的事件。损失函数的一阶和二阶导数为_j个ℓ( $\hat{β}$ ) = −(x个^j个)^T型(年− $\hat{π}$ )/n个和 $\nabla_{j个}^{2} ℓ (β̂) = {n个}^{- 1} \sum {w个}_{我} {x个}_{我 j个}^{2}$ ，使用w个_我= $\hat{π}$ _我(1 − $\hat{π}$ _我)和 $\hat{π}$ _我是我根据目前的估计，第次观察 $\hat{β}$ ，即。 ${ππ}_{我} = 1 / (1 + 经验 (- {x个}_{我}^{T型} β̂))$ 由于对于任何0≤π≤1，π（1−π）≤1/4(Böhning和Lindsay（1988）)，然后是的上限 $\nabla_{j个}^{2} ℓ (β̂)$ 是M（M）=1/4（标准化）x个^j个.相应的τ_j个= 4⁻¹ $\hat{β}$ _j个+n个⁻¹(x个^j个)^T型(年− $\hat{π}$ )的j个= 0, …,第页根据条件（b），我们要求SCAD惩罚γ>5，MCP惩罚γ>4。

4.1溶液表面的计算

应用SCAD和MCP惩罚的一种常见做法是计算固定值γ的λ解路径。例如，对于具有标准化变量的线性回归模型，建议对SCAD使用γ≈3.7(范和李（2001）)和γ≈2.7(张（2010）)用于MCP。然而，在包括逻辑回归的GLM模型中，这些值可能不合适。因此，我们使用数据驱动程序来选择γ和λ。这需要在（λ，γ）的二维网格上计算溶液表面。我们重新设定了κ=1/γ的参数，以便于描述计算溶液表面的方法。

定义[0，κ中的网格值_最大值) × [λ_最小值, λ_最大]0=κ₁≤ κ₂≤ ⋯ ≤ κ_K（K）< κ_最大值和λ_最大值= λ₁≥ λ₂≥ … ≥ λ_V（V）= λ_最小值.网格点的数量K（K）和V（V）是预先指定的。在我们的实现中，κ-网格点在标准尺度上是一致的，而λ的网格点在对数尺度上是统一的。κ_最大值是κ的最大值，因此MMCD算法的条件（b）有效。我们有κ_最大值=SCAD和κ的1/5_最大值=MCP的1/4。注意，当κ=0时，SCAD和MCP惩罚都成为拉索惩罚。λ_最大值是λ的最小值，使得 $\hat{β}$ _j个= 0,j个= 1, …,第页对于logistic回归模型，λ_最大值=n个⁻¹最大值_j个|(x个^j个)^T型(年− $\hat{π}$ )|对于每个κ_k个具有 $\hat{π}$ =ȳJ和J型元素都等于1的向量。我们设置λ_最小值= ελ_最大值，ε=0.0001，如果n个>第页ε=0.01。然后在矩形[0，κ上计算溶液表面_最大值) × [λ_最小值, λ_最大值]. 表示给定（κ）的MMCD解决方案_k个, λ_v（v）)由 $\hat{β}$ _{κ_k个,λ_v（v）}.

我们遵循的方法是Mazumder、Friedman和Hastie（2011）通过在λ值网格上的Lasso解处初始化MMCD算法来计算解曲面。拉索解对应于κ=0。然后，对于λ值网格中的每个点，我们从κ=0开始计算κ值网格上的解，将前一点的解用作当前点的初始值。该方法的细节如下。

该方法首先沿λ计算拉索解。计算时 $\hat{β}$ _{κ₀, λ_v（v）+1}，它使用 $\hat{β}$ _{κ₀, λ_v（v）}作为MMCD算法中的初始值。
对于给定λ_v（v）该方法沿κ计算解。这就是 $\hat{β}$ _{κ_k个, λ_v（v）}用作计算以下解的初始值 $\hat{β}$ _{κ_k个+1, λ_v（v）}.
然后，该方法循环通过v（v）=1…，V（V）用于步骤(2)以完成溶液表面。

定义一个变量为因果变量，如果其系数β≠0；否则，将其调用为null变量。图（1）使用MCP惩罚给出了β=2沿κ的因果变量的解路径。观察到，当κ与某些值交叉时，估计值可能会发生重大变化。这证明我们将κ视为一个调节参数是合理的，因为预先指定的κ可能不会给出最佳结果。这就是为什么我们更喜欢数据驱动程序来选择κ和λ的原因。

在单独的窗口中打开

图1

沿着κ的溶液路径图。它显示了β=2的因果变量的路径。观察到，当κ超过某些阈值时，估计值可能会发生实质性变化。

4.2模拟研究设计

让Z是协变量的设计矩阵，也就是说，它是X（X）删除了第一列。让A类₀≡ {1 ≤j个≤第页: β_j个≠0}是具有维数的因果变量集第页₀.我们修复第页₀= 10, β₀=0.0和系数A类₀为（±0.6，±1.2，2.4，∓0.6，\8723；1.2，±2.4）^T型信噪比（SNR）定义为 $S公司 N个 R（右） = \sqrt{β^{T型} {X（X）}^{T型} X（X） β / n个}$ ，大约在（3,4）的范围内。协变量是由均值和方差为零的多元正态分布∑生成的，∑是一个维数为的正定矩阵第页×第页响应变量的值由贝努利分布生成年_我~伯努利（1，第页_我)、和第页_我=经验（β^T型 x个_我)/（1+exp（β^T型 x个_我))的我= 1, …,n个.

我们考虑∑的五种相关结构。

（a）
独立结构（IN）第页惩罚变量。这里∑=我_第页，使用我_第页是维数的单位矩阵第页×第页.
（a）
独立结构（SP）。因果变量和零变量是独立的。设∑₀和∑₁分别为因果变量和零变量的协方差矩阵，则∑=块对角线（∑₀, Σ₁).在每一组变量中，我们假设一个复合对称结构，即ρ(x个_ij公司，x个_伊克)=ρj个≠k个.
（c）
部分相关结构（PC），即部分因果变量与部分零变量相关。具体来说，∑=块对角线（∑_一, Σ_b条, Σ_c（c）)，带∑_一是前5个因果变量的协方差矩阵；Σ_b条为剩余5个因果变量和5个零变量的协方差矩阵；Σ_c（c）剩余零变量的协方差矩阵。我们还假设∑内有一个复合对称结构_一, Σ_b条, Σ_c（c）.
（d）
一阶自回归（AR）结构，即ρ(x个_ij公司，x个_伊克) = ρ^{(|j个−k个|)}，用于j个≠k个.
（e）
复合对称（CS）结构第页变量。

4.3 LLA算法的数值实现

LLA的基本思想是近似凹罚ρ（|β_j个|; γ、 λ）依据 $ρ̇ (| {β̂}_{j个}^{米} |; γ, λ) | β_{j个} |$ 根据目前的估计 $\hat{β}$ ^米对于逻辑回归，我们也使用二次近似(7)损失依据 $\hat{β}$ ^米.要计算 $\hat{β}$ ^米+1使用LLA，我们最小化了Lasso类型标准

ℓ^{秒} (β | {β̂}^{米}) + \sum_{j个 = 1}^{第页} ρ̇ ({β̂}_{j个}^{米}; γ, λ) | β_{j个} | .

(18)

为了比较MMCD和LLA，我们用两种方法实现了LLA算法。第一个实现严格遵循中的描述邹和李（2008）这使用基于当前估算的工作数据，并将设计矩阵分为两部分， $U型 = {j个 : ρ̇ (| {β̂}_{j个}^{米} |; γ, λ) = 0}$ 和 $V（V） = {j个 : ρ̇ (| {β̂}_{j个}^{米} |; γ, λ) \neq 0}$ 用于当前估计 $\hat{β}$ ^米，带ρ̇(t吨)是ρ（·）的导数。计算 $\hat{β}$ ^米+1涉及 ${({X（X）}_{U型}^{* T型} {X（X）}_{U型}^{*})}^{- 1}$ ，使用 ${X（X）}_{U型}^{*} = ({X（X）}_{j个} : j个 \in U型)$ 是中变量的设计矩阵U型。因此，如果n个<第页_U型具有第页_U型是中的变量数U型因此，此方法通常仅适用于具有n个>第页.

在第二个实现中，我们使用CDA来最小化(18)。此实现可以使用第页≫n个。我们在下面将此实现称为LLA-CD算法。

因为这两种实现都需要初始估计 $\hat{β}$ 为了近似惩罚，我们使用拉索解沿κ启动LLA和LLA-CD算法的计算。LLA、自适应重缩放、LLA-CD和MMCD算法用类似的编程结构在FORTRAN中编程，以便进行公平比较。我们观察到，如果κ_最大值很大。因此，我们设置κ_最大值在我们的计算中，自适应缩放算法的值为0.25。我们使用收敛准则‖ $\hat{β}$ ^米+1− $\hat{β}$ ^米‖₂/(‖ $\hat{β}$ ^米‖₂+ δ) < ε. 我们选择δ=0.01和ε=0.001，如果n个>第页且ε=0.01，如果n个<第页.我们设置相关系数ρ=0.5，网格数K（K）= 10,V（V）= 100.

4.4计算效率比较

我们首先简要描述了可用的算法，即LQA、扰动LQA、LLA和MIST。这四个算法在使用代理函数优化惩罚ρ的意义上具有相同的精神(|t吨|; λ, γ). LQA使用以下近似值，

ρ (| t吨 |; λ, γ) \approx ρ (| {t吨}_{0} |; λ, γ) + \frac{ρ' (| {t吨}_{0} |; λ, γ)}{2 | {t吨}_{0} |} ({t吨}^{2} - {t吨}_{0}^{2}), t吨 \approx {t吨}_{0} .

(19)

然后采用Newton-Raphson型迭代，利用代理惩罚函数最小化惩罚准则。什么时候？t吨₀接近于零，算法不稳定。范和李（2001）建议如果 $\hat{β}$ _j个就足够小了| $\hat{β}$ _j个|<ε（预设值），一组 $\hat{β}$ _j个=0并删除j个迭代中的第个变量。因此，LQA算法的一个缺点是，一旦在任何迭代中删除了变量，就必须将其从最终模型中排除。亨特和李（2005）提出了LQA的扰动版本，以优化LQA。

ρ (| t吨 |; λ, γ) \approx ρ (| {t吨}_{0} |; λ, γ) + \frac{ρ' (| {t吨}_{0} |; λ, γ)}{2 | {t吨}_{0} + τ_{0} |} ({t吨}^{2} - {t吨}_{0}^{2}), t吨 \approx {t吨}_{0} .

(20)

实际上，如何确定τ的大小₀由于τ的大小₀可能会影响收敛速度和解的稀疏性。LLA算法（邹、李（2008）)将罚款近似为

ρ(|t吨|; λ, γ) ≈ ρ(|t吨₀|； λ, γ) + ρ′(|t吨₀|; λ, γ)(|t吨|−|t吨₀|), t吨 ≈ t吨₀.

(21)

水雾(Schifano、Strawderman和Wells（2010年）)该算法将LLA扩展到其他几个惩罚。

由于自适应重缩放方法只能应用于MCP惩罚，因此我们重点比较了MCP惩罚logistic回归模型的LLA、自适应重缩放、LLA-CD和MMCD的计算效率。计算是在Inter Xeon上完成的(W3540@2.93GHZ)使用Ubuntu 10.04操作系统（内核版本2.6）的机器。我们考虑两种设置n个>第页和n个<第页.

图2显示了基于100次复制的平均运行时间（以秒为单位）第页=100200和500，样本量固定n个= 1, 000. 观察到，当n个=1000和第页= 500. 这表明第页/n个对自适应缩放算法的效率影响较大。LLA-CD的速度也受到第页/n个比例，尽管程度较低。MMCD和LLA对于第页/n个比率。总的来说，MMCD是最快的。值得注意的是，在设置n个= 1,000,第页在一些复制中，当收敛标准ε=0.001时，自适应重缩放和LLA-CD在1000次迭代内没有收敛。

在单独的窗口中打开

图2

固定样本量的LLA、自适应缩放、LLA-CD和MMCD算法的计算效率(n个= 1, 000). 实心、虚线、虚线和带黑圈的虚线分别表示MMCD、LLA-CD、自适应缩放和LLA的平均消逝时间。这里，IN、SP、PC、AR和CS是指第4.2小节中描述的协变量的五种相关结构。

对于高维数据第页≫n个，我们重点比较了自适应缩放、LLA-CD和MMCD。我们使用了两种样本量n个=100和n个= 300. 变量的数量，第页，设置为500、1000、2000、5000和10000。图3显示基于100次复制的平均运行时间（秒）n个= 300. 它表明，作为第页增加，MMCD的优势变得更加明显。对于固定第页当预测因子相互关联时，市场管理部、商业险核保部和分销部的效率更高n个此外，MMCD的计算时间标准误差最小，其次是LLA-CD和自适应缩放。这表明MMCD是高维设置中考虑的三种算法中最稳定的一种。

在单独的窗口中打开

图3

自适应缩放、LLA-CD和MMCD算法的计算效率第页≫n个模型。样本大小为n个= 300. 实线、虚线和虚线分别表示MMCD、LLA-CD和自适应缩放的平均消逝时间。这里，IN、SP、PC、AR和CS是指第4.2小节中描述的协变量的五种相关结构。

4.5选择性能比较

我们进一步比较了MCP惩罚logistic模型的LLA、自适应重缩放、LLA-CD和MMCD的选择性能。由于本文中我们没有讨论调整参数选择的问题，因此将基于预测性能最佳的模型而不是特定调整参数选择方法选择的模型来比较算法。具体操作如下。我们首先计算[0，κ上的溶液表面_最大值) × [λ_最小值, λ_最大]基于训练数据集的每种算法。给出解决方案 $\hat{β}$ _{κ_k个,λ_v（v）}，我们计算ROC曲线下的预测面积（PAUC）AUC公司_{(κ_k个,λ_v（v）)}对于每个 $\hat{β}$ _{κ_k个,λ_v（v）}基于带有的验证集n个* = 3, 000. 最大预测值对应的模型AUC公司_{（κ_k个,λ_v（v）)}选择作为最终模型进行比较。我们比较了定义为所选变量总数的模型大小（MS）的结果；假发现率（FDR），定义为假阳性变量在所选总变量中的比例；验证数据集ROC曲线（PAUC）下的最大预测面积。以下报告的结果基于1000个重复。

表1介绍了中四种算法的比较n个>第页设置n个=1000和第页= 100. 结果表明，四种方法选择的模型具有相似的PAUC。就模型大小和FDR而言，MMCD和LLA-CD具有相似的性能，都具有比自适应重新缩放和LLA更小的模型大小和FDR。表2介绍了自适应重缩放、LLA-CD和MMCD在高维设置中与n个=100和第页= 2, 000. 与低维情况类似，三种方法的PAUC几乎相同。就模型大小和FDR而言，MMCD和LLA-CD具有相似的结果。

表1

四种算法的选择性能比较n个=1000和第页= 100. PAUC是指验证数据集ROC曲线（PAUC）下的最大预测区域。MS是模型尺寸。FDR是错误发现率。SE是基于1000个复制的标准错误。这里，IN、SP、PC、AR和CS指的是第4.2小节中描述的协变量的五种相关结构。

Σ （信噪比）	算法	PAUC公司（SE*10⁵)	微软（SE*10¹)	财务总监（SE*10^三)
英寸	拉	0.947(4.58)	10.96(0.55)	0.07(3.32)
(4.32)	广告资源	0.948(4.35)	16.28(1.21)	0.36(4.23)
	LLA-CD（LLA-CD）	0.948(3.45)	10.79(0.57)	0.06(3.09)
	市场管理部、商业险核保部和分销部	0.948(3.37)	10.90(0.56)	0.07(3.31)
服务提供商	拉	0.915(7.74)	11.39（0.77）	0.10(3.96)
(3.05)	广告资源	0.916(6.93)	14.14(0.87)	0.27(3.90)
	LLA-CD（LLA-CD）	0.917(7.24)	11.35(0.84)	0.10(4.10)
	市场管理部、商业险核保部和分销部	0.917(6.67)	11.27(0.64)	0.10(3.57)
个人计算机	拉	0.945(5.95)	14.25(1.50)	0.24(5.72)
(3.89)	广告资源	0.947(5.25)	15.55(1.09)	0.33(3.97)
	拉丁美洲和加勒比地区	0.947(5.61)	11.61(1.07)	0.11(4.46)
	市场管理部、商业险核保部和分销部	0.947(5.07)	11.41(0.79)	0.10(3.93)
AR公司	拉	0.921(8.83)	13.83(1.28)	0.24(5.34)
（3.20）	广告资源	0.924(6.73)	18.76(1.34)	0.44(3.94)
	LLA-CD（LLA-CD）	0.924(7.88)	11.29(0.76)	0.10(3.49)
	市场管理部、商业险核保部和分销部	0.925(5.98)	12.11(0.82)	0.15(4.45)
反恐精英	拉	0.919（8.13）	12.42(1.07)	0.16(4.70)
(3.06)	广告资源	0.921(7.02)	14.15(0.90)	0.27(3.98)
	LLA-CD（LLA-CD）	0.922(6.61)	10.64(0.54)	0.05（2.80）
	市场管理部、商业险核保部和分销部	0.922(6.60)	10.94(0.58)	0.07(3.32)

在单独的窗口中打开

表2

自适应缩放、LLA-CD和MMCD的选择性能比较n个=100和第页= 2, 000. PAUC是指验证数据集ROC曲线（PAUC）下的最大预测区域。MS是模型尺寸。FDR是错误发现率。SE是基于1000次复制的标准错误。这里，IN、SP、PC、AR和CS是指第4.2小节中描述的协变量的五种相关结构。

Σ （信噪比）	算法	PAUC公司（SE*10⁵)	微软（数量*10¹)	财务总监（SE*10^三)
英寸	广告资源	0.828(1.30)	12.25(3.01)	0.60(6.95)
（4.33）	拉丁美洲和加勒比地区	0.842(1.28)	5.56(2.02)	0.25(8.42)
	市场管理部、商业险核保部和分销部	0.844(1.27)	6.41(2.09)	0.28(9.06)
服务提供商	广告资源	0.778(1.96)	12.06(3.77)	0.62(6.05)
（3.05）	LLA-CD（LLA-CD）	0.795(1.74)	5.25(2.15)	0.26(8.16)
	市场管理部、商业险核保部和分销部	0.797(1.76)	5.75(2.25)	0.28(8.34)
个人计算机	广告资源	0.872(0.64)	7.12(1.37)	0.42（6.43）
(3.87)	LLA-CD（LLA-CD）	0.877(0.54)	5.19(1.29)	0.24(7.48)
	市场管理部、商业险核保部和分销部	0.877(0.54)	5.37(1.27)	0.26(7.46)
AR公司	广告资源	0.812(1.21)	6.21(1.69)	0.49（8.53）
(3.19)	LLA-CD（LLA-CD）	0.830(1.10)	3.02(0.71)	0.17(7.73)
	市场管理部、商业险核保部和分销部	0.831(1.07)	3.21(0.94)	0.18(8.10)
反恐精英	广告资源	0.770(1.79)	11.89(3.58)	0.64(6.32)
(3.04)	LLA-CD（LLA-CD）	0.776(1.80)	6.99(3.09)	0.37(9.27)
	市场管理部、商业险核保部和分销部	0.781(1.80)	7.39（2.99）	0.39(9.49)

在单独的窗口中打开

4.6癌症基因表达数据集的应用

本研究的目的是发现与乳腺癌预后相关的生物标志物(范特维尔等(2002);范德维杰等(2002)). 使用微阵列对约25000个基因进行扫描n个=295名患者。五年内转移被建模为结果。惩罚模型中使用了1000个与结果具有最高Spearman相关性的基因子集来稳定计算。出于与仿真研究中相同的原因，我们不采用任何调整参数选择程序来选择模型进行比较。相反，我们随机地对整个数据集进行分区n个=295分为训练数据集（约1/3的观察值）和验证数据集（大约2/3的观测值）。模型拟合仅基于训练数据集；选择与验证数据集的最大预测AUC对应的解决方案作为最终模型进行比较。我们重复这个随机划分过程900次。

表3显示了使用MMCD的SCAD惩罚结果，以及使用自适应缩放、LLA-CD和MMCD的MCP惩罚结果。来自不同方法的PAUC彼此接近。LLA-CD算法的MCP惩罚的模型大小恰好是最大的。

表3

SCAD和MCP在微阵列数据集中的应用。平均误差和标准误差是基于900个分割过程计算的。预测AUC计算为随机分割过程创建的验证数据集的最大预测AUC。在每个拆分过程中，大约n个=100个样本分配给训练数据集n个=200个样本进入验证数据集。

溶液表面	贫困人口（人口*10^三)	理学硕士（SE）
SCAD（市场管理部、商业险核保部）	0.7567 (0.99)	35.50(0.47)
MCP（适配器）	0.7565 (1.15)	39.06 (0.68)
MCP（LLA-CD）	0.7537 (0.99)	43.07 (0.63)
MCP（毫米CD）	0.7570(0.99)	35.66 (0.49)

在单独的窗口中打开

最后，我们使用ROC曲线下交叉验证区域（CV-AUC）作为调谐参数选择方法，展示了乳腺癌研究的结果。该方法结合了交叉验证和ROC方法。有关在惩罚logistic回归中使用CV-AUC调整参数选择的详细信息，请参阅江、黄、张（2011）我们使用5倍交叉验证来计算CV-AUC。对于该数据集，使用MMCD的SCAD、使用自适应缩放的MCP和LLA-CD选择相同的模型，其中包含67个变量，CV-AUC=0.7808。使用MMCD的MCP选择16个CV-AUC=0.8024的变量。

5讨论

在本文中，我们提出了计算凹惩罚解的MMCD算法。我们的仿真研究和数据示例表明，该算法在凹惩罚logistic回归模型中是有效的第页≫n个不同于现有的算法，例如近似惩罚函数的LQA、LLA和MIST，MMCD使用精确的惩罚项为每个坐标寻找闭合形式的解。优化仅适用于损失函数。这种方法提高了CDA在高维环境中的效率。在合理的条件下证明了MMCD的收敛性。

LLA、自适应重缩放、LLA-CD和MMCD之间的比较表明，MMCD是更有效的方法，尤其是对于大型第页和相关协变量，尽管LLA-CD在某些情况下具有竞争力。LLA-CD实现了相邻启动思想以降低计算成本，即使用 $\hat{β}$ _{κ_k个,λ_v（v）}作为要计算的初始值 $\hat{β}$ _{κ_k个+1, λ_v（v）}。在CDA组件中，解决方案按顺序更新，即使用 ${β̂}_{j个}^{米 + 1}$ 计算 ${β̂}_{j个 + 1}^{米 + 1}$ ，而不是矢量形式，它使用 $\hat{β}$ ^米计算 $\hat{β}$ ^米+1这与LLA-LARS的实施不同Breheny和Huang（2011）相邻的启动和顺序更新方案可能是LLA的两种实现表现如此不同的主要原因。

MMCD在logistic回归中的应用是因为可以为logistic似然构造一个简单有效的优化函数。其他GLM模型也可以找到优化函数，例如基线类别逻辑模型。在这种模型中，市场管理部、商业险核保部和分销部可以以类似的方式实施。然而，在GLM家族中的其他一些重要模型（如对数线性模型）中，似乎不存在简单的优化函数。一种可能的方法是根据每次迭代的解设计一系列优化函数。这是一个需要进一步调查的有趣问题。

致谢

作者感谢审稿人和副主编的有益评论，这些评论使论文有了很大的改进。Huang的研究得到了美国国立卫生研究院拨款R01CA120988、R01CA142774和美国国家科学基金会拨款DMS 1208225的部分支持。

附录

在附录中，我们证明了定理1。证明遵循以下基本思想Mazumder、Friedman和Hastie（2011）然而，也有一些重要的区别。特别是，我们需要注意引理1和定理1中的截距、损失函数的二次逼近和定理1的协调优化。

引理1 假设数据(y、 X（X）)位于紧集上，且以下条件成立：

损失函数ℓ(β)是（完全）可微的w.r.t。β对于任何β ∈ ℝ^第页+1.
惩罚函数ρ(t吨)关于0是对称的，并且在t上是可微的≥ 0; ρ′(|t吨|)非负、连续且一致有界，其中ρ′(|t吨|)是的导数ρ(|t吨|)水反应堆。|t吨|.
序列{β^k个}是有界的.
对于每个收敛子序列{β^n个_k个}⊂｛β^n个},连续差值收敛到零：β^n个_k个− β^{n个_k个− 1}→ 0

那么如果β^∞ 是序列的任何极限点{β^k个},然后β^∞ 是函数Q的最小值(β);即

\underset{α ↓ 0 +}{lim信息} {\frac{问 (β^{\infty} + α δ) - 问 (β^{\infty})}{α}} \geq 0,

(22)

对于任何δ = (δ₀, …, δ_第页) ∈ ℝ^第页+1.

证明对于任何β=（β₀，…，β_第页)^T型和δ_j个= (0, …, δ_j个, …, 0) ∈ ℝ^第页+1，我们有

\underset{α ↓ 0 +}{lim信息} {\frac{问 (β + α δ_{j个}) - 问 (β)}{α}} = \nabla_{j个} ℓ (β) δ_{j个} + \underset{α ↓ 0 +}{lim信息} {\frac{ρ (| β_{j个} + α δ_{j个} |) - ρ (| β_{j个} |)}{α}} = \nabla_{j个} ℓ (β) δ_{j个} + \partial ρ (β_{j个}; δ_{j个}),

(23)

对于j个∈ {1, …,第页}，使用

\partial ρ (β_{j个}; δ_{j个}) = {\begin{matrix} ρ' (| β_{j个} |) sgn公司 (β_{j个}) δ_{j个}, & | β_{j个} | > 0; \\ ρ' (0) | δ_{j个} |, & | β_{j个} | = 0, \end{matrix}

(24)

假设 $β^{{n个}_{k个}} \to β^{\infty} = (β_{0}^{\infty}, \dots, β_{第页}^{\infty})$ ，根据假设4k个→ ∞

β_{j个}^{{n个}_{k个} - 1} = (β_{0}^{{n个}_{k个}}, \dots, β_{j个 - 1}^{{n个}_{k个}}, β_{j个}^{{n个}_{k个}}, β_{j个 + 1}^{{n个}_{k个} - 1}, \dots, β_{第页}^{{n个}_{k个} - 1}) \to (β_{0}^{\infty}, \dots, β_{j个 - 1}^{\infty}, β_{j个}^{\infty}, β_{j个 + 1}^{\infty}, \dots, β_{第页}^{\infty})

(25)

由(24)和(25)，我们有以下结果j个∈ {1, …,第页}.

\partial ρ (β_{j个}^{{n个}_{k个}}; δ_{j个}) \to \partial ρ (β_{j个}^{\infty}; δ_{j个}) 如果 β_{j个}^{\infty} \neq 0; \partial ρ (β_{j个}^{\infty}; δ_{j个}) \geq \underset{k个}{lim信息} \partial ρ (β_{j个}^{{n个}_{k个}}; δ_{j个}) 如果 β_{j个}^{\infty} = 0 .

(26)

按协调最小值j个第th坐标j个∈ {1, …,第页}，我们有

\nabla_{j个} ℓ (β_{j个}^{{n个}_{k个} - 1}) δ_{j个} + \partial ρ (β_{j个}^{{n个}_{k个}}; δ_{j个}) \geq 0, 为所有人 k个 .

(27)

因此(26,27)意味着所有人j个∈ {1, …,第页},

\nabla_{j个} ℓ (β^{\infty}) δ_{j个} + \partial ρ (β_{j个}^{\infty}; δ_{j个}) \geq \underset{k个}{lim信息} {\nabla_{j个} ℓ (β_{j个}^{{n个}_{k个} - 1}) δ_{j个} + \partial ρ (β_{j个}^{{n个}_{k个}}; δ_{j个})} \geq 0 .

(28)

由(23,28)，用于j个∈ {1, …,第页}，我们有

\underset{α ↓ 0 +}{lim信息} {\frac{问 (β^{\infty} + α δ_{j个}) - 问 (β^{\infty})}{α}} \geq 0 .

(29)

根据上述论点，很容易看出j个= 0

Ş₀ℓ(β^∞)δ₀ ≥ 0.

(30)

因此，对于δ=（δ₀, …, δ_第页) ∈ ℝ^第页+1，根据的可微性ℓ（β），我们有

\underset{α ↓ 0 +}{lim信息} {\frac{问 (β^{\infty} + α δ) - 问 (β^{\infty})}{α}} = \nabla_{0} ℓ (β^{\infty}) δ_{0} + \sum_{j个 = 1}^{第页} [\nabla_{j个} ℓ (β^{\infty}) δ_{j个} + \underset{α ↓ 0 +}{lim信息} {\frac{ρ (| β_{j个}^{\infty} + α δ_{j个} |) - ρ (| β_{j个}^{\infty} |)}{α}}] = \nabla_{0} ℓ (β^{\infty}) δ_{1} + \sum_{j个 = 1}^{第页} \underset{α ↓ 0 +}{lim信息} {\frac{问 (β^{\infty} + α δ_{j个}) - 问 (β^{\infty})}{α}} \geq 0,

（31）

通过(29,30)这就完成了证明。

定理1的证明

证明要简化符号，请写 $χ_{β_{0}, \dots, β_{j个 - 1}, β_{j个 + 1}, \dots, β_{第页}}^{j个} \equiv χ (u个)$ 对于问（β）作为的函数j个th坐标与β_我,我≠j个正在修复。我们首先处理j个∈ {1, …,第页}坐标，然后是以下参数中的截距（第0个坐标）。

对于j个∈ {1, …,第页}坐标，观察

χ(u个+ δ)−χ(u个) = ℓ(β₀, …, β_j个−1, u个+ δ, β_j个+1, …, β_第页)−ℓ(β₀, …, β_j个−1, u个, β_j个+1, …, β_第页) + ρ(|u个+ δ|)−ρ(|u个|)

(32)

= \nabla_{j个} ℓ (β_{0}, \dots, β_{j个 - 1}, u个, β_{j个 + 1}, \dots, β_{第页}) δ + \frac{1}{2} \nabla_{j个}^{2} ℓ (β_{0}, \dots, β_{j个 - 1}, u个, β_{j个 + 1}, \dots, β_{第页}) δ^{2} + o（o） (δ^{2}) + ρ' (| u个 |) (| u个 + δ | - | u个 |) + \frac{1}{2} ρ ″ (| {u个}^{*} |) {(| u个 + δ | - | u个 |)}^{2},

(33)

带有|u个*|介于|u个+δ|和|u个|。符号+_j个ℓ(β₀，…，β_j个−1,u个, β_j个+1, …, β_第页)和 $\nabla_{j个}^{2} ℓ (β_{0}, \dots, β_{j个 - 1}, u个, β_{j个 + 1}, \dots, β_{第页})$ 表示函数的一阶导数和二阶导数ℓ w.r.t.的j个th坐标（假设根据条件存在（1）).

我们重新编写的RHS(33)如下：

R（右） H（H） S公司 (属于 33) = \nabla_{j个} ℓ (β_{0}, \dots, β_{j个 - 1}, u个, β_{j个 + 1}, \dots, β_{第页}) δ + (\nabla_{j个}^{2} ℓ (β_{0}, \dots, β_{j个 - 1}, u个, β_{j个 + 1}, \dots, β_{第页}) - M（M）) δ^{2} + ρ' (| u个 |) sgn公司 (u个) δ + ρ' (| u个 |) (| u个 + δ | - | u个 |) - ρ' (| u个 |) sgn公司 (u个) δ + \frac{1}{2} ρ ″ (| {u个}^{*} |) {(| u个 + δ | - | u个 |)}^{2} + (M（M） - \frac{1}{2} \nabla_{j个}^{2} ℓ (β_{0}, \dots, β_{j个 - 1}, u个, β_{j个 + 1}, \dots, β_{第页})) δ^{2} + o（o） (δ^{2}) .

(34)

另一方面j个th坐标(j个∈{1…，第页})是最小化以下功能，

问_{j个} (u个 | β) = ℓ (β) + \nabla_{j个} ℓ (β) (u个 - β_{j个}) + \frac{1}{2} \nabla_{j个}^{2} ℓ (β) {(u个 - β_{j个})}^{2} + ρ (| u个 |),

(35)

通过控制，我们 $\nabla_{j个}^{2} ℓ (β)$ 通过一个常数M（M）用于标准化变量。所以实际最小化的函数是

\tilde{问}

_j个(u个|β) = ℓ(β) + ∇_j个ℓ(β)(u个 − β_j个)+½M（M）(u个−β_j个)²+ ρ(|u个|).

（36）

自u个就是最小化（36），我们已经为j个第个(j个∈ {1, …,第页})坐标，

∇_j个ℓ(β) + M（M）(u个 − β_j个) + ρ′(|u个|)sgn公司(u个)＝0时

(37)

因为χ(u个)最小化为u个₀，由(37)，我们有

0 = \nabla_{j个} ℓ (β_{0}, \dots, β_{j个 - 1}, {u个}_{0} + δ, β_{j个 + 1}, \dots, β_{第页}) + M（M） ({u个}_{0} - {u个}_{0} - δ) + ρ' (| {u个}_{0} |) sgn公司 ({u个}_{0}) = \nabla_{j个} ℓ (β_{0}, \dots, β_{j个 - 1}, {u个}_{0}, β_{j个 + 1}, \dots, β_{第页}) + \nabla_{j个}^{2} ℓ (β_{0}, \dots, β_{j个 - 1}, {u个}_{0}, β_{j个 + 1}, \dots, β_{第页}) δ + o（o） (δ) - M（M） δ + ρ' (| {u个}_{0} |) sgn公司 ({u个}_{0}),

(38)

如果u个₀=0，则上述对sgn的某些值成立(u个₀) ∈ [−1, 1].

观察ρ′(|x个|)≥0，则

ρ′(|u个|)(|u个+ δ|−|u个|)负极ρ′(|u个|)sgn公司(u个)δ= ρ′(|u个|)[(|u个+ δ|−|u个|)−sgn(u个)δ] ≥ 0

(39)

因此使用(38,39)在里面(34)在u个₀，我们有，因为j个∈{1…，第页},

χ ({u个}_{0} + δ) - χ ({u个}_{0}) \geq \frac{1}{2} ρ ″ (| {u个}^{*} |) {(| u个 + δ | - | u个 |)}^{2} + δ^{2} (M（M） - \frac{1}{2} \nabla_{j个}^{2} ℓ (β_{0}, \dots, β_{j个 - 1}, {u个}_{0}, β_{j个 + 1}, \dots, β_{第页})) + o（o） (δ^{2}) \geq \frac{1}{2} M（M） δ^{2} + \frac{1}{2} ρ ″ (| {u个}^{*} |) {(| u个 + δ | - | u个 |)}^{2} + o（o） (δ^{2}) .

(40)

根据MMCD算法inf的条件（b）_t吨ρ″(|t吨|;λ, γ) > −M（M）和(|u个+ δ| − |u个|)²≤ δ²。因此存在θ₂= ½(M（M）+inf公司_x个ρ”(|x个|) + o（o）（1））>0，因此对于j个第个坐标，j个∈ {1, …,第页},

χ(u个₀+ δ)−χ(u个₀) ≥ θ₂δ².

(41)

现在考虑β₀，观察到

χ (u个 + δ) - χ (u个) = ℓ (u个 + δ, β_{1}, \dots, β_{第页}) - ℓ (u个, β_{1}, \dots, β_{第页}) = \nabla_{1} ℓ (u个, β_{1}, \dots, β_{第页}) δ + \frac{1}{2} \nabla_{1}^{2} ℓ (u个, β_{1}, \dots, β_{第页}) δ^{2} + o（o） (δ^{2}) = \nabla_{1} ℓ (u个, β_{1}, \dots, β_{第页}) δ + (\nabla_{1}^{2} (ℓ (u个, β_{1}, \dots, β_{第页}) - M（M）) δ^{2} + (M（M） - \frac{1}{2} Ş_{1}^{2} ℓ (u个, β_{1}, \dots, β_{第页})) δ^{2} + o（o） (δ^{2}),

(42)

通过类似的论据(38)，我们有

0 = \nabla_{1} ℓ ({u个}_{0} + δ, β_{1}, \dots, β_{第页}) + M（M） ({u个}_{0} + δ - {u个}_{0}) = \nabla_{1} ℓ ({u个}_{0}, β_{1}, \dots, β_{第页}) + \nabla_{1}^{2} ℓ ({u个}_{0}, β_{1}, \dots, β_{第页}) δ + o（o） (δ) - M（M） δ .

(43)

因此，通过(42,43)，对于β的第一个坐标

χ ({u个}_{0} + δ) - χ ({u个}_{0}) = (M（M） - \frac{1}{2} \nabla_{1}^{2} ℓ ({u个}_{0}, β_{1}, \dots, β_{第页})) δ^{2} + o（o） (δ^{2}) = \frac{1}{2} M（M） δ^{2} + \frac{1}{2} (M（M） - \nabla_{1}^{2} ℓ ({u个}_{0}, β_{1}, \dots, β_{第页})) δ^{2} + o（o） (δ^{2}) \geq \frac{1}{2} δ^{2} (M（M） + o（o） (1)) .

(44)

因此存在一个θ₁= ½(M（M）+ o（o）（1））>0，因此对于β的第一个坐标

χ(u个₀+ δ)−χ(u个₀) ≥ θ₁δ².

(45)

设θ=min（θ₁, θ₂)，使用(41,45)，我们有β的所有坐标，

χ(u个₀+ δ)−χ(u个₀) ≥ θδ²,

(46)

由(46)我们有

问 (β_{j个}^{米 - 1}) - 问 (β_{j个 + 1}^{米 - 1}) \geq θ {(β_{j个 + 1}^{米} - β_{j个 + 1}^{米 - 1})}^{2} = θ {‖ β_{j个}^{米 - 1} - β_{j个 + 1}^{米 - 1} ‖}_{2}^{2},

(47)

哪里 $β_{j个}^{米 - 1} = (β_{1}^{米}, \dots, β_{j个}^{米}, β_{j个 + 1}^{米 - 1}, \dots, β_{第页}^{米 - 1})$ . The(47)确定序列{β的有界性^米}对于每个米>从{β的起点开始¹} ∈ ℝ^第页+1.

应用(47)在所有坐标上，我们都有米

问 (β^{米}) - 问 (β^{米 + 1}) \geq θ {‖ β^{米 + 1} - β^{米} ‖}_{2}^{2} .

(48)

由于（递减）顺序问(β^米)汇聚，(48)显示序列{β^k个}有一个唯一的极限点。这就完成了{β的收敛性证明^k个}.

假设(3)和（4）在引理1中由(48)因此，{β的极限点^k个}最小值为问（β）通过引理1。这就完成了定理的证明。

脚注

补充材料

MMCD算法的R包：R包“cvplogistic”可在网址：www.r-project.org(R开发核心团队（2011）针对具有凹惩罚的logistic回归，实现了自适应重标度、LLA-CD和MMCD算法。

参与者信息

蒋定峰，探索统计、数据和统计科学，AbbVie Inc。moc.liamg@gnaijgnefgnid公司.

黄健，爱荷华大学统计与精算科学系和生物统计学系。

参考文献

Böhning D，Lindsay B.二次逼近算法的单调性。Ann.Inst.Stat.数学。1988;40(4):641–663. [谷歌学者]
Breheny P，Huang J.非凸惩罚回归的坐标下降算法，及其在生物特征选择中的应用。附录申请。斯达。2011;5(1):232–253. [PMC免费文章][公共医学][谷歌学者]
Donoho DL，Johnstone JM。小波收缩的理想空间自适应。生物特征。1994年；81(3):425–455. [谷歌学者]
Efron B，Hastie T，Johnstone I，Tibshirani R.最小角回归。Ann.统计。2004;32(2):407–451. [谷歌学者]
Fan J，Li R.通过非凹惩罚似然的变量选择及其预言性质。美国统计协会。2001;96(456):1348–13608. [谷歌学者]
Friedman J，Hastie T，Höfling H，Tibshirani R.路径坐标优化。附录申请。斯达。2007;1(2):302–332. [谷歌学者]
Friedman J，Hastie T，Tibshirani R.通过坐标下降法实现广义线性模型的正则化路径。J.统计软件。2010;33(1):1–22. [PMC免费文章][公共医学][谷歌学者]
Hunter DR，Lange K。MM算法教程。美国统计局。2004;58(1):30–37. [谷歌学者]
Hunter DR，Li R.使用MM算法选择变量。Ann.统计。2005;33(4):1617–1642. [PMC免费文章][公共医学][谷歌学者]
蒋D，黄J，张勇。高维数据MCP-Logistic回归的交叉验证AUC。统计方法。医学研究。2011认可的。[公共医学][谷歌学者]
兰格·K。优化。纽约：Springer；2004[谷歌学者]
Lange K，Hunter D，Yang I.使用替代目标函数优化转移（含讨论）J.计算。图表。斯达。2000;9(1):1–59. [谷歌学者]
Mazumder R、Friedman J、Hastie T。备用网络：协调下降与非凸惩罚。美国统计协会。2011;106(495):1125–1138. [PMC免费文章][公共医学][谷歌学者]
Ortega JM，Rheinbold WC。多元非线性方程的迭代解法。纽约州纽约市：学术出版社；1970[谷歌学者]
Osborne MR，Presnell B，Turlach BA。最小二乘问题中变量选择的新方法。IMA公司。J.数字。分析。2000;20（3）：389–403。 [谷歌学者]
Schifano ED，Strawderman RL，Wells MT。非光滑惩罚目标函数的优化最小算法。电子统计杂志。2010;4:1258–1299. [谷歌学者]
Tibshirani R.通过Lasso回归收缩和选择。J.R.统计社会服务。B。1996;58(1):267–288. [谷歌学者]
无微分极小化的块坐标下降法的Tseng P.收敛性。J.Optimiz。理论。应用程序。2001;109(3):475–494. [谷歌学者]
van t Veer LJ，Dai H，van de Vijver MJ等。基因表达谱预测乳腺癌的临床结局。自然。2002;415(31):530–536.[公共医学][谷歌学者]
van de Vijver MJ、He YD、van t Veer LJ等。基因表达特征作为乳腺癌生存预测因子。北英格兰。医学杂志。2002;347（25）：1999-2009年。[公共医学][谷歌学者]
Warge J.最小化某些凸函数。SIAM J.应用。数学。1963;11(3):588–593. [谷歌学者]
Wu TT，Lange K。拉索惩罚回归的坐标下降算法。附录申请。斯达。2008年；2(1):224–244. [PMC免费文章][公共医学][谷歌学者]
张川。极小极大凹罚下的几乎无偏变量选择。Ann.统计。2010;38(2):894–942. [谷歌学者]
Zou H，Li R.非洞穴惩罚似然模型中的一步稀疏估计。Ann.统计。2008年；36(4):1509–1533. [PMC免费文章][公共医学][谷歌学者]
R开发核心团队。R统计计算基金会。奥地利维也纳：ISBN；R：统计计算的语言和环境。3-900051-07-0,网址：http://www.R-project.org.[谷歌学者]