本文主要是介绍优化|求解非凸和无梯度lipschitz连续性的一阶算法在二次规划反问题中的应用(代码分享),希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!
原文信息(包括题目、发表期刊、原文链接等):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 { Ψ ( x ) : = g ( x ) + θ f ( x ) : x ∈ R d } \min\Big\{ \Psi(x):=g(x) + \theta f(x): x\in \mathbb{R}^{d}\Big\} min{Ψ(x):=g(x)+θf(x):x∈Rd}
其中 g ( x ) = 1 4 ∑ i = 1 m ( x T A i x − b i ) 2 , f ( x ) = ∥ x ∥ 1 g(x) = \frac{1}{4}\sum_{i=1}^{m}(x^{T}A_{i}x - b_{i})^2, f(x) = \|x\|_{1} g(x)=41∑i=1m(xTAix−bi)2,f(x)=∥x∥1,而且 A i A_{i} Ai是对称矩阵。
2 求解方法
在给出求解方法之前,我们首先定义
p λ ( x ) = λ ∇ g ( x ) − ∇ h ( x ) p_{\lambda}(x)=\lambda \nabla g(x)-\nabla h(x) pλ(x)=λ∇g(x)−∇h(x)
和软阈值算子
S τ ( y ) = max { ∣ y ∣ − τ , 0 } sgn ( y ) = max ( y − τ , 0 ) − max ( − y − τ , 0 ) ( 5.1 ) S_{\tau}(y)=\max\{|y|-\tau, 0\}\text{sgn}(y)=\max(y-\tau,0) - \max(-y-\tau,0) \qquad (5.1) Sτ(y)=max{∣y∣−τ,0}sgn(y)=max(y−τ,0)−max(−y−τ,0)(5.1)
为保证函数 g ( x ) , f ( x ) g(x),f(x) g(x),f(x)是L-smad,我们令
h ( x ) = 1 4 ∥ x ∥ 2 4 + 1 2 ∥ x ∥ 2 2 , h(x) = \frac{1}{4} \| x \|_2^4 + \frac{1}{2} \| x \|_2^2, h(x)=41∥x∥24+21∥x∥22,
具体见原文引理5.1。
本文的求解方法主要根据原文的命题5.1,如下所示
命题5.1 ( l 1 l_{1} l1范数正则化的Bregman近似公式) 令 f = ∥ ⋅ ∥ 1 f=\|\cdot\|_{1} f=∥⋅∥1且对 x ∈ R d x\in \mathbb{R}^{d} x∈Rd,令 v ( x ) : = S λ θ ( p λ ( x ) ) v(x):=S_{\lambda \theta}(p_{\lambda}(x)) v(x):=Sλθ(pλ(x))。那么,可得 x + = T λ ( x ) x^{+}=T_{\lambda}(x) x+=Tλ(x)为
x + = − t ∗ v ( x ) = t ∗ S λ θ ( ∇ h ( x ) − λ ∇ g ( x ) ) ( 5.2 ) x^{+}=-t^{*}v(x)=t^{*}S_{\lambda\theta}(\nabla h(x)-\lambda\nabla g(x)) \qquad (5.2) x+=−t∗v(x)=t∗Sλθ(∇h(x)−λ∇g(x))(5.2)
是显示公式,其中 t ∗ t^{*} t∗是下面方程的唯一正实根,
t 3 ∥ v ( x ) ∥ 2 2 + t − 1 = 0. ( 5.3 ) t^{3}\|v(x)\|_{2}^{2}+t-1=0. \qquad (5.3) t3∥v(x)∥22+t−1=0.(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 ∗ t^{*} 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) 计算 g ( x ) = 1 4 ∑ i = 1 m ( x T A i x − b i ) 2 g(x) = \frac{1}{4}\sum_{i=1}^{m}(x^{T}A_{i}x - b_{i})^2 g(x)=41∑i=1m(xTAix−bi)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);endreturn der;
end
(5) 计算 h ( x ) = 1 4 ∥ x ∥ 2 4 + 1 2 ∥ x ∥ 2 2 h(x)=\frac{1}{4}\|x\|_{2}^{4}+\frac{1}{2}\|x\|_{2}^{2} h(x)=41∥x∥24+21∥x∥22的导数
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
endfor k in range(1,m)b[k] = rand(d,1)[1];
end
(8) 根据引理5.1的结果可知 L ≥ ∑ i = 1 m 3 ∥ A i ∥ 2 + ∥ A i ∥ ∣ b i ∣ L\geq \sum_{i=1}^{m}3\|A_{i}\|^{2}+\|A_{i}\||b_{i}| L≥∑i=1m3∥Ai∥2+∥Ai∥∣bi∣。另外,根据定理 4.1 成立的条件 0 < λ L < 1 0<\lambda L<1 0<λL<1,可得 0 < λ < 1 L 0<\lambda<\frac{1}{L} 0<λ<L1。
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")
这篇关于优化|求解非凸和无梯度lipschitz连续性的一阶算法在二次规划反问题中的应用(代码分享)的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!