【凸优化】邻近点梯度法、交替方向乘子法、次梯度法 求解稀疏矩阵问题 python实现

1 问题描述

考虑一个包含 10 个节点的分布式系统。每个节点 i i i 有线性测量如下:

b i = A i x + e i b_i = A_i x + e_i bi=Aix+ei

其中:

  • b i b_i bi 是一个 5 维的测量值
  • A i A_i Ai 是一个 5 × 200 维的测量矩阵
  • x x x 是一个 稀疏度为 5 的 200 维的未知稀疏向量
  • $e_i $是 一个 5 维的测量噪声

我们需要从所有 b i b_i bi A i A_i Ai 中恢复 x x x,其目标函数为:

m i n x ( ∥ A 1 x − b 1 ∥ 2 2 + ⋯ + ∥ A 10 x − b 10 ∥ 2 2 + λ ∥ x ∥ 1 ) min_x \left( \|A_1x - b_1\|_2^2 + \dots + \|A_{10}x - b_{10}\|_2^2 + \lambda \|x\|_1 \right) minx(A1xb122++A10xb1022+λx1)

其中 λ λ λ > 0 是正则化参数。请设计下述算法求解该问题:

  1. 邻近点梯度法(Proximal Gradient Method)
  2. 交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)
  3. 次梯度法(Subgradient Method)

在实验中,假设:

  • x 的真实值中的非零元素服从均值为 0,方差为 1 的高斯分布。
  • A i A_i Ai 中的元素服从均值为 0,方差为 1 的高斯分布。
  • e i e_i ei 中的元素服从均值为 0,方差为 0.1 的高斯分布。

对于每种算法,给出以下结果:

  1. 每步计算结果与真值的距离
  2. 每步计算结果与最优解的距离。
  3. 讨论正则化参数 λ λ λ 对计算结果的影响。

2 算法设计与数值实验

m i n x ( ∥ A 1 x − b 1 ∥ 2 2 + ⋯ + ∥ A 10 x − b 10 ∥ 2 2 + λ ∥ x ∥ 1 ) min_x \left( \|A_1x - b_1\|_2^2 + \dots + \|A_{10}x - b_{10}\|_2^2 + \lambda \|x\|_1 \right) minx(A1xb122++A10xb1022+λx1)

  • 我们定义一个类,所有算法实现都放在方法里:

    class OptimizedAlgorithm:...
    
  • 依照题目要求初始化数据:

    def __init__(self):
            """
            self.num = 10  # 节点数
            self.ASize = (5, 200)  # 测量矩阵大小
            self.XSize = 200  # x向量维度
            self.X # x向量
            self.A = np.array(all_As) # 十个节点的A矩阵数组
            self.b = np.array(all_bs) # 十个节点的b矩阵数组
            """
            # 初始化数据
            self.num = 10  # 节点数
            self.ASize = (5, 200)  # 测量矩阵大小
            self.XSize = 200  # x向量维度
            # 生成x向量
            sparsity = 5
            self.X = np.zeros(self.XSize)
            XIndex = random.sample(list(range(self.XSize)), sparsity)
            for index in XIndex:
                self.X[index] = np.random.randn()  # x向量
            # 存放所有矩阵的列表
            all_As = []
            all_bs = []
            # 生成Ai ei bi
            for _ in range(self.num):
                Ai = np.random.normal(0, 1, self.ASize)
                ei = np.random.normal(0, 0.1, 5)
                bi = np.dot(Ai, self.X) + ei
                all_As.append(Ai)
                all_bs.append(bi)
    
            self.A = np.array(all_As)  # 十个节点的A矩阵数组
            self.b = np.array(all_bs)  # 十个节点的b矩阵数组
    

  • 邻近点梯度法

    1. 算法设计

      对于复合凸优化问题 m i n x ( f 0 ( x ) : = g ( x ) + r ( x ) ) min_x(f_0(x):=g(x)+r(x)) minx(f0(x):=g(x)+r(x)),其中:

      g ( x ) = ( ∥ A 1 x − b 1 ∥ 2 2 + ⋯ + ∥ A 10 x − b 10 ∥ 2 2 ) , r ( x ) = λ ∣ ∣ x ∣ ∣ 1 g(x)=\left( \|A_1x - b_1\|_2^2 + \dots + \|A_{10}x - b_{10}\|_2^2 \right),r(x)=\lambda||x||_1 g(x)=(A1xb122++A10xb1022),r(x)=λ∣∣x1

      定义邻近算子:

      p r o x α r ( x ) = argmin u ( r ( u ) + 1 2 α ∣ ∣ u − x ∣ ∣ 2 2 ) prox_{\alpha r}(x) = \text{argmin}_u \left( r(u) + \frac{1}{2\alpha} ||u - x||_2^2 \right) proxαr(x)=argminu(r(u)+2α1∣∣ux22)

      1. 针对光滑部分 g g g 进行梯度下降更新:

        x k + 1 2 = x k − α k ∗ ∇ g ( x k ) x^{k+{1\over2}} = x^k - \alpha_k * ∇g(x^k) xk+21=xkαkg(xk)

      2. 对非光滑部分 r r r 使用邻近算子转换:

        x k + 1 = p r o x α k r ( x k + 1 2 ) x^{k+1} = prox_{\alpha_k r}(x^{k+{1\over2}}) xk+1=proxαkr(xk+21)

      3. 将1,2步结合起来得到迭代式:

        x k + 1 = p r o x α k r ( x k − α k ∗ ∇ g ( x k ) ) x^{k+1} = prox_{\alpha_k r}(x^k - \alpha_k * ∇g(x^k)) xk+1=proxαkr(xkαkg(xk))

      4. 将原优化问题代入迭代式,可得:

        x k + 1 = argmin u ( λ ∣ ∣ x ∣ ∣ 1 + 1 2 α ∣ ∣ u − x k − 2 α k ∗ ∑ i = 1 10 A i T ( A i x k − b i ) ∣ ∣ 2 2 ) x^{k+1} = \text{argmin}_u \left(\lambda||x||_1 + \frac{1}{2\alpha} ||u - x^k - 2\alpha_k * \sum_{i=1}^{10}A_i^T(A_ix^k-b_i)||_2^2 \right) xk+1=argminu(λ∣∣x1+2α1∣∣uxk2αki=110AiT(Aixkbi)22)

      求解 a r g m i n argmin argmin 用软门限法。由于在计算中无法准确找到次梯度为0的点,我们近似取两步之间的二范数 < 1 0 − 5 <10^{-5} <105 的点。

    2. 数值实验

      def Proximal_Gradient_Method(self):
              A, X, b = self.A, self.X, self.b
      
              alpha = 0.001  # 步长
              P_half = 0.01  # 软门限算子系数
              Xk = np.zeros(self.XSize)  # 初始化变量 Xk
              zero = np.zeros(self.XSize)  # 创建和 Xk 相同形状的零向量
      
              X_dst_from_real = []  # 每步与真实值的距离
              X_dst_steps = []  # 每步的值
              while True:
                  # 计算梯度下降的一步
                  Xk_half = Xk
                  for i in range(self.num):
                      Xk_half -= 2 * alpha * np.dot(A[i].T, np.dot(A[i], Xk) - b[i])
      
                  # 软门限算子
                  Xk_new = zero.copy()  # 创建一个和 Xk 形状相同的零向量,用于存储更新后的 Xk
                  for i in range(self.XSize):
                      if Xk_half[i] < - alpha * P_half:
                          Xk_new[i] = Xk_half[i] + alpha * P_half
                      elif Xk_half[i] > alpha * P_half:
                          Xk_new[i] = Xk_half[i] - alpha * P_half
                      # else 0
      
                  # 记录每步计算结果与真值的距离
                  X_dst_from_real.append(np.linalg.norm(Xk_new - self.X, ord=2))
                  # 记录每一步的值
                  X_dst_steps.append(Xk_new)
      
                  # 判断迭代是否终止
                  if np.linalg.norm(Xk_new - Xk, ord=2) < 10e-5:  # 如果变化小于阈值,停止迭代
                      break
                  else:
                      Xk = Xk_new.copy()  # 更新 Xk 继续迭代
      
              self.__draw(X_dst_from_real, X_dst_steps, "Proximal_Gradient_Method")  # 绘制结果
      
  • 交替方向乘子法

    1. 算法设计

      对于两块可分结构的线性约束凸优化问题:

      m i n x , y   f ( x ) + g ( y ) s . t .   A x + B y = b ,   x ∈ R n , y ∈ R m . min_{x,y}\space f(x) + g(y) \\ s.t.\space Ax + By = b,\space x \in \mathbb{R}^n , y \in \mathbb{R}^m. minx,y f(x)+g(y)s.t. Ax+By=b, xRn,yRm.

      增广拉格朗日函数为:

      L c ( x , y , λ ) = f 1 ( x ) + f 2 ( y ) + λ T ( A x + B y − b ) + c 2 ∣ ∣ A x + B y − b ∣ ∣ 2 2   ,   c > 0. L_c(x, y, λ) = f_1(x) + f_2(y) + λ^T( A x + B y - b) + {c\over 2} ||A x + B y-b||_2^2\space ,\space c>0. Lc(x,y,λ)=f1(x)+f2(y)+λT(Ax+Byb)+2c∣∣Ax+Byb22 , c>0.

      我们对 x x x y y y 依次作极小化:

      x k + 1 = arg min x L c ( x , y k , λ k ) , y k + 1 = arg min y L c ( x k + 1 , y , λ k ) , λ k + 1 = λ k + τ c ( A x k + 1 + B y k + 1 − b ) x^{k+1} = \text{arg min}_x L_c(x, y^k, \lambda^k),\\ y^{k+1} = \text{arg min}_y L_c(x^{k+1}, y, \lambda^k),\\ \lambda^{k+1} = \lambda_k + \tau c(Ax^{k+1} + By^{k+1} - b) xk+1=arg minxLc(x,yk,λk),yk+1=arg minyLc(xk+1,y,λk),λk+1=λk+τc(Axk+1+Byk+1b)

      其中 c > 0 c > 0 c>0为罚因子, τ ∈ ( 0 , 1 + 2 5 ) \tau \in (0, 1+2\sqrt{5}) τ(0,1+25 )为对偶步长。

      对于原问题,我们引入 y = x y=x y=x,一致性约束为 x − y = 0 x − y = 0 xy=0,原问题转换为:

      m i n x ( ∑ i = 1 10 ∥ A i x − b i ∥ 2 2 + p ∥ x ∥ 1 ) s . t .   x − y = 0 min_x \left(\sum_{i=1}^{10} \|A_ix - b_i\|_2^2 +p \|x\|_1 \right)\\ s.t.\space x-y=0 minx(i=110Aixbi22+px1)s.t. xy=0

      其拉格朗日增广函数为:

      L c ( x , y , λ ) = ∑ i = 1 10 ∣ ∣ A i x − b i ∣ ∣ 2 2 + p ∣ ∣ y ∣ ∣ 1 + λ T ( y − x ) + c 2 ∣ ∣ y − x ∣ ∣ 2 2   ,   c > 0. L_c(x, y, λ) = \sum_{i=1}^{10}||A_ix-b_i||_2^2 + p ||y||_1 + λ^T( y-x) + {c\over 2} ||y-x||_2^2\space ,\space c>0. Lc(x,y,λ)=i=110∣∣Aixbi22+p∣∣y1+λT(yx)+2c∣∣yx22 , c>0.

      于是可以得到迭代公式:

      x k + 1 = arg min x ∑ i = 1 10 ∣ ∣ A i x − b i ∣ ∣ 2 2 + p ∣ ∣ y k ∣ ∣ 1 + λ k T ( y k − x ) + c 2 ∣ ∣ y k − x ∣ ∣ 2 2   , y k + 1 = arg min y ∑ i = 1 10 ∣ ∣ A i x k + 1 − b i ∣ ∣ 2 2 + p ∣ ∣ y ∣ ∣ 1 + λ k T ( y − x k + 1 ) + c 2 ∣ ∣ y − x k + 1 ∣ ∣ 2 2   , λ k + 1 = λ k + τ c ( x k + 1 − y k + 1 ) x^{k+1} = \text{arg min}_x\sum_{i=1}^{10}||A_ix-b_i||_2^2 + p ||y^k||_1 + λ{^k}^T( y^k-x) + {c\over 2} ||y^k-x||_2^2\space ,\\ y^{k+1} = \text{arg min}_y\sum_{i=1}^{10}||A_ix^{k+1}-b_i||_2^2 + p ||y||_1 + λ{^k}^T( y-x^{k+1}) + {c\over 2} ||y-x^{k+1}||_2^2\space ,\\ \lambda^{k+1} = \lambda^k + \tau c(x^{k+1} - y^{k+1} ) xk+1=arg minxi=110∣∣Aixbi22+p∣∣yk1+λkT(ykx)+2c∣∣ykx22 ,yk+1=arg minyi=110∣∣Aixk+1bi22+p∣∣y1+λkT(yxk+1)+2c∣∣yxk+122 ,λk+1=λk+τc(xk+1yk+1)

      x x x 部分光滑可导,可以得到 x x x 的梯度; y y y 部分用软门限法来求解其邻近点投影。

      x k + 1 = ( ∑ i = 1 10 A i T A i + c I ) − 1 ( ∑ i = 1 10 A i T b i + c y k − λ k ) , y k + 1 = argmin y   p ∥ y ∥ 1 + c 2 ∥ y − x k + 1 − λ k c ∥ 2 2 x^{k+1} = (\sum_{i=1}^{10} A_i^T A_i + cI)^{-1} (\sum_{i=1}^{10} A_i^T b_i + cy^k - \lambda^k),\\ y^{k+1} = \text{argmin}_y \ p\|y\|_1 + \frac{c}{2} \|y - x^{k+1} - \frac{\lambda^k}{c} \|^2_2 xk+1=(i=110AiTAi+cI)1(i=110AiTbi+cykλk),yk+1=argminy py1+2cyxk+1cλk22

      数值技巧:在更新 x x x 时,可先缓存矩阵 A T A + c I A^TA + cI ATA+cI 的分解以及 A T b A^Tb ATb,以减小后续迭代的计算量。

    2. 数值实验

      def Alternating_Direction_Method_of_Multipliers(self):
              A, X, b = self.A, self.X, self.b
      
              P_half = 0.01
              c = 0.05
              Xk = np.zeros(self.XSize)
              Yk = np.zeros(self.XSize)
              Vk = np.zeros(self.XSize)
      
              X_dst_from_real = []  # 存储每步与真实值的距离
              X_dst_steps = []  # 存储每步的值
      
              while True:
                  # 计算 Xk_new
                  Xk1 = 0
                  Xk2 = 0
                  for i in range(self.num):
                      Xk1 += np.dot(A[i].T, A[i])
                      Xk2 += np.dot(A[i].T, b[i])
                  Xk_new = np.dot(
                      np.linalg.inv(Xk1 + c * np.eye(self.XSize, self.XSize)),
                      c * Yk + Vk + Xk2
                  )
      
                  # 软门限算子更新 Yk_new
                  Yk_new = np.zeros(self.XSize)
                  for i in range(self.XSize):
                      if Xk_new[i] - Vk[i] / c < - P_half / c:
                          Yk_new[i] = Xk_new[i] - Vk[i] / c + P_half / c
                      elif Xk_new[i] - Vk[i] / c > P_half / c:
                          Yk_new[i] = Xk_new[i] - Vk[i] / c - P_half / c
      
                  # 更新 Vk_new
                  Vk_new = Vk + c * (Yk_new - Xk_new)
      
                  # 计算距离并存储步骤
                  X_dst_from_real.append(np.linalg.norm(Xk_new - X, ord=2))
                  X_dst_steps.append(Xk_new)
      
                  # 判断是否收敛
                  if np.linalg.norm(Xk_new - Xk, ord=2) < 10e-5:
                      break
                  else:
                      Xk = Xk_new.copy()
                      Yk = Yk_new.copy()
                      Vk = Vk_new.copy()
      
              self.__draw(X_dst_from_real, X_dst_steps, "Alternating_Direction_Method_of_Multipliers")
      
  • 次梯度法

    1. 算法设计

      1. 选择下降方向:在 x k x^k xk处选择一个负的次梯度 ν k ∈ ∂ f 0 ( x k ) \nu_k \in \partial f_0(x_k) νkf0(xk) 作为下降方向: d k = − ν k d_k = -\nu_k dk=νk ν k ∈ ∂ f 0 ( x k ) \nu_k \in \partial f_0(x_k) νkf0(xk)
      2. 更新下一个点:使用步长 α k \alpha_k αk,迭代更新 x k + 1 = x k + α k d k x_{k+1} = x_k + \alpha_k d_k xk+1=xk+αkdk

      我们的不可微函数的次梯度通过软门限法求得,可微函数的次梯度即为其梯度。

    2. 数值实验

      def Subgradient_Method(self):
              A, X, b = self.A, self.X, self.b
      
              alpha = 0.001
              p_half = 0.01
      
              Xk = np.zeros(self.XSize)
              X_dst_from_real = []  # 每步与真实值的距离
              X_dst_steps = []  # 每步的值
      
              while True:
                  gradient = 0
                  # 求可导部分的梯度
                  for i in range(self.num):
                      gradient += 2 * np.dot(A[i].T, (np.dot(A[i], Xk) - b[i]))
                  # 软门限法 求一范数规范化的次梯度
                  Xtemp = Xk.copy()
                  for i, data in enumerate(Xk):
                      if data == 0:
                          Xtemp[i] = 2 * np.random.random() - 1
                      else:
                          Xtemp[i] = np.sign(Xk[i])
                  subgradient = p_half * Xtemp
                  gradient += subgradient
                  # 迭代更新
                  Xk_new = Xk - alpha * gradient
                  X_dst_from_real.append(np.linalg.norm(Xk_new - X, ord=2))
                  X_dst_steps.append(Xk_new)
                  if np.linalg.norm(Xk_new - Xk, ord=2) < 10e-5:
                      break
                  else:
                      Xk = Xk_new.copy()
      
              self.__draw(X_dst_from_real, X_dst_steps, "Subgradient_Method")
      

3 结果分析

我通过画图来直观地显示结果,结果如下图所示:

  • 蓝色线表示每步计算结果与真值的距离。
  • 橙色线表示每步计算结果与最优解的距离。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

对结果进行分析如下:

  1. 每步计算结果与真值的距离:

    • 邻近点梯度法:该方法在逼近稀疏解时较为有效,收敛速度相对快,在较少的迭代步骤内逼近真实值,但在迭代的后期收敛速度放缓。
    • 交替方向乘子法:一开始迅速收敛,而后逐渐收敛,是三个算法中最逼近真实值的。
    • 次梯度法:收敛速度快,但是与真实值距离是三个算法中最大的。对于稀疏性强的问题,可能需要更多的迭代次数。
  2. 每步计算结果与最优解的距离:

    • 每步计算结果与最优解的距离不断减小,但由于我们的稀疏性较强(200维5稀疏度),邻近点梯度法和次梯度法并不完全接近真实值,而交替方向乘子法表现更好。
    • 对于每种方法,可能需要更多的迭代步骤、较好的正则化 λ λ λ 参数和更优的步长更新算法(我的算法设计里都使用了固定步长,因此某些算法效果可能不是很好)才能得到更接近真实值的结果。
  3. 尝试使用不同的正则化参数 λ λ λ 值,观察算法在不同 λ λ λ下的表现。由算法设计中可以知道, λ λ λ 是软门限法的“门限”,因此,较小的 λ λ λ 会产生更稀疏的解,而较大的 λ λ λ 会得到更多的非零元素,更接近真实值。

完整代码:

import numpy as np
import random

from matplotlib import pyplot as plt


class OptimizedAlgorithm:
    def __init__(self):
        """
        self.num = 10  # 节点数
        self.ASize = (5, 200)  # 测量矩阵大小
        self.XSize = 200  # x向量维度
        self.X # x向量
        self.A = np.array(all_As) # 十个节点的A矩阵数组
        self.b = np.array(all_bs) # 十个节点的b矩阵数组
        """
        # 初始化数据
        self.num = 10  # 节点数
        self.ASize = (5, 200)  # 测量矩阵大小
        self.XSize = 200  # x向量维度
        # 生成x向量
        sparsity = 5
        self.X = np.zeros(self.XSize)
        XIndex = random.sample(list(range(self.XSize)), sparsity)
        for index in XIndex:
            self.X[index] = np.random.randn()  # x向量
        # 存放所有矩阵的列表
        all_As = []
        all_bs = []
        # 生成Ai ei bi
        for _ in range(self.num):
            Ai = np.random.normal(0, 1, self.ASize)
            ei = np.random.normal(0, 0.1, 5)
            bi = np.dot(Ai, self.X) + ei
            all_As.append(Ai)
            all_bs.append(bi)

        self.A = np.array(all_As)  # 十个节点的A矩阵数组
        self.b = np.array(all_bs)  # 十个节点的b矩阵数组

    def __draw(self, X_dst_from_real, X_dst_steps, name):
        X_opt = X_dst_steps[-1]
        for i, data in enumerate(X_dst_steps):
            X_dst_steps[i] = np.linalg.norm(data - X_opt, ord=2)
        plt.title(name + " " + "Distance")
        plt.plot(X_dst_from_real, label='distance-from-real')
        plt.plot(X_dst_steps, label='distance-from-best')
        plt.grid()
        plt.legend()
        plt.show()

    def Proximal_Gradient_Method(self):
        A, X, b = self.A, self.X, self.b

        alpha = 0.001  # 步长
        P_half = 0.01  # 软门限算子系数
        Xk = np.zeros(self.XSize)  # 初始化变量 Xk
        zero = np.zeros(self.XSize)  # 创建和 Xk 相同形状的零向量

        X_dst_from_real = []  # 每步与真实值的距离
        X_dst_steps = []  # 每步的值
        print("Proximal_Gradient_Method is running...")
        while True:
            # 计算梯度下降的一步
            Xk_half = Xk
            for i in range(self.num):
                Xk_half -= 2 * alpha * np.dot(A[i].T, np.dot(A[i], Xk) - b[i])

            # 软门限算子
            Xk_new = zero.copy()  # 创建一个和 Xk 形状相同的零向量,用于存储更新后的 Xk
            for i in range(self.XSize):
                if Xk_half[i] < - alpha * P_half:
                    Xk_new[i] = Xk_half[i] + alpha * P_half
                elif Xk_half[i] > alpha * P_half:
                    Xk_new[i] = Xk_half[i] - alpha * P_half
                # else 0

            # 记录每步计算结果与真值的距离
            X_dst_from_real.append(np.linalg.norm(Xk_new - self.X, ord=2))
            # 记录每一步的值
            X_dst_steps.append(Xk_new)

            # 判断迭代是否终止
            if np.linalg.norm(Xk_new - Xk, ord=2) < 10e-5:  # 如果变化小于阈值,停止迭代
                break
            else:
                Xk = Xk_new.copy()  # 更新 Xk 继续迭代

        self.__draw(X_dst_from_real, X_dst_steps, "Proximal_Gradient_Method")  # 绘制结果

    def Alternating_Direction_Method_of_Multipliers(self):
        A, X, b = self.A, self.X, self.b

        P_half = 0.01
        c = 0.05
        Xk = np.zeros(self.XSize)
        Yk = np.zeros(self.XSize)
        Vk = np.zeros(self.XSize)

        X_dst_from_real = []  # 存储每步与真实值的距离
        X_dst_steps = []  # 存储每步的值
        print("Alternating_Direction_Method_of_Multipliers is running...")

        while True:
            # 计算 Xk_new
            Xk1 = 0
            Xk2 = 0
            for i in range(self.num):
                Xk1 += np.dot(A[i].T, A[i])
                Xk2 += np.dot(A[i].T, b[i])
            Xk_new = np.dot(
                np.linalg.inv(Xk1 + c * np.eye(self.XSize, self.XSize)),
                c * Yk + Vk + Xk2
            )

            # 软门限算子更新 Yk_new
            Yk_new = np.zeros(self.XSize)
            for i in range(self.XSize):
                if Xk_new[i] - Vk[i] / c < - P_half / c:
                    Yk_new[i] = Xk_new[i] - Vk[i] / c + P_half / c
                elif Xk_new[i] - Vk[i] / c > P_half / c:
                    Yk_new[i] = Xk_new[i] - Vk[i] / c - P_half / c

            # 更新 Vk_new
            Vk_new = Vk + c * (Yk_new - Xk_new)

            # 计算距离并存储步骤
            X_dst_from_real.append(np.linalg.norm(Xk_new - X, ord=2))
            X_dst_steps.append(Xk_new)

            # 判断是否收敛
            if np.linalg.norm(Xk_new - Xk, ord=2) < 10e-5:
                break
            else:
                Xk = Xk_new.copy()
                Yk = Yk_new.copy()
                Vk = Vk_new.copy()

        self.__draw(X_dst_from_real, X_dst_steps, "Alternating_Direction_Method_of_Multipliers")

    def Subgradient_Method(self):
        A, X, b = self.A, self.X, self.b

        alpha = 0.001
        p_half = 0.01

        Xk = np.zeros(self.XSize)
        X_dst_from_real = []  # 每步与真实值的距离
        X_dst_steps = []  # 每步的值
        print("Subgradient_Method is running...")

        while True:
            gradient = 0
            # 求可导部分的梯度
            for i in range(self.num):
                gradient += 2 * np.dot(A[i].T, (np.dot(A[i], Xk) - b[i]))
            # 软门限法 求一范数规范化的次梯度
            Xtemp = Xk.copy()
            for i, data in enumerate(Xk):
                if data == 0:
                    Xtemp[i] = 2 * np.random.random() - 1
                else:
                    Xtemp[i] = np.sign(Xk[i])
            subgradient = p_half * Xtemp
            gradient += subgradient
            # 迭代更新
            Xk_new = Xk - alpha * gradient
            X_dst_from_real.append(np.linalg.norm(Xk_new - X, ord=2))
            X_dst_steps.append(Xk_new)
            if np.linalg.norm(Xk_new - Xk, ord=2) < 10e-5:
                break
            else:
                Xk = Xk_new.copy()

        self.__draw(X_dst_from_real, X_dst_steps, "Subgradient_Method")


if __name__ == '__main__':
    oa = OptimizedAlgorithm()
    oa.Proximal_Gradient_Method()
    oa.Alternating_Direction_Method_of_Multipliers()
    oa.Subgradient_Method()

参考

  • 23
    点赞
  • 26
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值