凸优化（4）——次梯度案例，加速梯度法，随机梯度下降法，近端梯度法引入

学弱猹

发布于 2021-08-09 17:27:12

1.4K0

发布于 2021-08-09 17:27:12

文章被收录于专栏：学弱猹的精品小屋学弱猹的精品小屋

上一节笔记：凸优化（3）——梯度与次梯度：方法，性质与比较

————————————————————————————————————————

大家好！

这一节我们开始把我们之前与梯度法和次梯度法有关的，但是还没有说完的部分说完。还有篇幅的话，就再谈一谈随机梯度下降方法。

那么我们开始吧。

Source

Hung-yi Lee, Machine Learning (videos on Bilibili)
CMU 10-725, Convex Optimization
Amir Beck (2017), FIrst-order Methods in Optimization
Michael W. Mahoney, et al. The Mathematics of Data
Boyd, Vandenberghe, Convex Optimization
L. Vandenberghe, Lecture notes for EE 236C, UCLA, Spring 2011-2012
D. Bersekas (2010), Incremental gradient, subgradient, and proximal methods for convex optimization: a survey
A. Nemirovski and A. Juditsky and G. Lan and A. Shapiro (2009), Robust stochastic optimization approach to stochastic programming
Ryan Tibshirani (2015), A general framework for fast stagewise algorithms
B. Polyak (1987), Introduction to optimization Chapter 5

梯度方法，次梯度方法的加速

看到这个主题，你可能第一眼会觉得，我们是要谈一谈方法的加速。但是实际上我们要先泼一盆冷水：我们不会提次梯度方法的加速。明明不会提，为什么要起这个标题？我们后面讲到这里大家就明白了。

我们先提一些梯度方法的延伸，也就是它们的一些可能的加速策略。

动量方法

我们要提到的梯度的加速方法都不可避免的和动量（momentum）这个概念相关，我们后面会解释这个概念。

波利亚重球法

对，就是那个上一节提到的波利亚（Polyak）。

波利亚重球法（Polyak Heavy-Ball Method）严格来说不是一个加速方法，但可以帮助我们很好的理解动量这样的概念。回顾一下梯度方法的迭代公式

x_{k + 1} = x_k - \alpha_k \nabla f(x_k)

那么对于重球法来说其实就是这么一个迭代公式

x_{k + 1} = x_k - \alpha_k \nabla f(x_k) + \beta_k (x_k - x_{k - 1})

你会发现，这个公式就多了一项，这一项是

x_k - x_{k - 1}

的倍数。所以实际上如果我们认为

\nabla f(x_k)

是当前时间点的物体运动的速度，那么

x_k - x_{k - 1}

事实上可以近似理解为物体在之前一个时间点的速度，也就是物理中动量的概念（毕竟物体的动量为

p = mv

嘛）。为什么要引入物理的概念呢？一个比较好的解释可以见下面这张图。

可以看出，如果只有梯度的话，球到达一个局部极小值点，就会停下来。但是对于这张图而言，其实球可以往下再走一点，“越过”这个山坡，然后到达更小的值。那么动量的作用就相当于认为球”不应该停留在那个局部极小值“，所以就完成了这个任务。

所以说重球法其实目的是希望函数迭代到更小的值，和加速没什么关系。但是我们到后面提到的方法，其实都和这个动量的概念有关。

加速梯度方法

好的，现在我们进入正题了。我们先说一句：共轭梯度法（Conjugate Gradient）本身也算是一种加速优化方法。列出它的原因在于共轭梯度法也只使用了一阶信息，但是因为没有直接使用梯度，所以它不属于加速梯度方法。这里我们贴出《数值优化》第3和4节的链接，大家可以点进去了解共轭梯度法。

数值优化（3）——线搜索中的步长选取方法，线性共轭梯度法

数值优化（4）——非线性共轭梯度法，信赖域法

Nesterov加速梯度法

加速梯度法最有名的就是Nesterov加速梯度法，但必须要承认它的原理是非常难以理解的，这也是我们不会在这里详细证明它的性质的原因（没啥必要……）。

对于Nesterov加速方法，它的迭代公式为

x_{k + 1} = x_k - \alpha_k \nabla f(x_k + \beta_k (x_k - x_{k - 1})) + \beta_k(x_k - x_{k - 1})

还有一种拆分而连贯的写法是

x_{k + 1} = y_k - \frac 1 L \nabla f(y_k)

y_{k + 1} = x_{k + 1} = \beta_{k + 1}(x_{k + 1} - x_k)

然后设初值

x_0 = y_0

，人工给定。

这一个公式的关键变化就在于梯度内添加了一个动量项

\beta_k (x_k - x_{k - 1})

，这个方法很强大的地方在于它可以达到

O(\frac1 {k^2})

的收敛速度，并且**对

f(x)

并没有强凸的要求**。这个就是下面这个定理揭示的内容。

Theorem 1: 设

f(x)

凸，

\nabla f(x)

Lipschitz连续，且Lipschitz常数为

，那么如果设

\lambda_0 = 0, \lambda_{k + 1}^2 - \lambda_{k + 1} = \lambda_k^2, \beta_k = \frac{\lambda_k - 1}{\lambda_{k + 1}}

，那么有

f(x_T) - f(x^*) \le \frac{2L\|x_0 - x^* \|^2_2}{(T + 1)^2}

，其中

是迭代步数，

f(x^*)

为最小值。

通过这个上界就可以得到收敛速率的结论。

一个自然的问题就是：如果有强凸，加速梯度法有没有作用？事实上也是一样，可以从次线性收敛速度推到线性收敛速度。当然具体的原因就不解释了，我也解释不明白这玩意。

为什么说这一个方法之后就不再继续说了呢？一方面当然是因为其他方法没有特别有名，另一方面，也是因为这是一阶方法能够达到的极限。换句话说，如果完全不考虑问题的特性，希望一个万金油的方法，那么最好的收敛速度也就是

O(\frac 1{k^2})

。下面这个性质就解释了这个现象。

Proposition 1: 对任意

k \le \frac {n -1}2

，任意出发点

x^{(0)}

，都存在一个函数

，使得任何一阶方法都有

f(x^{(k)}) - f(x^*) \ge \frac{3L\|x_0 - x^*\|_2^2}{32(k + 1)^2}

，这里

表示数据的维数

不过这个定理严格来说不是很严谨，有的人可能会问

k > \frac{n - 1}2

怎么办。事实上这就要交给其他的性质了，我们这里不详谈，感兴趣的可以参考Source中的The Mathematics of Data相关的章节。

当然了这里为了严谨，多说一句，一阶方法的严格定义就是每一步迭代，值都落在

x^{(0)} + \operatorname{span} \{\nabla f(x_0), \nabla f(x_1), \ldots, \nabla f(x_{k - 1})\}

这个线性空间内。

最后我们再简单多说几句与次梯度方法相关的加速情况。我们上一节说对于次梯度方法，它的收敛速度一般是

O(\frac 1{\sqrt {k}})

。但是这个收敛速度能改进吗？比方说我们添加对函数凸的要求什么的。但很遗憾，次梯度方法无法进行加速，也是有下面这个性质保证的。

Proposition 2: 对任意

k \le n - 1

和出发点

x^{(0)}

，存在一个

，使得任何一个非光滑（梯度不存在）的一阶方法，都有

f(x^{(k)}) - f(x^*) \ge \frac{RG}{2(1 + \sqrt{k + 1})}

如果你忘记了

R, G

的含义，记得查看上一节。

好的，关于两种不同的方法的延伸和一些与加速有关的核心性质，我们这里就介绍完了。

次梯度方法的案例分析

这一部分事实上就是上一节跳过的，与次梯度有关的一些部分。因为和梯度不太一样的事，次梯度在分析上其实难度是大很多的，实际的应用场景也广泛很多，所以值得多提一些相关的案例。

首先我们来看一下机器学习里面非常有名的LASSO。

Example 1: LASSO（Source：CMU）考虑问题

y \in \mathbb R^n, X \in \mathbb R^{n \times p}

，观察优化问题

\min_\beta \frac 12 \|y - X \beta \|_2^2 + \lambda \|\beta \|_1

，

\lambda \ge 0

的最优性条件。

注意到

\lambda \|\beta \|_1

是一个非光滑项，因此我们才考虑了次梯度方法。事实上我们在这里，我们只需要计算

\partial f(\beta)

，也就是计算次微分就可以了。

对光滑项计算次微分没什么好说的，因为就是简单的向量求导，《数值优化》第1节（数值优化（1）——引入，线搜索：步长选取条件）有说。所以我们可以算出

向量一范数的次梯度不难求，我们留给读者，这里直接给出结论为

X^T ( y - X \beta) = \lambda v

其中

v_i \in \begin{cases}\{1\} & \beta_i > 0 \\ \{-1\} & \beta_i < 0 \\ [-1, 1] & \beta_i = 0 \end{cases}

，

i = 1,\ldots, p

再结合一下做一个化简归类，就能得到

\begin{cases}X_i^T (y - X \beta) = \lambda \operatorname{sign} (\beta_i) & \beta_i \ne 0 \\ |X_i^T (y- X \beta) | \le \lambda & \beta_i = 0\end{cases}

这个东西的好处就是可以看出，对于LASSO来说，什么时候

\beta

的分量为0。可以看出，对于某一个特征（这里要注意一下，

X_i

指的是

的第

列）来说，如果它与估计值的残差的关系不大（用关系一词，是因为它是相关系数的分子），那么这个特征就不能算是一个重要的特征。更通俗一点来说，如果与残差关系不大，那么即使对这一个特征多加一点分量，也不会对残差的优化带来什么影响，那不如就不加了。

第二个例子其实还是LASSO，我们提它，只是因为之后的一些问题中可能会用到这个工具。

Example 2：Soft-thresholding（Source：CMU）考虑问题

\min_\beta \frac 12 \|y - \beta \|_2^2 + \lambda \|\beta \|_1

的最优性条件。

这个就是

X = I

的简化版LASSO，最优性条件也非常简单，就是

\beta = S_\lambda (y)

其中

[S_\lambda (y)]_i = \begin{cases}y_i - \lambda & y_i > \lambda \\ 0 & y_i \in [-\lambda, \lambda] \\ y_i + \lambda & y_i < -\lambda \end{cases}

i = 1, \ldots, n

你懂了Example 1，相信就能自己算出Example 2为什么是这个结果。这里要多说一句，

S_\lambda

是软阈值（soft-thresholding）算子，有些优化问题中会有一些作用。

接下来再说一个非常有意思的优化问题。

Example 3: Intersection of sets（Source：CMU）考虑问题

f_i(x) = \operatorname{dist}(x, C_i), i = 1, \ldots, m

，

f(x) = \max_{i = 1, \ldots, m}f_i(x)

，考虑优化问题

\min_x f(x)

这个问题其实就是一个寻找某一个属于一系列集合的交的元素的问题。我们可以看出，如果找到了一个点

x^*

满足

f(x^*) = 0

，那么就完成了任务。

我们注意到

\operatorname{dist}(x, C) = \min_{y \in C}\|y - x\|_2

是点到集合的最小距离，并且我们有

\nabla \operatorname{dist}(x, C) = \frac{x - P_C(x)}{\|x - P_C(x)\|_2}

其中

P_C(x)

表示

到

的投影。具体的定义可以参考《数值优化》第8节

数值优化（8）——带约束优化：引入，梯度投影法

注意区分一下标记，《数值优化》中的标记有所不同。

这个结论不是特别容易，我们首先要说明的是，如果有

\varphi_C(x) = \frac 12 \|x - P_C(x)\|_2^2

，那么

\nabla \varphi_C(x) = x - P_C(x)

。然后才能利用次梯度方法的链式法则来得到这个结论。次梯度方法的链式法则和梯度的链式法则非常类似，也是针对复合函数求次梯度使用的。这里虽然我们使用了梯度的符号，但是事实上验证这个函数是一个可微函数，次微分只有一个元素等等，也并不是一件容易的事情。碍于篇幅我们这里不细说，大家可以参考Source中的Amir (2017)中的Example 3.49。

考虑一个单独的距离函数，是因为我们的目标函数是求一系列距离函数的最大值。这里牵涉到了下面一个性质。

Lemma 1: 若

f(x) = \max_{i = 1, \ldots, m} f_i(x)

，那么

\partial f(x) = \operatorname{conv}(\cup_{i: f_i(x) = f(x)} \partial f_i(x))

。

这个乍一看可能不是特别好理解，但假如说我们考虑平面上的两个线性函数的最大值，画一张图来看，这个含义还是很清晰的。

那么这个性质就告诉我们：只要

f_i(x) = f(x), g_i \in \partial f_i(x)

，就会有

g_i \in \partial f(x)

。再多解释一些，就是如果

C_i

是距离

最远的点（否则就不是

f_i(x) = f(x)

了），那么就可以得到一个次梯度

g_i(x) = \nabla f_i(x) = \frac{x - P_{C_i}(x)}{\|x - P_{C_i}(x)\|_2}

，因此就可以使用次梯度方法了。

对于次梯度方法，我们这里使用波利亚步长（见《凸优化》第3节），这是因为波利亚步长计算简单，可以直接得到

t_k = f(x^{(k - 1)})

。这样的话，经过简单的计算，代入次梯度表达式，就可以得到

x^{(k)} = P_{C_i}(x^{(k - 1)})

简单到不可思议，当然这里的

C_i

一定要是点距离最远的集合。

下面一张图就是当集合只有2个的时候模型的结果。可以看出模型就是不停的在两个集合间去“跳转”，最终收敛到一个点，这个点属于两个集合的交集。这一个方法也有自己的名字，叫作交替投影法（Alternating Projections Algorithm）。

那么关于次梯度方法的几个案例，我们就说这么多。

随机梯度下降法

随机梯度下降法（Stochastic Gradient Descent，SGD）也是一个非常有名的机器学习和深度学习中的优化算法。它的名字带有“随机”二字，也正是因为它利用到了统计中的一些知识和概念。

随机梯度下降法一般是解决问题

\min_{x} \frac 1 m \sum_{i = 1}^m f_i(x)

如果我们使用一般的梯度下降法，可以得到迭代公式

x^{(k)} = x^{(k - 1)} - t_k \cdot \frac 1m \sum_{i =1}^m \nabla f_i(x^{(k - 1)})

可以发现，如果

取得很大，那么这意味着需要用到所有的

\nabla f_i(x^{(k -1)})

。这样的话一方面，对于一系列梯度的求和往往比较费时。因此相比较梯度下降法，随即梯度下降法考虑下面的迭代公式

x^{(k)} = x^{(k - 1)} - t_k \cdot \nabla f_{i_k}(x^{(k - 1)})

这个

\nabla f_{i_k}(x^{(k - 1)})

是啥？要说明这一点，我们要先说明我们的目标。我们之前的目标是计算梯度

\nabla f(x)

，那么如果使用随机的视角，我们就是希望找到一个替代品，这个替代品满足它的期望是

\nabla f(x)

。所以我们这里的

i_k

的选择，就要保证在指标上达成这个目标，因此选取的时候，要对每一个分量雨露均沾。这就要求

i_k

在选取的时候，必须要使得它的概率分布为均匀分布。简单来说，可以按均匀分布进行抽样，也可以说按照顺序，第一步迭代使用第一个分量函数

f_1(x)

的梯度，第二步迭代使用

f_2(x)

的梯度，以此类推循环下去。

那么关于随机梯度下降方法的步长，就没有固定步长的说法了，一般来说都考虑下降步长序列（见《凸优化》第3节）。为什么？假如说我们考虑随机梯度下降和梯度下降法，并且使用固定步长，就会得到下面两个迭代公式

x^{(k + m)} = x^{(k)} - t \sum_{i = 1}^m \nabla f_i(x^{(k + i - 1)})

x^{(k + 1)} = x^{(k)} - t \sum_{i = 1}^m \nabla f_i(x^{(k)})

如果快收敛了，这两个迭代点应该不会差距太大。然而右边的差距是

t\sum_{i = 1}^m (\nabla f_i(x^{(k + i - 1)}) - \nabla f_i(x^{(k)}))

，如果大家跑个实验的话会发现，这个式子一般是不会收敛的，因为随机梯度下降方法在接近最小值的时候，会有很大的徘徊的可能。这个具体的可以看下面这一张图。

步长说完了，我们来看看它的一些收敛性结果。在我录制的b站视频

https://www.bilibili.com/video/BV1ZK4y1b7Xt/

中，有对次梯度方法的收敛性给出证明，大体上来说，如果函数

f(x)

是一个凸函数，那么无论对于梯度是否有连续性的优化（是否Lipschitz连续），它的收敛速度都是

O(\frac 1 {\sqrt k})

，和次梯度方法是相同的。如果说有了强凸的条件，那么收敛速度会加快，但也只是

O(\frac 1 k)

。所以随机梯度下降方法虽然节省了对梯度做大量运算的时间，但是收敛所需要的迭代步数就不那么尽人意了。因此在很多时候，什么样的方法好，至少在这一方面要做一个trade-off。

当然了，对于这个方法本身，也有可以考虑优化的点。比方说小批量下降方法（Mini-batch Gradient Descent），它的思路是随机选取一个小批量的指标

I_k \subset \{1, \ldots, m\}

，且

|I_k| = b \ll m

，然后考虑迭代公式

x^{(k)} = x^{(k - 1)} - t_k \cdot \frac 1 b \sum_{i \in I_k} \nabla f_{i_k}(x^{(k - 1)})

同样的，我们希望期望就是

\nabla f(x)

，也就是

\mathbb{E} \left[ \frac 1 b \sum_{i \in I_k} \nabla f_i(x) \right] = \nabla f(x)

这个方法相比较随机梯度下降方法，因为多抽了几个梯度的方法，所以不稳定性会有减弱（在靠近极小值的时候不会反复横跳的特别厉害hhh）最好的一个收敛速率是

O(1/\sqrt{bk} + \frac 1k)

，当然这是Dekel (2012)的结果，具体现在咋样我就不知道了233。

下面这两张图是比较了不同的情况下的收敛速度，其中Full就表示正常的梯度下降法。第一张图的横坐标就是迭代步数，可以看出，在相同的迭代步数下，每一步使用到的信息越多，下降的速度就越快。但是真正的速度其实要对着第二张图看，因为第二张图暴露说，在相同的浮点运算次数下，其实反而是随机梯度下降方法最好，这潜在也说明了，其实保证随机梯度是一个原始梯度的无偏估计，在效率上还是有所保证的，毕竟浮点运算次数才是真正体现了计算机内部运算的复杂性，不然的话随机梯度下降方法明明怎么看收敛速度都不如梯度下降法，为什么反而在大规模问题中是主流算法呢？

再下一张图是观察次梯度方法与梯度方法在迭代进行过程中的极小值精度。我们上一节（也就是这一节反复出现过的《凸优化》第3节）有提到说，次梯度方法在精度上其实是比较一般的。而在随机梯度下降方法中，这个现象也得到了传承，并且可以看出，信息越多，差距越小。不过这也只是经验性的判断，不能完全的当成一个准则。

大规模问题下的随机梯度下降方法使用技巧

这一部分会非常的不数学，因为有非常多的依赖经验的技巧。这些技巧多见于机器学习，深度学习中。虽说不够数学，但它们都有一些相对比较合理的解释，所以倒也不算是空穴来风。

先说一个比较次要的。

Trick 1: 在使用随机梯度下降进行大规模数据的训练时，可以先在一个小范围的数据集内训练，实验，调参。

这一句话看似是一个废话，但是也是有很多可以说的点。因为实际的数据集之大，往往是撑不住一次训练的时间成本的，所以我们要考虑随机筛选一部分足够泛化的数据。这主要还是为了参数的选取，因为如果我们选的泛化性很好，那么在真正的大规模的数据集下，就没有必要再考虑调参这一步了。

Trick 2: 若使用优化算法解决机器学习中的优化问题，那么一般使用固定步长。

这个要看情况来说，比方说SGD就是一个反例。但如果说去掉SGD，为什么其他方法推荐用固定步长呢？这是因为，对于机器学习来说，它更关心最优解带来的性质，那么这个最优解误差0.01，和误差0.000001，有什么本质区别呢？其实是没有的。这也是为什么次梯度方法虽然有精度的劣势，但是在机器学习中会被应用的原因。

Trick 3: 优化方法无处不在，建议积攒一些调参经验。

emmmm……这个我相信搞深度学习的会很有感触。比方说深度学习中有相当多的优化器，比如说Adam，比如说SVRG，当然SGD自己就是一个被广为使用的优化器。但是究竟对于各种问题中，什么样的形式是效果最好的，这是一个很难回答的问题，可以说没有什么定数，全看命了。

Trick 4: 对于一些带约束的优化问题，提前终止是一个好的方法。

这个可能有点难理解，我们拿逻辑回归来举例子。考虑带罚项的逻辑回归问题

\min_\beta \frac 1n \sum_{i=1}^n (-y_i x_i^T \beta + \log(1 + e^{x_i ^T \beta})) \quad s.t. \|\beta \|_2 \le t

根据优化问题对偶性，这个和带罚项的问题的解是一致的。关于对偶性会在《凸优化》之后的章节重点提到，这里大家知道这个结论就可以了。

事实上，我们可以忽略这个约束，只对目标函数进行一个优化，然后提前终止。这个阴间操作真的能行吗？答案可能会有些出乎意料，因为下面这张图确确实实说明了，两种方法对应的得到的不同的参数的解是非常相近的。图中纵坐标就是分量值，横坐标在这个问题中就是约束值

。

所以我们在实际操作的时候，我们可以观察我们函数值

\beta

。如果

\beta

一直很小，那没什么好说的，但等到

\beta

差不多到达

，或者要超过

的时候，就考虑终止它。这样做的理论基础尚未建立，但有一个经验化的解释是：在梯度下降方法中，我们每一步都会降低函数值，但同时会增大自变量值。这个行为其实和罚项所希望满足的事情是相似的。

那么到此为止，关于梯度方法，我们就算说完了。

近端梯度方法引入

你以为到此就结束了？8太可能的，我这字数的kpi还没到呢……还剩一点字数，我们简单说一下近端梯度方法（Proximal Gradient Method）的设计原理。

在这一节我们提到说，对于次梯度方法，它的最好的收敛速度也只是

O(\frac 1 {\sqrt k})

。因此如果再考虑所有可能的问题，也就没有再进一步改进算法的意义。因此近端梯度方法只关心一个特定的问题，也就是

f(x) =g(x) + h(x)

这里只要求

g(x)

是一个凸和可微的函数，

h(x)

虽然要求凸，但是它可以不光滑，另外我们希望它尽量简单。

近端梯度方法的思路其实也来源于梯度下降法。注意到梯度下降法的更新公式为

x^+ = x - t \nabla f(x)

但是它同时也有表达式

x^+ = \arg \min_{z}f (x) + \nabla f(x)^T ( z -x ) + \frac 1 {2t} \|z - x \|_2^2

如果我们设

\bar f_t(z) = f (x) + \nabla f(x)^T ( z -x ) + \frac 1 {2t} \|z - x \|_2^2

，那么我们面对这个求和的问题，自然会想到，如果不考虑

，只考虑对

进行这样的替换，那么迭代公式会变成什么？所以我们考虑下面这样的变换

x^+ = \arg \min_z \bar g_t(z) + h(z)

然后考虑变换

x^+ = \arg\min_{z} \bar g_t(z) + h(z) = \arg\min_z g(x) + \nabla g(x)^T (z -x ) + \frac 1 {2t}\|z - x\|_2^2 + h(z)

= \arg\min_z \frac 1{2t} \|z - (x - t\nabla g(x))\|_2^2 + h(z)

注意到变换中，我们只需要求一个关于

的极小值。所以配方前后可能式子的值不一样，但是关于

的那一部分一定是没变的。

所以说现在我们将这个式子拆分成为了两个部分。第一个部分是

\frac 1{2t} \|z - (x - t\nabla g(x))\|_2^2

，表示希望下一个迭代点

与

的梯度下降的更新结果差距尽可能的小。另外一方面，还希望

h(z)

能够尽可能的小。这也就是近端梯度方法希望解决的问题。

那么至于如何解决呢？我们就放到下一节再继续说吧。

小结

这一节我们走马观花的介绍了很多算法的延伸和一些经验性的实操技巧。它们都或多或少的与机器学习，深度学习等等有些关系。当然了这些技巧和实战经验都有一些还算合理的解释和实操方案，了解这些也会对各种有趣的优化算法有更加深层次的理解。

在下一节我们会继续介绍近端梯度方法。如果篇幅允许，我们会开始再回到理论，介绍运筹学中的一些对偶理论。

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2020-09-12，如有侵权请联系 cloudcommunity@tencent.com 删除

机器学习

神经网络

深度学习

人工智能

本文分享自学弱猹的精品小屋微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度