Optimization problem. The optimization problem is
$$ \begin{array}{ll} \textrm{minimize} & f(x)\\ \textrm{subject to} & x\in\mathbf{R}^{d}, \end{array} $$where $f\in\mathcal{F}_{0,L}$ is the class of $L$-smooth convex functions.
Fixed-step first-order algorithm. The optimization algorithm in consideration is:
$$ \begin{align*} \left(\forall_{i\in[0:N]}\right)\;w_{i} & =w_{i-1}-\sum_{j\in[0:i-1]}\frac{h_{i,j}}{L}\nabla f(w_{j})\\ & =w_{0}-\sum_{i\in[0:i-1]}\frac{\alpha_{i,j}}{L}\nabla f(w_{j}), \end{align*} $$where we use the notation that $\alpha_{0,i}=h_{0,i}=0$ and $\{\alpha_{i,j}\}$ and $\{h_{i,j}\}$ are related via the following relationship:
$$ \left(\forall i\in[1:N]\right)\;\left(\forall j\in[0:i-1]\right)\quad h_{i,j}=\begin{cases} \alpha_{i,j}, & \textrm{if }j=i-1\\ \alpha_{i,j}-\alpha_{i-1,j}, & \textrm{if }j\in[0:i-1]. \end{cases} $$Notation. Denote,
$$
\begin{align*}
J & =\{(i,j)\mid j=i+1,i\in[0:N-1]\}\cup\{(i,j)\mid i=\star,j\in[0:N]\},
\end{align*}
$$
and
$$
\mathbf{w}_{0}=e_{1}\in\mathbf{R}^{N+2},\mathbf{g}_{i}=e_{i+2}\in\mathbf{R}^{N+2},\mathbf{f}_{i}=e_{i+1}\in\mathbf{R}^{N+1},
$$
where $e_{i}$ is the unit vector with $i$th component equal to $1$ and the rest being zero. Next, define
$$
\begin{align*}
& S\left(\tau,\{\lambda_{i,j}\},\{\alpha_{i,j}^{\prime}\}\right)\\
= & c_{w}\tau\mathbf{w}_{0}\mathbf{w}_{0}^{\top}+\\
& \frac{1}{2L}\Bigg[\sum_{i\in[0:N-1]}\lambda_{i,i+1}\left\{ (\mathbf{g}_{i}-\mathbf{g}_{i+1})\odot(\mathbf{g}_{i}-\mathbf{g}_{i+1})\right\} +\\
& \sum_{j\in[0:N]}\lambda_{\star,j}\left\{ (\mathbf{g}_{\star}-\mathbf{g}_{j})\odot(\mathbf{g}_{\star}-\mathbf{g}_{j})\right\} \Bigg]\\
& -\lambda_{\star,0}\{\mathbf{g}_{0}\odot\mathbf{w}_{0}\}\\
& -\sum_{i\in[1:N-1]}\Bigg[\lambda_{i,i+1}\left\{ \mathbf{g}_{i}\odot\mathbf{w}_{0}\right\} -\sum_{j\in[0:i-1]}\frac{\alpha_{i,j}^{\prime}}{L}\left\{ \mathbf{g}_{i}\odot\mathbf{g}_{j}\right\} \Bigg]\\
& -\left((\mathbf{g}_{N}\odot\mathbf{w}_{0})-\sum_{j\in[0:N-1]}\frac{\alpha_{N,j}^{\prime}}{L}\left\{ \mathbf{g}_{N}\odot\mathbf{g}_{j}\right\} \right)\\
& +\sum_{i\in[0:N-1]}\Bigg[\lambda_{i,i+1}\left\{ \mathbf{g}_{i+1}\odot\mathbf{w}_{0}\right\} -\sum_{j\in[0:i-1]}\frac{\alpha_{i,j}^{\prime}}{L}\left\{ \mathbf{g}_{i+1}\odot\mathbf{g}_{j}\right\} \Bigg],
\end{align*}
$$
where
$$
\alpha_{i,j}^{\prime}=\begin{cases}
\lambda_{i,i+1}\alpha_{i,j}, & \textrm{if }i\in[1:N-1],j\in[0:i-1]\\
\alpha_{N,j}, & \textrm{if }i\in N,j\in[0:N-1]\\
\lambda_{0,1}\underbrace{\alpha_{0,j}}_{=0}=0, & \textrm{if }i=0.
\end{cases}
$$
Then, the performance estimation optimization algorithm is:
$$
\begin{array}{ll}
\textrm{minimize} & \tau\\
\textrm{subject to} & -\mathbf{f}_{N}+\mathbf{f}_{\star}+\sum_{(i,j)\in J}\lambda_{i,j}\left(\mathbf{f}_{j}-\mathbf{f}_{i}\right)+c_{f}\tau\left(\mathbf{f}_{0}-\mathbf{f}_{\star}\right)=0\\
& S\left(\tau,\{\lambda_{i,j}\},\{\alpha_{i,j}^{\prime}\}\right)\succeq0\\
& \left(\forall(i,j)\in J\right)\quad\lambda_{i,j}\geq0\\
& \tau\geq0,
\end{array}
$$
where the decision variables are $\tau,\{\lambda_{i,j}\},$ and $\{\alpha_{i,j}^{\prime}\}.$ Let us see, how we can solve this problem in Julia
.
# Load the packages
using JuMP, MosekTools, Mosek, LinearAlgebra, SCS, COSMO, Literate, OffsetArrays
# %% Some helper functions
# %% construct e_i in R^n
function e_i(n, i)
e_i_vec = zeros(n, 1)
e_i_vec[i] = 1
return e_i_vec
end
# %% construct symmetric outer product
function ⊙(a,b)
return ((a*b')+(b*a')) ./ 2
end
⊙ (generic function with 1 method)
# %% Parameters to be tuned
N = 5
L = 1
μ = 0
c_w = 1
c_f = 0
0
Next, define the bold vectors:
$$ \mathbf{w}_{0}=e_{1}\in\mathbf{R}^{N+2},\mathbf{g}_{i}=e_{i+2}\in\mathbf{R}^{N+2},\mathbf{f}_{i}=e_{i+1}\in\mathbf{R}^{N+1}, $$where $e_{i}$ is the unit vector with $i$th component equal to $1$ and the rest being zero.
# define all the bold vectors
# --------------------------
# define 𝐰_0
𝐰_0 = e_i(N+2, 1)
𝐰_star = zeros(N+2, 1)
𝐠_star = zeros(N+2,1)
# define 𝐠_0, 𝐠_1, …, 𝐠_N
# first we define 𝐠_Julia vectors and then 𝐠 vectors
𝐠 = OffsetArray(zeros(N+2, N+1), 1:N+2, 0:N)
# 𝐠= [𝐠_0 𝐠_1 𝐠_2 ... 𝐠_N]
for i in 0:N
𝐠[:,i] = e_i(N+2, i+2)
end
𝐟_star = zeros(N+1,1)
# time to define 𝐟_0, 𝐟_1, …, 𝐟_N
𝐟 = OffsetArray(zeros(N+1, N+1), 1:N+1, 0:N)
# 𝐟 = [𝐟_0, 𝐟_1, …, 𝐟_N]
for i in 0:N
𝐟[:,i] = e_i(N+1, i+1)
end
Next, we define our decision variables using JuMP
.
pep_model = Model(optimizer_with_attributes(Mosek.Optimizer, "MSK_DPAR_INTPNT_CO_TOL_PFEAS" => 1e-10))
# time to define all the variables
# -------------------------------
# define α′ (can be typed as \alpha[TAB]\prime[TAB])
@variable(pep_model, α′[1:N, 0:N-1])
# define λ variables
@variable(pep_model, λ_i_ip1[0:N-1] >= 0)
# this defines (λ_{i,i+1})_{i∈[0:N-1]} in Julia indexing
@variable(pep_model, λ_star_i[0:N] >= 0)
# this defines (λ_{⋆,i})_{i∈[0:N]} in Julia indexing
# define τ
@variable(pep_model, τ >= 0)
Define the objective next, which is to minimize $\tau$.
# define objective
@objective(
pep_model,
Min,
τ
)
Time to add the linear constraint $$ \begin{align*} & -\mathbf{f}_{N}+\mathbf{f}_{\star}+\sum_{(i,j)\in J}\lambda_{i,j}\left(\mathbf{f}_{j}-\mathbf{f}_{i}\right)+c_{f}\tau\left(\mathbf{f}_{0}-\mathbf{f}_{\star}\right)=0\\ \Leftrightarrow & c_{f}\tau\left(\mathbf{f}_{0}-\mathbf{f}_{\star}\right)+\sum_{i\in[0:N-1]}\lambda_{i,i+1}(\mathbf{f}_{i+1}-\mathbf{f}_{i})+\sum_{i\in[0:N]}\lambda_{\star,i}(\mathbf{f}_{i}-\mathbf{f}_{\star})=\mathbf{f}_{N}-\mathbf{f}_{\star} \end{align*} $$ first, which are much simpler.
# Add the linear equality constraint
# ----------------------------------
@constraint(pep_model,
τ * c_f .* 𝐟[:,0] + sum(λ_i_ip1[i] .* (𝐟[:,i+1]-𝐟[:,i]) for i in 0:N-1)
+ sum(λ_star_i[i] .* (𝐟[:,i] - 𝐟_star) for i in 0:N)
.== 𝐟[:,N] - 𝐟_star
)
6×1 OffsetArray(::Matrix{ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape}}, 1:6, 1:1) with eltype ConstraintRef{Model, MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64}, MathOptInterface.EqualTo{Float64}}, ScalarShape} with indices 1:6×1:1: -λ_i_ip1[0] + λ_star_i[0] == 0.0 λ_i_ip1[0] - λ_i_ip1[1] + λ_star_i[1] == 0.0 λ_i_ip1[1] - λ_i_ip1[2] + λ_star_i[2] == 0.0 λ_i_ip1[2] - λ_i_ip1[3] + λ_star_i[3] == 0.0 λ_i_ip1[3] - λ_i_ip1[4] + λ_star_i[4] == 0.0 λ_i_ip1[4] + λ_star_i[5] == 1.0
Now let us construct the giant sdp constraint step by step. It has 6 summands, which we add one by one.
# Add the giant LMI constraint step by step
# ----------------------------------------
Term 1 is $\texttt{term}_1 = c_{w}\tau\mathbf{w}_{0}\mathbf{w}_{0}^{\top}.$
term_1 = c_w * τ * ⊙(𝐰_0,𝐰_0)
7×7 Matrix{AffExpr}: τ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Term 2 is $$ \texttt{term}_2 = \frac{1}{2L}\Bigg[\sum_{i\in[0:N-1]}\lambda_{i,i+1}\left\{ (\mathbf{g}_{i}-\mathbf{g}_{i+1})\odot(\mathbf{g}_{i}-\mathbf{g}_{i+1})\right\} +\\\sum_{j\in[0:N]}\lambda_{\star,j}\left\{ (\mathbf{g}_{\star}-\mathbf{g}_{j})\odot(\mathbf{g}_{\star}-\mathbf{g}_{j})\right\} \Bigg]. $$
term_2_part_1 = sum(λ_i_ip1[i] .* ⊙(𝐠[:,i]-𝐠[:,i+1],𝐠[:,i]-𝐠[:,i+1]) for i in 0:N-1)
7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7: 0 0 0 … 0 0 λ_i_ip1[0] -λ_i_ip1[0] 0 0 -λ_i_ip1[0] λ_i_ip1[0] + λ_i_ip1[1] 0 0 0 -λ_i_ip1[1] 0 0 0 0 0 0 0 0 … -λ_i_ip1[4] 0 0 0 λ_i_ip1[4]
term_2_part_2 = sum(λ_star_i[i] .* ⊙(𝐠_star - 𝐠[:,i],𝐠_star - 𝐠[:,i]) for i in 0:N)
7×7 Matrix{AffExpr}: 0 0 0 0 … 0 0 0 λ_star_i[0] 0 0 0 0 0 0 λ_star_i[1] 0 0 0 0 0 0 λ_star_i[2] 0 0 0 0 0 0 0 0 0 0 0 0 … λ_star_i[4] 0 0 0 0 0 0 λ_star_i[5]
term_2 = (term_2_part_1 + term_2_part_2) ./ (2*L)
7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7: 0 0 … 0 0 0.5 λ_i_ip1[0] + 0.5 λ_star_i[0] 0 0 -0.5 λ_i_ip1[0] 0 0 0 0 0 0 0 0 0 … -0.5 λ_i_ip1[4] 0 0 0.5 λ_i_ip1[4] + 0.5 λ_star_i[5]
Term 3 is $\texttt{term}_3 = -\lambda_{\star,0}\{\mathbf{g}_{0}\odot\mathbf{w}_{0}\}$.
term_3 = - λ_star_i[0] .* ⊙(𝐠[:,0],𝐰_0)
7×7 Matrix{AffExpr}: 0 -0.5 λ_star_i[0] 0 0 0 0 0 -0.5 λ_star_i[0] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Term 4 is $$ \texttt{term}_4 = -\sum_{i\in[1:N-1]}\Bigg[\lambda_{i,i+1}\left\{ \mathbf{g}_{i}\odot\mathbf{w}_{0}\right\} -\sum_{j\in[0:i-1]}\frac{\alpha_{i,j}^{\prime}}{L}\left\{ \mathbf{g}_{i}\odot\mathbf{g}_{j}\right\} \Bigg]. $$
term_4_part_1 = - sum( λ_i_ip1[i] .* ⊙(𝐠[:,i],𝐰_0) for i in 1:N-1)
7×7 Matrix{AffExpr}: 0 0 -0.5 λ_i_ip1[1] … -0.5 λ_i_ip1[3] -0.5 λ_i_ip1[4] 0 0 0 0 0 0 0 -0.5 λ_i_ip1[1] 0 0 0 0 0 -0.5 λ_i_ip1[2] 0 0 0 0 0 -0.5 λ_i_ip1[3] 0 0 0 0 0 -0.5 λ_i_ip1[4] 0 0 … 0 0 0 0 0 0 0 0 0
term_4_part_2 = (1/L) .*
sum( sum(α′[i,j] .* ⊙(𝐠[:,i],𝐠[:,j]) for j in 0:i-1) for i in 1:N-1)
7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7: 0 0 0 0 0 0 0 0 0 0.5 α′[1,0] 0.5 α′[2,0] 0.5 α′[3,0] 0.5 α′[4,0] 0 0 0.5 α′[1,0] 0 0.5 α′[2,1] 0.5 α′[3,1] 0.5 α′[4,1] 0 0 0.5 α′[2,0] 0.5 α′[2,1] 0 0.5 α′[3,2] 0.5 α′[4,2] 0 0 0.5 α′[3,0] 0.5 α′[3,1] 0.5 α′[3,2] 0 0.5 α′[4,3] 0 0 0.5 α′[4,0] 0.5 α′[4,1] 0.5 α′[4,2] 0.5 α′[4,3] 0 0 0 0 0 0 0 0 0
term_4 = term_4_part_1 + term_4_part_2
7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7: 0 0 … -0.5 λ_i_ip1[3] -0.5 λ_i_ip1[4] 0 0 0 0.5 α′[3,0] 0.5 α′[4,0] 0 -0.5 λ_i_ip1[1] 0.5 α′[1,0] 0.5 α′[3,1] 0.5 α′[4,1] 0 -0.5 λ_i_ip1[2] 0.5 α′[2,0] 0.5 α′[3,2] 0.5 α′[4,2] 0 -0.5 λ_i_ip1[3] 0.5 α′[3,0] 0 0.5 α′[4,3] 0 -0.5 λ_i_ip1[4] 0.5 α′[4,0] … 0.5 α′[4,3] 0 0 0 0 0 0 0
Term 5 is $$ \texttt{term}_5 = -\left((\mathbf{g}_{N}\odot\mathbf{w}_{0})-\sum_{j\in[0:N-1]}\frac{\alpha_{N,j}^{\prime}}{L}\left\{ \mathbf{g}_{N}\odot\mathbf{g}_{j}\right\} \right). $$
term_5 = - ( ⊙(𝐠[:,N],𝐰_0) - (1/L) .* sum(α′[N,j] .* ⊙(𝐠[:,N],𝐠[:,j]) for j in 0:N-1))
7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7: 0 0 0 0 … 0 -0.5 0 0 0 0 0 0.5 α′[5,0] 0 0 0 0 0 0.5 α′[5,1] 0 0 0 0 0 0.5 α′[5,2] 0 0 0 0 0 0.5 α′[5,3] 0 0 0 0 … 0 0.5 α′[5,4] -0.5 0.5 α′[5,0] 0.5 α′[5,1] 0.5 α′[5,2] 0.5 α′[5,4] 0
Okay, we are almost there. One more term. Term 6 is: $$ \texttt{term}_6 = \sum_{i\in[0:N-1]}\Bigg[\lambda_{i,i+1}\left\{ \mathbf{g}_{i+1}\odot\mathbf{w}_{0}\right\} -\sum_{j\in[0:i-1]}\frac{\alpha_{i,j}^{\prime}}{L}\left\{ \mathbf{g}_{i+1}\odot\mathbf{g}_{j}\right\} \Bigg]. $$
term_6_part_1 = sum(λ_i_ip1[i] .* ⊙(𝐠[:,i+1],𝐰_0) for i in 0:N-1)
7×7 Matrix{AffExpr}: 0 0 0.5 λ_i_ip1[0] … 0.5 λ_i_ip1[3] 0.5 λ_i_ip1[4] 0 0 0 0 0 0.5 λ_i_ip1[0] 0 0 0 0 0.5 λ_i_ip1[1] 0 0 0 0 0.5 λ_i_ip1[2] 0 0 0 0 0.5 λ_i_ip1[3] 0 0 … 0 0 0.5 λ_i_ip1[4] 0 0 0 0
term_6_part_2 = - (1/L) .*
sum( sum( α′[i,j] .* ⊙(𝐠[:,i+1],𝐠[:,j]) for j in 0:i-1) for i in 1:N-1)
7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7: 0 0 0 0 … 0 0 0 0 0 -0.5 α′[1,0] -0.5 α′[3,0] -0.5 α′[4,0] 0 0 0 0 -0.5 α′[3,1] -0.5 α′[4,1] 0 -0.5 α′[1,0] 0 0 -0.5 α′[3,2] -0.5 α′[4,2] 0 -0.5 α′[2,0] -0.5 α′[2,1] 0 0 -0.5 α′[4,3] 0 -0.5 α′[3,0] -0.5 α′[3,1] -0.5 α′[3,2] … 0 0 0 -0.5 α′[4,0] -0.5 α′[4,1] -0.5 α′[4,2] 0 0
term_6 = term_6_part_1 + term_6_part_2
7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7: 0 0 … 0.5 λ_i_ip1[3] 0.5 λ_i_ip1[4] 0 0 -0.5 α′[3,0] -0.5 α′[4,0] 0.5 λ_i_ip1[0] 0 -0.5 α′[3,1] -0.5 α′[4,1] 0.5 λ_i_ip1[1] -0.5 α′[1,0] -0.5 α′[3,2] -0.5 α′[4,2] 0.5 λ_i_ip1[2] -0.5 α′[2,0] 0 -0.5 α′[4,3] 0.5 λ_i_ip1[3] -0.5 α′[3,0] … 0 0 0.5 λ_i_ip1[4] -0.5 α′[4,0] 0 0
Time to add all these terms to construct $S\left(\tau,\{\lambda_{i,j}\},\{\alpha_{i,j}^{\prime}\}\right)$: $$ S\left(\tau,\{\lambda_{i,j}\},\{\alpha_{i,j}^{\prime}\}\right) = \sum_{i\in [1:6]}\texttt{term}_i $$
# oof, okay constructed the terms for the LMI constraint, let us hope that there is no mistake and add them together
S_mat = term_1 + term_2 + term_3 + term_4 + term_5 + term_6
7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7: τ … 0.5 λ_i_ip1[4] - 0.5 -0.5 λ_star_i[0] 0.5 α′[5,0] - 0.5 α′[4,0] -0.5 λ_i_ip1[1] + 0.5 λ_i_ip1[0] 0.5 α′[5,1] - 0.5 α′[4,1] -0.5 λ_i_ip1[2] + 0.5 λ_i_ip1[1] 0.5 α′[5,2] - 0.5 α′[4,2] -0.5 λ_i_ip1[3] + 0.5 λ_i_ip1[2] 0.5 α′[5,3] - 0.5 α′[4,3] -0.5 λ_i_ip1[4] + 0.5 λ_i_ip1[3] … -0.5 λ_i_ip1[4] + 0.5 α′[5,4] 0.5 λ_i_ip1[4] - 0.5 0.5 λ_i_ip1[4] + 0.5 λ_star_i[5]
Let's add the LMI cosntraint $S\left(\tau,\{\lambda_{i,j}\},\{\alpha_{i,j}^{\prime}\}\right)\succeq0$.
# add the LMI constraint now
@constraint(pep_model,
S_mat in PSDCone()
)
Time to optimize the code now!
# time to optimize!
# -----------------
optimize!(pep_model)
println("termination status =", termination_status(pep_model) )
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 34 Cones : 0 Scalar variables : 37 Matrix variables : 1 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator started. Freed constraints in eliminator : 5 Eliminator terminated. Eliminator started. Freed constraints in eliminator : 0 Eliminator terminated. Eliminator - tries : 2 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - number : 0 Presolve terminated. Time: 0.00 Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 34 Cones : 0 Scalar variables : 37 Matrix variables : 1 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 29 Optimizer - Cones : 1 Optimizer - Scalar variables : 23 conic : 16 Optimizer - Semi-definite variables: 1 scalarized : 28 Factor - setup time : 0.00 dense det. time : 0.00 Factor - ML order time : 0.00 GP order time : 0.00 Factor - nonzeros before factor : 423 after factor : 435 Factor - dense dim. : 0 flops : 1.21e+04 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 1.0e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.00 1 2.7e-01 2.7e-01 1.2e-01 2.14e+00 -2.378627944e-02 -6.117274954e-03 2.7e-01 0.00 2 1.2e-01 1.2e-01 2.9e-02 1.32e+00 7.175923608e-02 7.139302023e-02 1.2e-01 0.00 3 4.2e-02 4.2e-02 6.9e-03 1.05e+00 4.477758503e-02 4.729110878e-02 4.2e-02 0.00 4 2.3e-02 2.3e-02 2.9e-03 7.37e-01 2.880293993e-02 3.057913998e-02 2.3e-02 0.00 5 6.7e-03 6.7e-03 4.0e-04 9.85e-01 3.356646000e-02 3.302278094e-02 6.7e-03 0.00 6 2.5e-03 2.5e-03 1.1e-04 6.92e-01 2.291866424e-02 2.288414204e-02 2.5e-03 0.00 7 6.6e-04 6.6e-04 1.5e-05 8.71e-01 2.066303721e-02 2.066465445e-02 6.6e-04 0.00 8 9.8e-05 9.8e-05 8.8e-07 8.83e-01 1.883068878e-02 1.882786667e-02 9.8e-05 0.00 9 1.7e-06 1.7e-06 1.9e-09 9.93e-01 1.859217137e-02 1.859212284e-02 1.7e-06 0.00 10 4.0e-08 4.0e-08 7.3e-12 1.00e+00 1.858820828e-02 1.858820710e-02 4.0e-08 0.00 11 2.1e-09 2.1e-09 8.5e-14 1.00e+00 1.858813924e-02 1.858813917e-02 2.1e-09 0.00 12 7.9e-11 9.5e-11 6.3e-16 1.00e+00 1.858813673e-02 1.858813673e-02 7.9e-11 0.00 Optimizer terminated. Time: 0.01 termination status =OPTIMAL
We now collect the optimal values of our decision variables.
# Collect the decision variables
# -----------------------------
α′_opt = OffsetArray(value.(α′), 1:N, 0:N-1)
λ_opt_i_ip1 = OffsetVector(value.(λ_i_ip1), 0:N-1)
λ_opt_star_i = OffsetVector(value.(λ_star_i), 0:N)
τ_opt = value(τ)
0.018588136733035998
Next, we recover $\alpha_{i,j}$ from $\alpha^\prime_{i,j}$ using the formula $$ \alpha_{i,j}=\begin{cases} \frac{\alpha_{i,j}^{\prime}}{\lambda_{i,i+1}}, & \textrm{if }i\in[1:N-1],j\in[0:i-1]\\ \alpha_{N,j}^{\prime}, & \textrm{if }i=N. \end{cases} $$
# Recover α_{i,j} from α′_{i,j}
# -----------------------------
α_opt = OffsetArray(zeros(size(α′_opt)), 1:N, 0:N-1)
for i in 1:N-1
for j in 0:i-1
α_opt[i,j] = (α′_opt[i,j] / λ_opt_i_ip1[i])
end
end
for j in 0:N-1
α_opt[N,j] = α′_opt[N,j]
end
We now, recover $h_{i,j}$ from $\alpha_{i,j}$ using the formula $$ \left(\forall i\in[1:N]\right)\;\left(\forall j\in[0:i-1]\right)\quad h_{i,j}=\begin{cases} \alpha_{i,j}, & \textrm{if }j=i-1\\ \alpha_{i,j}-\alpha_{i-1,j}, & \textrm{if }j\in[0:i-1]. \end{cases} $$
# Recover h_{i,j} from α_{i,j}
h_opt = OffsetArray(zeros(size(α_opt)), 1:N, 0:N-1)
# set h(1,0)
h_opt[1,0] = α_opt[1,0]
for i in 2:N
for j in 0:i-1
if j == i-1
h_opt[i,j] = α_opt[i,j]
else
h_opt[i,j] = α_opt[i,j] - α_opt[i-1,j]
end
end
end
Putting everything in a function. Now, let us put everything in a function. We need to add the packages first.
function full_pep_solver(N, L, μ, c_w, c_f)
# define all the bold vectors
# --------------------------
# define 𝐰_0
𝐰_0 = e_i(N+2, 1)
𝐰_star = zeros(N+2, 1)
𝐠_star = zeros(N+2,1)
# define 𝐠_0, 𝐠_1, …, 𝐠_N
# first we define 𝐠_Julia vectors and then 𝐠 vectors
𝐠 = OffsetArray(zeros(N+2, N+1), 1:N+2, 0:N)
# 𝐠= [𝐠_0 𝐠_1 𝐠_2 ... 𝐠_N]
for i in 0:N
𝐠[:,i] = e_i(N+2, i+2)
end
𝐟_star = zeros(N+1,1)
# time to define 𝐟_0, 𝐟_1, …, 𝐟_N
𝐟 = OffsetArray(zeros(N+1, N+1), 1:N+1, 0:N)
# 𝐟 = [𝐟_0, 𝐟_1, …, 𝐟_N]
for i in 0:N
𝐟[:,i] = e_i(N+1, i+1)
end
# Define JuMP model
# ----------------
pep_model = Model(optimizer_with_attributes(Mosek.Optimizer, "MSK_DPAR_INTPNT_CO_TOL_PFEAS" => 1e-10))
#define all the variables
# -------------------------------
# define α′ (can be typed as \alpha[TAB]\prime[TAB])
@variable(pep_model, α′[1:N, 0:N-1])
# define λ variables
@variable(pep_model, λ_i_ip1[0:N-1] >= 0)
# this defines (λ_{i,i+1})_{i∈[0:N-1]} in Julia indexing
@variable(pep_model, λ_star_i[0:N] >= 0)
# this defines (λ_{⋆,i})_{i∈[0:N]} in Julia indexing
# define τ
@variable(pep_model, τ >= 0)
# define objective
# ------------------
@objective(
pep_model,
Min,
τ
)
# Add the linear equality constraint
# ----------------------------------
@constraint(pep_model,
τ * c_f .* 𝐟[:,0] + sum(λ_i_ip1[i] .* (𝐟[:,i+1]-𝐟[:,i]) for i in 0:N-1)
+ sum(λ_star_i[i] .* (𝐟[:,i] - 𝐟_star) for i in 0:N)
.== 𝐟[:,N] - 𝐟_star
)
# Add the giant LMI constraint step by step
# ----------------------------------------
# Define all the terms one by one
term_1 = c_w * τ * ⊙(𝐰_0,𝐰_0)
term_2_part_1 = sum(λ_i_ip1[i] .* ⊙(𝐠[:,i]-𝐠[:,i+1],𝐠[:,i]-𝐠[:,i+1]) for i in 0:N-1)
term_2_part_2 = sum(λ_star_i[i] .* ⊙(𝐠_star - 𝐠[:,i],𝐠_star - 𝐠[:,i]) for i in 0:N)
term_2 = (term_2_part_1 + term_2_part_2) ./ (2*L)
term_3 = - λ_star_i[0] .* ⊙(𝐠[:,0],𝐰_0)
term_4_part_1 = - sum( λ_i_ip1[i] .* ⊙(𝐠[:,i],𝐰_0) for i in 1:N-1)
term_4_part_2 = (1/L) .*
sum( sum(α′[i,j] .* ⊙(𝐠[:,i],𝐠[:,j]) for j in 0:i-1) for i in 1:N-1)
term_4 = term_4_part_1 + term_4_part_2
term_5 = - ( ⊙(𝐠[:,N],𝐰_0) - (1/L) .* sum(α′[N,j] .* ⊙(𝐠[:,N],𝐠[:,j]) for j in 0:N-1))
term_6_part_1 = sum(λ_i_ip1[i] .* ⊙(𝐠[:,i+1],𝐰_0) for i in 0:N-1)
term_6_part_2 = - (1/L) .*
sum( sum( α′[i,j] .* ⊙(𝐠[:,i+1],𝐠[:,j]) for j in 0:i-1) for i in 1:N-1)
term_6 = term_6_part_1 + term_6_part_2
# oof, okay constructed the terms for the LMI constraint, let us hope that there is no mistake and add them together
S_mat = term_1 + term_2 + term_3 + term_4 + term_5 + term_6
# add the LMI constraint now
@constraint(pep_model,
S_mat in PSDCone()
)
# time to optimize!
# -----------------
optimize!(pep_model)
println("termination status =", termination_status(pep_model) )
α′_opt = OffsetArray(value.(α′), 1:N, 0:N-1)
λ_opt_i_ip1 = OffsetVector(value.(λ_i_ip1), 0:N-1)
λ_opt_star_i = OffsetVector(value.(λ_star_i), 0:N)
τ_opt = value(τ)
# Recover α_{i,j} from α′_{i,j}
# -----------------------------
α_opt = OffsetArray(zeros(size(α′_opt)), 1:N, 0:N-1)
for i in 1:N-1
for j in 0:i-1
α_opt[i,j] = (α′_opt[i,j] / λ_opt_i_ip1[i])
end
end
for j in 0:N-1
α_opt[N,j] = α′_opt[N,j]
end
# Recover h_{i,j} from α_{i,j}
# ----------------------------
h_opt = OffsetArray(zeros(size(α_opt)), 1:N, 0:N-1)
# set h(1,0)
h_opt[1,0] = α_opt[1,0]
# set the rest of the variables
for i in 2:N
for j in 0:i-1
if j == i-1
h_opt[i,j] = α_opt[i,j]
else
h_opt[i,j] = α_opt[i,j] - α_opt[i-1,j]
end
end
end
# return all the outputs
# ----------------------
return α′_opt, λ_opt_i_ip1, λ_opt_star_i, τ_opt, α_opt, h_opt
end
full_pep_solver (generic function with 1 method)
Time to run and test.
α′_opt, λ_opt_i_ip1, λ_opt_star_i, τ_opt, α_opt, h_opt = full_pep_solver(N, L, μ, c_w, c_f)
Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 34 Cones : 0 Scalar variables : 37 Matrix variables : 1 Integer variables : 0 Optimizer started. Presolve started. Linear dependency checker started. Linear dependency checker terminated. Eliminator started. Freed constraints in eliminator : 5 Eliminator terminated. Eliminator started. Freed constraints in eliminator : 0 Eliminator terminated. Eliminator - tries : 2 time : 0.00 Lin. dep. - tries : 1 time : 0.00 Lin. dep. - number : 0 Presolve terminated. Time: 0.00 Problem Name : Objective sense : min Type : CONIC (conic optimization problem) Constraints : 34 Cones : 0 Scalar variables : 37 Matrix variables : 1 Integer variables : 0 Optimizer - threads : 4 Optimizer - solved problem : the primal Optimizer - Constraints : 29 Optimizer - Cones : 1 Optimizer - Scalar variables : 23 conic : 16 Optimizer - Semi-definite variables: 1 scalarized : 28 Factor - setup time : 0.00 dense det. time : 0.00 Factor - ML order time : 0.00 GP order time : 0.00 Factor - nonzeros before factor : 423 after factor : 435 Factor - dense dim. : 0 flops : 1.21e+04 ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME 0 1.0e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.00 1 2.7e-01 2.7e-01 1.2e-01 2.14e+00 -2.378627944e-02 -6.117274954e-03 2.7e-01 0.00 2 1.2e-01 1.2e-01 2.9e-02 1.32e+00 7.175923608e-02 7.139302023e-02 1.2e-01 0.00 3 4.2e-02 4.2e-02 6.9e-03 1.05e+00 4.477758503e-02 4.729110878e-02 4.2e-02 0.00 4 2.3e-02 2.3e-02 2.9e-03 7.37e-01 2.880293993e-02 3.057913998e-02 2.3e-02 0.00 5 6.7e-03 6.7e-03 4.0e-04 9.85e-01 3.356646000e-02 3.302278094e-02 6.7e-03 0.00 6 2.5e-03 2.5e-03 1.1e-04 6.92e-01 2.291866424e-02 2.288414204e-02 2.5e-03 0.00 7 6.6e-04 6.6e-04 1.5e-05 8.71e-01 2.066303721e-02 2.066465445e-02 6.6e-04 0.00 8 9.8e-05 9.8e-05 8.8e-07 8.83e-01 1.883068878e-02 1.882786667e-02 9.8e-05 0.00 9 1.7e-06 1.7e-06 1.9e-09 9.93e-01 1.859217137e-02 1.859212284e-02 1.7e-06 0.02 10 4.0e-08 4.0e-08 7.3e-12 1.00e+00 1.858820828e-02 1.858820710e-02 4.0e-08 0.02 11 2.1e-09 2.1e-09 8.5e-14 1.00e+00 1.858813924e-02 1.858813917e-02 2.1e-09 0.02 12 7.9e-11 9.5e-11 6.3e-16 1.00e+00 1.858813673e-02 1.858813673e-02 7.9e-11 0.02 Optimizer terminated. Time: 0.02 termination status =OPTIMAL
([0.3149624409684461 0.0 … 0.0 0.0; 0.641151089740302 0.7224418141396977 … 0.0 0.0; … ; 1.5400244548738986 2.176849466814128 … 1.9095083874281356 0.0; 1.9256474462514417 2.800800572131717 … 2.969891147070692 2.0777698575076413], [0.19465749455910497, 0.3577518196167857, 0.5622058085120105, 0.8071885034449748, #undef], [0.12030494777139428, 0.16309432505765467, 0.20445398889522481, 0.24498269493296432, 0.19281149655502516, #undef], 0.018588136733035998, [1.618033981593307 0.0 … 0.0 0.0; 1.7921672360103886 2.019393821430617 … 0.0 0.0; … ; 1.9078870032232569 2.696829126683073 … 2.365628820626909 0.0; 1.9256474462514417 2.800800572131717 … 2.969891147070692 2.0777698575076413], [1.618033981593307 0.0 … 0.0 0.0; 0.17413325441708172 2.019393821430617 … 0.0 0.0; … ; 0.04013848423238109 0.23497477381004428 … 2.365628820626909 0.0; 0.017760443028184802 0.10397144544864378 … 0.604262326443783 2.0777698575076413])