用于非负最小二乘回归的约束贪婪随机化高斯-塞德尔方法

《Journal of Computational and Applied Mathematics》：Constrained greedy randomized Gauss-Seidel method for nonnegative least-squares regression

　　

## 摘要
为求解带非负约束的大规模线性回归问题，我们提出了一种基于贪婪概率准则的随机迭代投影方法。进行了理论分析并开展了数值实验。数值结果表明，我们的方法在迭代步数和计算时间上均显著优于现有最先进的方法。

## 引言
本文考虑如下带非负约束的大规模线性系统
$$Ax=b, \quad \text{s.t. } x \ge 0,$$
其中 $A \in \mathcal{R}^{m \times n}$ 是系数矩阵，$b \in \mathcal{R}^m$ 是右端向量，$x \in \mathcal{R}^n$ 是未知非负向量。由于非负约束，线性系统 (1.1) 实际上通常是不相容的。此时，我们可以转而考虑如下非负最小二乘回归问题
$$\min_{x \ge 0} f(x) = \min_{x \ge 0} \frac{1}{2} \|b - Ax\|_2^2,$$
其中 $\|\cdot\|_2$ 表示向量的欧几里得范数。线性系统 (1.1) 或优化问题 (1.2) 出现在广泛的应用中，例如图像重建 [1], [2] 和非负学习 [3], [4]，并在过去几十年中受到了广泛关注。由于其在实际数据集上的优异表现，许多用于求解其解的数值确定性方法 [5], [6], [7], [8], [9], [10] 已相继被开发出来。其中，投影梯度 (PG) 方法 [7] 可能最具吸引力。

PG 的迭代形式最初由 Goldstein [11] 以及 Levitin 和 Polyak [12] 提出，形式如下
$$x^{k+1} = \max(x^k - \alpha_k g^k, 0),$$
其中 $x^k$ 表示第 $k$ 个近似解，$\alpha_k \ge 0$ 表示搜索步长，$g^k = A^T(Ax^k - b)$ 表示目标函数 $f(x)$ 在 $x^k$ 处的梯度，$\max(\cdot, 0)$ 表示非负投影，该投影按分量定义。PG 的关键在于高效地选择合适的步长 $\alpha_k$，使得由 (1.3) 生成的向量序列 $\{x^k\}$ 的任何极限点 $x^*$ 满足优化问题 (1.2) 的一阶必要条件，即
$$\begin{cases} g^*(i) \ge 0, & \text{if } x^*(i) = 0, \\ g^*(i) = 0, & \text{if } x^*(i) > 0, \end{cases} \quad \text{对于任意 } i \in [n] = \{1, 2, \dots, n\},$$
其中 $g^*(i)$ 和 $x^*(i)$ 分别表示向量 $g^* = g(x^*) = A^T(Ax^* - b)$ 和 $x^*$ 的第 $i$ 个分量。可以轻易证实，当且仅当 $p^* = p(x^*) = 0$ 时，$x^*$ 满足条件 (1.4)，其中 $p(x) = \min(Dx, g(x))$ 按分量定义，且 $D$ 是任何具有正对角线元素的对角矩阵。

为高效选择 $\alpha_k$，Bertsekas [5] 提出了一种名为广义 Armijo 步长规则的实际方法。对于任意给定的 $\beta, \sigma \in (0, 1)$ 和 $\gamma > 0$，步长 $\alpha_k$ 被选为 $\alpha_k = \beta \gamma^{l_k}$，其中 $l_k$ 是第一个满足以下条件的非负整数 $l$：
$$f_k - f_{k+1} \ge \sigma g_k^T(x^k - x^{k+1}),$$
其中 $f_k = f(x^k), f_{k+1} = f(x^{k+1}), x^{k+1}$ 由公式 (1.3) 计算得到，且 $\alpha_k = \beta \gamma^{l_k}$，上标 $T$ 表示对应矩阵或向量的转置。然而，如 (1.5) 所述确定合适的步长 $\alpha_k$ 是耗时的，特别是对于大规模系统，因为我们需要在每次迭代中不断评估几次矩阵 - 向量乘积来计算 $f_{k+1}$ 以测试规则 (1.5)。

据我们所知，以 Kaczmarz 方法 [13] 为代表的行投影迭代方法（也被称为代数重建技术 (ART)）对于大规模无约束线性系统非常有效，因为它们不需要访问整个数据集，而在每次迭代中仅需要访问单行或几行。当线性系统 $Ax=b$ 是相容的时，Kaczmarz 型方法具有以下形式
$$x^{k+1} = x^k + \frac{b(i_k) - A(i_k, :)x^k}{\|A(i_k, :)\|_2^2} A(i_k, :)^T,$$
其中 $A(i_k, :)$ 表示 $A$ 的第 $i_k$ 行，且 $i_k$ 从行索引集 $[m]$ 中按特定准则选择。在经典的 Kaczmarz [13] 方法中，$i_k$ 在第 $k$ 步被选为 $(k \mod m) + 1$。

基于这一点，Popa [1] 为约束线性系统 (1.1) 提出了一种约束 Kaczmarz 扩展 (CKE) 方法，其形式如下
$$\begin{cases} z^k = z^k - \frac{A(:, j)^T z^k}{\|A(:, j)\|_2^2} A(:, j), & j = 1, 2, \dots, n, \\ x^k = x^k + \frac{b(i) - z^k(i) - A(i, :)x^k}{\|A(i, :)\|_2^2} A(i, :)^T, & i = 1, 2, \dots, m, \\ x^k = \max(x^k, 0), \\ k := k + 1, \end{cases}$$
其中 $z^0 = b$，且 $A(:, j)$ 表示 $A$ 的第 $j$ 列。可以注意到，(1.7) 中的前两个方程对应于 Kaczmarz 扩展过程 [14]，因为系统 (1.1) 通常是不相容的。

代替给定的循环顺序选择 $i_k$，Strohmer 和 Vershynin [15] 建议随机选择 $i_k$，概率为 $\|A(i, :)\|_2^2 / \|A\|_F^2$，其中 $\|\cdot\|_F$ 表示 Frobenius 矩阵范数，并提出了随机 Kaczmarz (RK) 方法。由于 RK 可以大大提高经典 Kaczmarz 方法的收敛率，Wu 和 Xiang [16] 为约束线性系统 (1.1) 提出了一种投影随机扩展 Kaczmarz 方法 (PREK)。PREK 的形式可以表示为
$$\begin{cases} z^{k+1} = z^k - \frac{A(:, j_k)^T z^k}{\|A(:, j_k)\|_2^2} A(:, j_k), \quad \Pr(j_k) = \frac{\|A(:, j_k)\|_2^2}{\|A\|_F^2}, \\ x^{k+1} = x^k + \frac{b(i_k) - z^{k+1}(i_k) - A(i_k, :)x^k}{\|A(i_k, :)\|_2^2} A(i_k, :)^T, \quad \Pr(i_k) = \frac{\|A(i_k, :)\|_2^2}{\|A\|_F^2}, \\ x^{k+1} = \max(x^{k+1}, 0), \end{cases}$$
其中 $z^0 = b$。

然而，RK 在选择工作行 $i_k$ 时存在一个明显的弱点，即如果 $A$ 被一个合适的对角矩阵缩放，使得其所有行的欧几里得范数归一化为相同的常数，则概率准则会退化为均匀采样。为克服这一弱点，Bai 和 Wu [17] 设计了一种贪婪概率准则，旨在每次迭代中抓住残差向量 $r^k = b - Ax^k$ 的较大分量，并提出了贪婪随机 Kaczmarz (GRK) 方法。具体来说，如果 $|r^k(i)| > |r^k(j)|$，则希望第 $i$ 行被选中的概率大于第 $j$ 行，以便尽可能优先消除 $r^k$ 的较大分量。回到线性系统 (1.1) 或优化问题 (1.2)，根据条件 (1.4)，我们显然希望 $p^k = \min(Dx^k, g^k)$ 的较大分量（而不是 $r^k$）能被优先尽可能消除，这类似于贪婪随机坐标下降 (GRCD) 方法 [18] 中处理 $g^k$ 而不是 $r^k$ 的情况，该方法也被称为贪婪随机 Gauss-Seidel (GRGS) 方法。受此思想启发，本文提出了一种约束贪婪随机 Gauss-Seidel (CGRGS) 方法，用于求解线性回归问题 (1.2) 的非负解。

本文其余部分组织如下。在第 2 节中，我们提出了约束贪婪随机 Gauss-Seidel 方法。在第 3 节中，我们对 CGRGS 方法进行了理论分析。在第 4 节中，我们实施了一些数值实验来说明所提出方法的表现。最后，我们在第 5 节中给出了一些结论以结束本文。

## 2 约束贪婪随机 Gauss-Seidel 方法
注意，通过选择不同的对角矩阵 $D$，我们可以得到不同的 $p^k$。本文中，我们将其选为 $D = \text{diag}(d_1^2, d_2^2, \dots, d_n^2)$，其中 $d_j = \|A(:, j)\|_2$，则 CGRGS 的形式可以表示为
$$x^{k+1}(j_k) = x^k(j_k) - \frac{1}{d_{j_k}^2} p^k(j_k),$$
其中 $p^k = \min(Dx^k, g^k)$，$g^k = -A^T r^k$，$r^k = b - Ax^k$，且工作列 $j_k$ 由以下贪婪概率准则选择
$$\Pr(j_k) = \frac{\tilde{p}_k^2(j_k)}{\|\tilde{p}_k\|_2^2},$$
其中
$$\begin{cases} \tilde{p}_k(j) = \begin{cases} p^k(j), & \text{if } j \in U_k, \\ 0, & \text{if } j \notin U_k, \end{cases} \\ U_k = \{j | p^2_k(j) \ge \delta_k \|p_k\|_2^2 d_j^2 \}, \\ \delta_k = \theta \max_{j \in [n]} \frac{p^2_k(j)}{d_j^2} \|p_k\|_2^2 + (1-\theta)\|A\|_F^2, \end{cases}$$
其中 $\theta \in [0, 1]$ 是一个参数。

## 3 理论分析
我们首先指出，根据贪婪概率准则 (2.2)，$p^k$ 的较大分量确实被优先尽可能消除了。

**引理 3.1** 如果 $j_k$ 是第 $k$ 步的工作列，则对所有 $k \ge 0$，都有 $p^{k+1}(j_k) = 0$ 成立。

**证明** 可以轻易发现 $p^{k+1}(j_k) = \min(d_{j_k}^2 x^{k+1}(j_k), g^{k+1}(j_k))$ 且 $g^{k+1}(j_k) = g^k(j_k) - p^k(j_k)$。然后我们在以下两种情况下证明 $p^{k+1}(j_k) = 0$。

**情况 I.** 如果 $x^{k+1}(j_k) = 0$，我们从公式 (2.1) 得到 $p^k(j_k) = d_{j_k}^2 x^k(j_k) \ge 0$，这意味着 $g^k(j_k) \ge p^k(j_k$。

## 4 数值实验
在本节中，我们将进行一些数值实验来展示我们所提出方法的有效性，即 CGRGS (算法 2.1) 用于求解非负最小二乘回归问题 (1.2)。作为比较，我们也给出了 CKE、PREK 和 PG 的数值行为。这里的数值行为主要指迭代步数，用"IT"表示，以及以秒为单位的计算时间，用"CPU"表示。此外，我们还可以比较函数值 $f_k$，用

## 5 结论
本文提出了一种约束贪婪随机 Gauss-Seidel 方法，用于求解大规模非负最小二乘问题 (1.2)，该方法旨在根据其最优条件 (1.4) 优先尽可能消除 $p^k$ 的较大分量。我们分析了 CGRGS 所有样本路径迭代的极限点平稳性性质，并展示了其数值行为以及 CKE、PREK 和 PG 的表现。实验结果表明，CGRGS 可以优于 CKE 和

## 利益冲突声明
作者声明他们没有任何已知的竞争性经济利益或个人关系，这些可能影响本文报告的工作。

## 致谢
本研究由中华人民共和国国家自然科学基金（批准号：12201302）和江苏省高等教育机构自然科学基金（批准号：21KJB110017）资助。

Wei Wei|Ming Yuan|Xiaoling Li|Xi Yang
南京林业大学应用数学系，南京，210037，中华人民共和国