优化｜求解非凸和无梯度lipschitz连续性的一阶算法在二次规划反问题中的应用(代码分享)

news2024/11/27 22:31:03

在这里插入图片描述

原文信息（包括题目、发表期刊、原文链接等）：First Order Methods Beyond Convexity and Lipschitz Gradient Continuity with Applications to Quadratic Inverse Problems
原文作者：Jérôme Bolte, Shoham Sabach, Marc Teboulle, and Yakov Vaisbourd
代码分享者：李朋

1 问题描述

考虑下面的二次规划反问题
$\min\Big\{ \Psi(x):=g(x) + \theta f(x): x\in \mathbb{R}^{d}\Big\}$

其中 $\frac{1}{4}\sum_{i=1}^{m}(x^{T}A_{i}x - b_{i})^2, f(x) = \|x\|_{1}$ ，而且 $A_{i}$ 是对称矩阵。

2 求解方法

在给出求解方法之前，我们首先定义
$p_{\lambda}(x)=\lambda \nabla g(x)-\nabla h(x)$
和软阈值算子
$S_{\tau}(y)=\max\{|y|-\tau, 0\}\text{sgn}(y)=\max(y-\tau,0) - \max(-y-\tau,0) \qquad (5.1)$
为保证函数 $g (x), f (x)$ 是L-smad，我们令
$\frac{1}{4} \| x \|_2^4 + \frac{1}{2} \| x \|_2^2,$
具体见原文引理5.1。

本文的求解方法主要根据原文的命题5.1，如下所示

命题5.1 ( $l_{1}$ 范数正则化的Bregman近似公式) 令 $f=\|\cdot\|_{1}$ 且对 $x\in \mathbb{R}^{d}$ ，令 $v(x):=S_{\lambda \theta}(p_{\lambda}(x))$ 。那么，可得 $x^{+}=T_{\lambda}(x)$ 为
$x^{+}=-t^{*}v(x)=t^{*}S_{\lambda\theta}(\nabla h(x)-\lambda\nabla g(x)) \qquad (5.2)$

是显示公式，其中 $t^{*}$ 是下面方程的唯一正实根，
$t^{3}\|v(x)\|_{2}^{2}+t-1=0. \qquad (5.3)$

3 代码实现

在本次仿真中，我们采用Julia语言编写一个求解二次规划反问题的算法 (5-2)。

(1) 用using 添加一些要用到的库。

using Roots
using LinearAlgebra
using SparseArrays
using Distributions
using Random
using Printf
using Plots
using Polynomials

(2) 根据公式 (5-1) 定义软阈值函数

function compute_softThreshold(y,τ)
    p = max.(y.-τ,0) - max.(-y.-τ,0);
    return p;
end

(3)根据公式(5-3) 计算 $t^{*}$

function find_positiveRoot(S)
    t = variable();
    v = sum(S.^2);
    
    f = t^3*v + t -1;
    t_opt = find_zero(f,(0,1));
    
    return t_opt;
end

(4) 计算 $\frac{1}{4}\sum_{i=1}^{m}(x^{T}A_{i}x - b_{i})^2$ 的导数

function derivative_g(A,b,x,m,n)
    # compute the derivative of g(x)
    der = zeros(n,1);
    for k in range(1,m)
        der = der + (transpose(x)*A[k]*x.-b[k]).*(A[k]*x);
    end
    return der;
end

(5) 计算 $h(x)=\frac{1}{4}\|x\|_{2}^{4}+\frac{1}{2}\|x\|_{2}^{2}$ 的导数

function derivative_h(x)
    # compute the derivative of h(x)
    der = (sum(x.^2) + 1).*x;
    return der;
end

(6) 全局参数

# Global Parameters
MAXITE = 500;
m =3;
n = 2;

(7) 生成问题数据

θ = 0.5;

Random.seed!(123);

A = Array{Matrix}(undef,m);
b = Array{Float64}(undef,m); 


d = Normal(2,2);
for k in range(1,m)
    A[k] = rand(d,n,n)
    A[k] = (transpose(A[k])+A[k])./2
end

for k in range(1,m)
    b[k] = rand(d,1)[1];
end

(8) 根据引理5.1的结果可知 $L\geq \sum_{i=1}^{m}3\|A_{i}\|^{2}+\|A_{i}\||b_{i}|$ 。另外，根据定理 4.1 成立的条件 $0<\lambda L<1$ ，可得 $0<\lambda<\frac{1}{L}$ 。

L = sum([3*norm(A[k]).^2 + norm(A[k])*norm(b[k]) for k =1:m])+1;
λ = 1/L;   #λ≤1/L

(9) 主程序

x = ones(n,1)
objval_vec = zeros(1,MAXITE);  #存储计算过程中目标函数值
x_vec = zeros(n,MAXITE);       #存储计算过程中变量值

for k in range(1,MAXITE)
    
    #计算、存储当前目标函数值
    objval = sum([1/4*(transpose(x)*A[k]*x.-b[k])^2 for k=1:m]) .+ θ.*norm(x,1); 
    objval_vec[1,k] = objval[1,1];  
    
    #存储当前变量值
    x_vec[:,k] = x; 
    
    #计算函数g(x)、h(x)当前时刻的导数值
    xold = x;
    der_h = derivative_h(xold);
    der_g = derivative_g(A,b,xold,m,n);
    
    y = λ*der_g - der_h;
    τ = λ * θ;
    v = compute_softThreshold(y,τ);   #计算公式(5-2)中的软阈值算子部分   
    
    topt = find_positiveRoot(v);  #计算公式(5-2)中的 t*
    
    x = -topt.*v; # 根据公式(5-2) 求出下一时刻 x 的值
end
print("最优解：",x,"\n");
print("最小目标值：",objval_vec[end]);

(10) 画出目标函数值随计算步数的变化

K = range(1, MAXITE);
plot(K, [objval_vec[k] for k=1:MAXITE], yaxis=:log10,label="object value")

(11) 画出变量值随计算步数的变化

plot(x_vec[1,1:MAXITE], x_vec[2,1:MAXITE], arrow = :arrow)
scatter!([x_vec[1,1]], [x_vec[2,1]], markshape=:rect, marksize = 5, markercolor= :red, legend = false)
xlabel!("x1")
ylabel!("x2")