凸优化实战1

近期工作主要在关注于优化问题求解，积累了一些常用套路

线性回归

设样本数 $n$, 样本特征矩阵 $X \in \mathbb{R}^{n \times m}$, 样本标签矩阵 $Y \in \mathbb{R}^{n \times q}$, 线性回归参数矩阵 $W \in \mathbb{R}^{m \times q}$.

则最基础的线性回归问题有如下形式

$$
\min_W \frac{1}{2}\|Y-XW\|_F^2
$$

由于函数是一个凸函数 (证明略），肯定有全局最优解。

令

$$
\nabla \ell(W) = X^T (XW-Y) = 0
$$

不难得到

$$
W^* = (X^T X)^{-1} X^T Y
$$

事实上，由于要考虑 $X$ 矩阵病态不可求逆的情况，
数值计算时大多数时候都会用 $\alpha$ 特别小（$1\times 10^{-8}$ ）的 Rigid Regression
来代替朴素的 Linear Regression。

F范数正则化 （Rigid Regression）

对$W$ 添加一个F范数正则项（正则项系数 $\alpha \in \mathbb{R}$），则有如下形式（也称岭回归问题）

$$
\min_W \frac{1}{2}\|Y-XW\|_F^2 + \frac{\alpha}{2} \|W\|_F^2
$$

函数是凸函数加凸函数，还是凸函数，依然有全局最优解。继续用上面的思路去找到闭式解。

令

$$
\nabla \ell(W) = X^T(XW-Y) + \alpha W = 0
$$

不难得到

$$
W^* = (X^T X + \alpha I)^{-1} X^T Y
$$

L1范数正则化 （Lasso Regression）

$$
\min_W \frac{1}{2} \|Y-XW\|_F^2 + \alpha \|W\|_1
$$

L1范数也是凸函数，所以这个问题依然有全局最优解。
然而，L1范数不好求导，你得不到闭式解。

已知，这个世界上有个技术叫Shrinkage Operator，它专克 L1范数。

记软阈值算子 $\mathcal{S}_\tau(\cdot)$ （soft-thresholding operator)

$$
\begin{aligned}
\mathcal{S}_\tau(x) & = \left\{
\begin{aligned}
& x - \tau, \quad & x > & \ \tau \\
& 0, \quad & |x| < & \ \tau \\
& x + \tau, \quad & x < & -\tau \\
\end{aligned} \right. \\
& = \min(x + \tau, 0) + \max(x - \tau, 0)
\end{aligned}
$$

对于

$$
\min_X \frac{1}{2} \|X-Y\|_F^2 + \alpha \|X\|_1
$$

其最优解为

$$
X^* = \mathcal{S}_\alpha(Y)
$$

那么，为了把这条引理利用起来，我们将原问题进行松弛

$$
\begin{aligned}
\min_{W,P} \quad & \frac{1}{2} \|Y-XW\|_F^2 + \alpha \|P\|_1 \\
s.t. \quad & P = W
\end{aligned}
$$

用增广拉格朗日函数（ADM）改写损失函数, 引入惩罚项消除等式约束
（其中 $\rho$ 是超参数, 一般在 $(0, 1)$ 之间取值）

$$
\min_{W,P,U} \quad \frac{1}{2} \|Y-XW\|_F^2 + \alpha \|P\|_1 +
\frac{\rho}{2} \|P-W+U\|_F^2
$$

然后使用ADMM算法以如下步骤迭代优化

1. $W \leftarrow \arg \min \ell(W) = \frac{1}{2} \|Y-XW\|_F^2$
2. $P \leftarrow \arg \min \ell(P) = \frac{\rho}{2} \|P-W+U\|_F^2 + \alpha \|P\|_1$
3. $U \leftarrow U + (P-W)$

$W$ 的更新就不赘述了，可以看到这里 $P$ 的更新恰好可以用上Shrinkage Operator的引理

$$
\begin{aligned}
\arg\min_P & \quad \frac{1}{2} \|P-(W-U)\|_F^2 + \frac{\alpha}{\rho} \|P\|_1 \\
= & \quad \mathcal{S}_{\alpha/\rho}(W-U)
\end{aligned}
$$

综上，这个算法以如下步骤反复迭代更新

1. $W \leftarrow X^T(XW-Y)$
2. $P \leftarrow \min(W-U+\alpha/\rho, 0) + \max(W-U-\alpha/\rho, 0)$
3. $U \leftarrow U + (P-W)$

直到 $W$ 和 $P$ 的更新幅度足够小，例如

$$
\frac{\|\Delta W\|_F^2}{\|W\|_F^2} < 1 \times 10^{-4}
$$

核范数正则化

$$
\min_W \frac{1}{2} \|Y-XW\|_F^2 + \alpha \|W\|_*
$$

核范数的定义为 矩阵奇异值的和

$$
\begin{aligned}
\|X\|_* & = \mathrm{sum}(\Sigma) \\
s.t. \quad & X = U \Sigma V^T
\end{aligned}
$$

其中 $\Sigma$ 为 $X$ 的奇异值。

或者, 矩阵 $X$ 的核范数也可以定义为

$$
\|X\|_* = \mathrm{tr}\bigl((X^T X)^{1/2}\bigr) \qquad
$$

其中 $(\cdot)^{1/2}$ 为矩阵的平方根分解, 即, 对于有特征值分解的矩阵 $A = P \Lambda P^{-1}$, 其平方根分解为 $A^{1/2} = P \Lambda^{1/2} P^{-1}$.

(感谢吕文龙同学的勘误!)

类似的，也有一个叫奇异值收缩（Singular Value Thresholding）的引理可以帮助我们解决这个问题。

记奇异值收缩算子 $\mathcal{D}_\tau(\cdot)$ (Singular Value Shrinkage Operator)

$$
\begin{aligned}
\mathcal{D}_\tau(X) & = U \tilde{\Sigma} V^T \\
s.t.\quad X & = U \Sigma V^T \\
\quad \tilde{\Sigma} & = \max(\Sigma - \tau, 0)
\end{aligned}
$$

对于

$$
\min_X \frac{1}{2} \|X-Y\|_F^2 + \alpha \|X\|_*
$$

其最优解为

$$
X^* = \mathcal{D}_{\alpha}(Y)
$$

类似地，我们引入松弛变量 $P$, 用增广拉格朗日函数(ADM)消除等式约束

$$
\min_{W,P,U} \quad \frac{1}{2} \|Y-XW\|_F^2 + \alpha \|P\|_* +
\frac{\rho}{2} \|P-W+U\|_F^2
$$

然后使用ADMM算法以如下步骤迭代优化

1. $W \leftarrow \arg \min \ell(W) = \frac{1}{2} \|Y-XW\|_F^2$
2. $P \leftarrow \arg \min \ell(P) = \frac{\rho}{2} \|P-W+U\|_F^2 + \alpha \|P\|_*$
3. $U \leftarrow U + (P-W)$

代入奇异值收缩方法，得到如下更新公式

1. $W \leftarrow X^T(XW-Y)$
2. $P \leftarrow \mathcal{D}_{\alpha/\rho}(W-U)$
3. $U \leftarrow U + (P-W)$

直到 $W$ 和 $P$ 的更新幅度足够小

低秩表示（Low-Rank Representaiton, LRR）

给定数据矩阵 $X \in \mathbb{R}^{n \times n}$ 和字典矩阵 $A \in \mathbb{R}^{n \times d}$，
低秩表示旨在求解如下优化问题

$$
\begin{aligned}
\min_{Z,E} \quad &
\|Z\|_* + \lambda \|E\|_1 \\
s.t. \quad &
X = A Z + E
\end{aligned}
$$

得到低秩表示系数 $Z \in \mathbb{R}^{d \times m}$ 和噪声 $E \in \mathbb{R}^{n \times m}$。

引入松弛变量 $J$, 使得（下一步增广之后）单独优化 $Z$ 时关于 $Z$ 损失函数足够简单.

$$
\begin{aligned}
\min_{Z, E, J} \quad & \|J\|_* + \lambda \|E\|_1 \\
s.t. \quad & X = A Z + E \\
\quad & Z = J
\end{aligned}
$$

使用ADM对目标函数进行增广，消除等式约束.

$$
\begin{aligned}
\min_{Z, E, J, L_1, L_2}
\quad & \|J\|_* + \lambda \|E\|_1 + \\
\quad & \mathrm{tr}\bigl( L_1^T (AZ+E-X) \bigr) + \\
\quad & \mathrm{tr}\bigl( L_2^T (AZ+E-X) \bigr) + \\
\quad & \frac{\rho}{2} \|AZ+E-X\|_F^2 + \frac{\rho}{2} \|J-Z\|_F^2 \\
\end{aligned}
$$

设

$$
U_1 = \frac{1}{\rho} L1 \\
U_2 = \frac{1}{\rho} L2
$$

代入后，ADMM可以得到更简洁的形式（Scaled Form）

$$
\begin{aligned}
\min_{Z, E, J, U_1, U_2}
\quad & \|J\|_* + \lambda \|E\|_1 + \frac{\rho}{2} \|AZ+E-X+U_1\|_F^2 + \\
\quad & \frac{\rho}{2} \|J-Z+U_2\|_F^2 \\
\end{aligned}
$$

然后以如下步骤迭代优化

1. $Z \leftarrow \arg \min \ell(Z)$
2. $E \leftarrow \arg \min \ell(E)$
3. $J \leftarrow \arg \min \ell(J)$
4. $U_1 \leftarrow U_1 + (AZ+E-X)$
5. $U_2 \leftarrow U_2 + (J-Z)$

每个公式都能得到闭式解，稍微推一下就能得到如下内容

1. $Z \leftarrow \bigl( A^TA+I \bigr)^{-1} \bigl( A^T(X-E-U_1)+J+U_2 \bigr)$
2. $E \leftarrow \mathcal{S}_{\lambda/\rho}(X-AZ-U_1)$
3. $J \leftarrow \mathcal{D}_{\lambda/\rho}(Z-U_2)$
4. $U_1 \leftarrow U_1 + (AZ+E-X)$
5. $U_2 \leftarrow U_2 + (J-Z)$

就这么迭代吧！

附：参考代码（matlab）

% low rank representation
function [rep, noise] = low_rank_representation(data, dict, varargin)
    % given 
    %       X - n x m matrix of data (data)
    %       A - n x d matrix of dictionary (dict)
    % low-rank represetation problem aims to solve optimization problem
    %       min_{Z, E}  |Z|_* + lambda |E|_1
    %       s.t.        X = AZ + E
    % where
    %       Z - d x m matrix of representation coefficient (rep)
    %       E - n x m matrix of noise (noise)
    [n, m] = size(data);
    [n, d] = size(dict);
    assert(size(data, 1) == size(dict, 1));

    args = inputParser;
    addParameter(args, 'lambda', 0.1);      % smaller lambda will extract more noise
    addParameter(args, 'rho', 0.1);         % ADMM hyper-parameter, in (0, 1)
    addParameter(args, 'max_iter', 200); 
    addParameter(args, 'tol', 1e-4);
    addParameter(args, 'show_progress', false);
    parse(args, varargin{:});

    lambda        = args.Results.lambda;
    rho           = args.Results.rho;
    max_iter      = args.Results.max_iter;
    tol           = args.Results.tol;
    show_progress = args.Results.show_progress;

    % by applying exact ADM method on the raw optimization problem, 
    % we have
    %       min_{Z, E, J, U1, U2}   |J|_* + lambda |E|_1
    %                               + rho/2 |AZ+E-X+U1|_F^2
    %                               + rho/2 |J-Z+U2|_F^2 + const, 
    % where
    %       J  - d x m matrix of slack variable
    %       U1 - n x m matrix of lagrange multiplier
    %       U2 - d x m matrix of lagrange multiplier.
    % and algorithm can be described as follow
    %       Z  <- (A'A+I)\(A'(X-E-U1)+(J+U2))           (1)
    %       E  <- soft(X-AZ-U1, lambda/rho)             (2)
    %       J  <- svt(Z-U2, 1/rho)                      (3)
    %       U1 <- U1 + (AZ+E-X)                         (4)
    %       U2 <- U2 + (J-Z)                            (5), 
    % where
    %       soft(X, a) = max(X - a, 0) + min(X + a, 0)
    %       svt(X, a)  = U S^ V'
    %                    (s.t. USV' = X, S^ = max(S - a, 0)).

    function loss = get_loss()
        l1   = sum(svd(Z));
        l2   = lambda * sum(abs(E), 'all');
        loss = [l1, l2];
    end

    X  = data; % n x m
    A  = dict; % n x m
    Z  = zeros(d, m);
    E  = zeros(n, m);
    J  = zeros(d, m);
    U1 = zeros(n, m);
    U2 = zeros(d, m);

    norm_X = norm(X, 'fro'); % relative number for tolerance setting

    if show_progress
        fprintf('     loss       |Z|_*      l|E|_1     soft-lim\n');
    end

    for t0 = 1:max_iter
        for t1 = 1:max_iter
            Z_last = Z;
            E_last = E;

            % update Z
            %   Z  <- (A'A+I)\(A'(X-E-U1)+(J+U2))           (1)
            Z = (A'*A+eye(m))\(A'*(X-E-U1)+(J+U2));

            % update E
            %   E  <- soft(X-AZ-U1, lambda/rho)             (2)
            E_hat = X-A*Z-U1;
            E = max(E_hat - lambda/rho, 0) + min(E_hat + lambda/rho, 0);

            % update J
            %   J  <- svt(Z-U2, 1/rho)                      (3)
            J_hat = Z-U2;
            [U_, Sigma_, V_] = svd(J_hat, 'econ');
            Sigma_ = max(Sigma_ - 1/rho, 0);
            J = U_ * Sigma_ * V_';

            % converage condition
            delta_Z = norm(Z - Z_last, 'fro');
            delta_E = norm(E - E_last, 'fro');
            if (delta_Z/norm_X < tol) && (delta_E/norm_X < tol)
                break;
            end
        end

        % update U1, U2
        %   U1 <- U1 + (AZ+E-X)                         (4)
        %   U2 <- U2 + (J-Z)                            (5)
        U1 = U1 + (A*Z+E-X);
        U2 = U2 + (J-Z);

        % if mannually scale up rho, it converage much faster
        rho = min(6.0 * rho, 1e10);

        if show_progress
            loss = get_loss();
            fprintf('%2d:  %8.2e   %8.2e   %8.2e   %8.2e\n', ..., 
                    t0, sum(loss), loss(1), loss(2), norm(A*Z+E-X, 'fro')^2;);
        end

        % converage condition
        error1 = norm(X-A*Z-E, 'fro');
        error2 = norm(Z-J, 'fro');
        if (error1/norm_X < tol) && (error2/norm_X < tol)
            break;
        end
    end

    rep   = Z;
    noise = E;
end

（还是讨厌matlab，讨厌讨厌讨厌！）

---

感谢一言不发的大野和笨蛋春雄！
封面来自押切莲介老师的《高分少女》。