{ "cells": [ { "cell_type": "markdown", "source": [ "**Optimization problem.** The optimization problem is\n", "\n", "\n", "$$\n", "\\begin{array}{ll}\n", "\n", "\\textrm{minimize} & f(x)\\\\\n", "\n", "\\textrm{subject to} & x\\in\\mathbf{R}^{d},\n", "\n", "\\end{array}\n", "$$\n", "\n", "\n", "where $f\\in\\mathcal{F}_{0,L}$ is the class of $L$-smooth convex functions.\n", "\n", "**Fixed-step first-order algorithm.** The optimization algorithm in consideration is:\n", "\n", "\n", "$$\n", "\\begin{align*}\n", "\n", "\\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})\\\\\n", "\n", " & =w_{0}-\\sum_{i\\in[0:i-1]}\\frac{\\alpha_{i,j}}{L}\\nabla f(w_{j}),\n", "\n", "\\end{align*}\n", "$$\n", "\n", "\n", "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:\n", "\n", "\n", "$$\n", "\\left(\\forall i\\in[1:N]\\right)\\;\\left(\\forall j\\in[0:i-1]\\right)\\quad h_{i,j}=\\begin{cases}\n", "\n", "\\alpha_{i,j}, & \\textrm{if }j=i-1\\\\\n", "\n", "\\alpha_{i,j}-\\alpha_{i-1,j}, & \\textrm{if }j\\in[0:i-1].\n", "\n", "\\end{cases}\n", "$$\n", "\n", "\n", "\n", "**Notation.** Denote,\n", "$$\n", "\\begin{align*}\n", "J & =\\{(i,j)\\mid j=i+1,i\\in[0:N-1]\\}\\cup\\{(i,j)\\mid i=\\star,j\\in[0:N]\\},\n", "\\end{align*}\n", "$$\n", " and\n", "$$\n", "\\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},\n", "$$\n", "where $e_{i}$ is the unit vector with $i$th component equal to $1$ and the rest being zero. Next, define\n", "$$\n", "\\begin{align*}\n", " & S\\left(\\tau,\\{\\lambda_{i,j}\\},\\{\\alpha_{i,j}^{\\prime}\\}\\right)\\\\\n", "= & c_{w}\\tau\\mathbf{w}_{0}\\mathbf{w}_{0}^{\\top}+\\\\\n", " & \\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\\} +\\\\\n", " & \\sum_{j\\in[0:N]}\\lambda_{\\star,j}\\left\\{ (\\mathbf{g}_{\\star}-\\mathbf{g}_{j})\\odot(\\mathbf{g}_{\\star}-\\mathbf{g}_{j})\\right\\} \\Bigg]\\\\\n", " & -\\lambda_{\\star,0}\\{\\mathbf{g}_{0}\\odot\\mathbf{w}_{0}\\}\\\\\n", " & -\\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]\\\\\n", " & -\\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)\\\\\n", " & +\\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],\n", "\\end{align*}\n", "$$\n", "where\n", "$$\n", "\\alpha_{i,j}^{\\prime}=\\begin{cases}\n", "\\lambda_{i,i+1}\\alpha_{i,j}, & \\textrm{if }i\\in[1:N-1],j\\in[0:i-1]\\\\\n", "\\alpha_{N,j}, & \\textrm{if }i\\in N,j\\in[0:N-1]\\\\\n", "\\lambda_{0,1}\\underbrace{\\alpha_{0,j}}_{=0}=0, & \\textrm{if }i=0.\n", "\\end{cases}\n", "$$\n", "Then, the performance estimation optimization algorithm is:\n", "$$\n", "\\begin{array}{ll}\n", "\\textrm{minimize} & \\tau\\\\\n", "\\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\\\\\n", " & S\\left(\\tau,\\{\\lambda_{i,j}\\},\\{\\alpha_{i,j}^{\\prime}\\}\\right)\\succeq0\\\\\n", " & \\left(\\forall(i,j)\\in J\\right)\\quad\\lambda_{i,j}\\geq0\\\\\n", " & \\tau\\geq0,\n", "\\end{array}\n", "$$\n", "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`." ], "metadata": {} }, { "cell_type": "code", "source": [ "# Load the packages\n", "using JuMP, MosekTools, Mosek, LinearAlgebra, SCS, COSMO, Literate, OffsetArrays" ], "outputs": [], "execution_count": 1, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:26:56.917Z", "iopub.execute_input": "2021-04-23T13:26:57.487Z", "iopub.status.idle": "2021-04-23T13:27:00.662Z" } } }, { "cell_type": "code", "source": [ "# %% Some helper functions\n", "\n", "# %% construct e_i in R^n\n", "function e_i(n, i)\n", " e_i_vec = zeros(n, 1)\n", " e_i_vec[i] = 1\n", " return e_i_vec\n", "end\n", "\n", "# %% construct symmetric outer product\n", "\n", "function ⊙(a,b)\n", " return ((a*b')+(b*a')) ./ 2\n", "end" ], "outputs": [ { "output_type": "execute_result", "execution_count": 2, "data": { "text/plain": "⊙ (generic function with 1 method)" }, "metadata": {} } ], "execution_count": 2, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:00.682Z", "iopub.execute_input": "2021-04-23T13:27:01.025Z", "iopub.status.idle": "2021-04-23T13:27:02.260Z" } } }, { "cell_type": "code", "source": [ "# %% Parameters to be tuned\n", "N = 5\n", "L = 1\n", "μ = 0\n", "c_w = 1\n", "c_f = 0" ], "outputs": [ { "output_type": "execute_result", "execution_count": 3, "data": { "text/plain": "0" }, "metadata": {} } ], "execution_count": 3, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:02.283Z", "iopub.execute_input": "2021-04-23T13:27:02.300Z", "iopub.status.idle": "2021-04-23T13:27:02.422Z" } } }, { "cell_type": "markdown", "source": [ "Next, define the bold vectors:\n", "\n", "$$\n", "\\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},\n", "$$\n", "where $e_{i}$ is the unit vector with $i$th component equal to $1$ and the rest being zero." ], "metadata": {} }, { "cell_type": "code", "source": [ "# define all the bold vectors\n", "# --------------------------\n", "\n", "# define 𝐰_0\n", "\n", "𝐰_0 = e_i(N+2, 1)\n", "\n", "𝐰_star = zeros(N+2, 1)\n", "\n", "𝐠_star = zeros(N+2,1)\n", "\n", "# define 𝐠_0, 𝐠_1, …, 𝐠_N\n", "\n", "# first we define 𝐠_Julia vectors and then 𝐠 vectors\n", "\n", "𝐠 = OffsetArray(zeros(N+2, N+1), 1:N+2, 0:N)\n", "# 𝐠= [𝐠_0 𝐠_1 𝐠_2 ... 𝐠_N]\n", "for i in 0:N\n", " 𝐠[:,i] = e_i(N+2, i+2)\n", "end\n", "\n", "\n", "𝐟_star = zeros(N+1,1)\n", "\n", "# time to define 𝐟_0, 𝐟_1, …, 𝐟_N\n", "\n", "𝐟 = OffsetArray(zeros(N+1, N+1), 1:N+1, 0:N)\n", "# 𝐟 = [𝐟_0, 𝐟_1, …, 𝐟_N]\n", "\n", "for i in 0:N\n", " 𝐟[:,i] = e_i(N+1, i+1)\n", "end" ], "outputs": [], "execution_count": 4, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:02.443Z", "iopub.execute_input": "2021-04-23T13:27:02.459Z", "iopub.status.idle": "2021-04-23T13:27:02.671Z" } } }, { "cell_type": "markdown", "source": [ "Next, we define our decision variables using `JuMP`." ], "metadata": {} }, { "cell_type": "code", "source": [ "pep_model = Model(optimizer_with_attributes(Mosek.Optimizer, \"MSK_DPAR_INTPNT_CO_TOL_PFEAS\" => 1e-10))\n", "\n", "# time to define all the variables\n", "# -------------------------------\n", "\n", "# define α′ (can be typed as \\alpha[TAB]\\prime[TAB])\n", "@variable(pep_model, α′[1:N, 0:N-1])\n", "\n", "# define λ variables\n", "\n", "@variable(pep_model, λ_i_ip1[0:N-1] >= 0)\n", "# this defines (λ_{i,i+1})_{i∈[0:N-1]} in Julia indexing\n", "\n", "@variable(pep_model, λ_star_i[0:N] >= 0)\n", "# this defines (λ_{⋆,i})_{i∈[0:N]} in Julia indexing\n", "\n", "# define τ\n", "@variable(pep_model, τ >= 0)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 5, "data": { "text/plain": "τ", "text/latex": "$$ τ $$" }, "metadata": {} } ], "execution_count": 5, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:02.692Z", "iopub.execute_input": "2021-04-23T13:27:02.706Z", "iopub.status.idle": "2021-04-23T13:27:07.172Z" } } }, { "cell_type": "markdown", "source": [ "Define the objective next, which is to minimize $\\tau$." ], "metadata": {} }, { "cell_type": "code", "source": [ "# define objective\n", "@objective(\n", " pep_model,\n", " Min,\n", " τ\n", " )" ], "outputs": [ { "output_type": "execute_result", "execution_count": 6, "data": { "text/plain": "τ", "text/latex": "$$ τ $$" }, "metadata": {} } ], "execution_count": 6, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:07.196Z", "iopub.execute_input": "2021-04-23T13:27:07.213Z", "iopub.status.idle": "2021-04-23T13:27:07.283Z" } } }, { "cell_type": "markdown", "source": [ "Time to add the linear constraint\n", "$$\n", "\\begin{align*}\n", " & -\\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\\\\\n", "\\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}\n", "\\end{align*}\n", "$$\n", " first, which are much simpler." ], "metadata": {} }, { "cell_type": "code", "source": [ "# Add the linear equality constraint\n", "# ----------------------------------\n", "@constraint(pep_model,\n", " τ * c_f .* 𝐟[:,0] + sum(λ_i_ip1[i] .* (𝐟[:,i+1]-𝐟[:,i]) for i in 0:N-1)\n", " + sum(λ_star_i[i] .* (𝐟[:,i] - 𝐟_star) for i in 0:N)\n", " .== 𝐟[:,N] - 𝐟_star\n", ")" ], "outputs": [ { "output_type": "execute_result", "execution_count": 7, "data": { "text/plain": "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:\n -λ_i_ip1[0] + λ_star_i[0] == 0.0\n λ_i_ip1[0] - λ_i_ip1[1] + λ_star_i[1] == 0.0\n λ_i_ip1[1] - λ_i_ip1[2] + λ_star_i[2] == 0.0\n λ_i_ip1[2] - λ_i_ip1[3] + λ_star_i[3] == 0.0\n λ_i_ip1[3] - λ_i_ip1[4] + λ_star_i[4] == 0.0\n λ_i_ip1[4] + λ_star_i[5] == 1.0" }, "metadata": {} } ], "execution_count": 7, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:07.309Z", "iopub.execute_input": "2021-04-23T13:27:07.322Z", "iopub.status.idle": "2021-04-23T13:27:11.778Z" } } }, { "cell_type": "markdown", "source": [ "Now let us construct the giant sdp constraint step by step. It has 6 summands, which we add one by one." ], "metadata": {} }, { "cell_type": "code", "source": [ "# Add the giant LMI constraint step by step\n", "# ----------------------------------------" ], "outputs": [], "execution_count": 8, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:11.799Z", "iopub.execute_input": "2021-04-23T13:27:11.815Z", "iopub.status.idle": "2021-04-23T13:27:11.835Z" } } }, { "cell_type": "markdown", "source": [ "Term 1 is $\\texttt{term}_1 = c_{w}\\tau\\mathbf{w}_{0}\\mathbf{w}_{0}^{\\top}.$" ], "metadata": {} }, { "cell_type": "code", "source": [ "term_1 = c_w * τ * ⊙(𝐰_0,𝐰_0)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 9, "data": { "text/plain": "7×7 Matrix{AffExpr}:\n τ 0 0 0 0 0 0\n 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0" }, "metadata": {} } ], "execution_count": 9, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:11.856Z", "iopub.execute_input": "2021-04-23T13:27:11.869Z", "iopub.status.idle": "2021-04-23T13:27:13.610Z" } } }, { "cell_type": "markdown", "source": [ "Term 2 is \n", "$$\n", "\\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].\n", "$$" ], "metadata": {} }, { "cell_type": "code", "source": [ "term_2_part_1 = sum(λ_i_ip1[i] .* ⊙(𝐠[:,i]-𝐠[:,i+1],𝐠[:,i]-𝐠[:,i+1]) for i in 0:N-1)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 10, "data": { "text/plain": "7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7:\n 0 0 0 … 0\n 0 λ_i_ip1[0] -λ_i_ip1[0] 0\n 0 -λ_i_ip1[0] λ_i_ip1[0] + λ_i_ip1[1] 0\n 0 0 -λ_i_ip1[1] 0\n 0 0 0 0\n 0 0 0 … -λ_i_ip1[4]\n 0 0 0 λ_i_ip1[4]" }, "metadata": {} } ], "execution_count": 10, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:13.631Z", "iopub.execute_input": "2021-04-23T13:27:13.645Z", "iopub.status.idle": "2021-04-23T13:27:14.553Z" } } }, { "cell_type": "code", "source": [ "term_2_part_2 = sum(λ_star_i[i] .* ⊙(𝐠_star - 𝐠[:,i],𝐠_star - 𝐠[:,i]) for i in 0:N)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 11, "data": { "text/plain": "7×7 Matrix{AffExpr}:\n 0 0 0 0 … 0 0\n 0 λ_star_i[0] 0 0 0 0\n 0 0 λ_star_i[1] 0 0 0\n 0 0 0 λ_star_i[2] 0 0\n 0 0 0 0 0 0\n 0 0 0 0 … λ_star_i[4] 0\n 0 0 0 0 0 λ_star_i[5]" }, "metadata": {} } ], "execution_count": 11, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:14.573Z", "iopub.execute_input": "2021-04-23T13:27:14.585Z", "iopub.status.idle": "2021-04-23T13:27:15.428Z" } } }, { "cell_type": "code", "source": [ "term_2 = (term_2_part_1 + term_2_part_2) ./ (2*L)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 12, "data": { "text/plain": "7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7:\n 0 0 … 0\n 0 0.5 λ_i_ip1[0] + 0.5 λ_star_i[0] 0\n 0 -0.5 λ_i_ip1[0] 0\n 0 0 0\n 0 0 0\n 0 0 … -0.5 λ_i_ip1[4]\n 0 0 0.5 λ_i_ip1[4] + 0.5 λ_star_i[5]" }, "metadata": {} } ], "execution_count": 12, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:15.450Z", "iopub.execute_input": "2021-04-23T13:27:15.463Z", "iopub.status.idle": "2021-04-23T13:27:15.764Z" } } }, { "cell_type": "markdown", "source": [ "Term 3 is $\\texttt{term}_3 = -\\lambda_{\\star,0}\\{\\mathbf{g}_{0}\\odot\\mathbf{w}_{0}\\}$." ], "metadata": {} }, { "cell_type": "code", "source": [ "term_3 = - λ_star_i[0] .* ⊙(𝐠[:,0],𝐰_0)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 13, "data": { "text/plain": "7×7 Matrix{AffExpr}:\n 0 -0.5 λ_star_i[0] 0 0 0 0 0\n -0.5 λ_star_i[0] 0 0 0 0 0 0\n 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0\n 0 0 0 0 0 0 0" }, "metadata": {} } ], "execution_count": 13, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:15.783Z", "iopub.execute_input": "2021-04-23T13:27:15.796Z", "iopub.status.idle": "2021-04-23T13:27:16.237Z" } } }, { "cell_type": "markdown", "source": [ "Term 4 is \n", "$$\n", "\\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].\n", "$$" ], "metadata": {} }, { "cell_type": "code", "source": [ "term_4_part_1 = - sum( λ_i_ip1[i] .* ⊙(𝐠[:,i],𝐰_0) for i in 1:N-1)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 14, "data": { "text/plain": "7×7 Matrix{AffExpr}:\n 0 0 -0.5 λ_i_ip1[1] … -0.5 λ_i_ip1[3] -0.5 λ_i_ip1[4] 0\n 0 0 0 0 0 0\n -0.5 λ_i_ip1[1] 0 0 0 0 0\n -0.5 λ_i_ip1[2] 0 0 0 0 0\n -0.5 λ_i_ip1[3] 0 0 0 0 0\n -0.5 λ_i_ip1[4] 0 0 … 0 0 0\n 0 0 0 0 0 0" }, "metadata": {} } ], "execution_count": 14, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:16.256Z", "iopub.execute_input": "2021-04-23T13:27:16.269Z", "iopub.status.idle": "2021-04-23T13:27:16.371Z" } } }, { "cell_type": "code", "source": [ "term_4_part_2 = (1/L) .*\n", "sum( sum(α′[i,j] .* ⊙(𝐠[:,i],𝐠[:,j]) for j in 0:i-1) for i in 1:N-1)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 15, "data": { "text/plain": "7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7:\n 0 0 0 0 0 0 0\n 0 0 0.5 α′[1,0] 0.5 α′[2,0] 0.5 α′[3,0] 0.5 α′[4,0] 0\n 0 0.5 α′[1,0] 0 0.5 α′[2,1] 0.5 α′[3,1] 0.5 α′[4,1] 0\n 0 0.5 α′[2,0] 0.5 α′[2,1] 0 0.5 α′[3,2] 0.5 α′[4,2] 0\n 0 0.5 α′[3,0] 0.5 α′[3,1] 0.5 α′[3,2] 0 0.5 α′[4,3] 0\n 0 0.5 α′[4,0] 0.5 α′[4,1] 0.5 α′[4,2] 0.5 α′[4,3] 0 0\n 0 0 0 0 0 0 0" }, "metadata": {} } ], "execution_count": 15, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:16.392Z", "iopub.execute_input": "2021-04-23T13:27:16.406Z", "iopub.status.idle": "2021-04-23T13:27:16.594Z" } } }, { "cell_type": "code", "source": [ "term_4 = term_4_part_1 + term_4_part_2" ], "outputs": [ { "output_type": "execute_result", "execution_count": 16, "data": { "text/plain": "7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7:\n 0 0 … -0.5 λ_i_ip1[3] -0.5 λ_i_ip1[4] 0\n 0 0 0.5 α′[3,0] 0.5 α′[4,0] 0\n -0.5 λ_i_ip1[1] 0.5 α′[1,0] 0.5 α′[3,1] 0.5 α′[4,1] 0\n -0.5 λ_i_ip1[2] 0.5 α′[2,0] 0.5 α′[3,2] 0.5 α′[4,2] 0\n -0.5 λ_i_ip1[3] 0.5 α′[3,0] 0 0.5 α′[4,3] 0\n -0.5 λ_i_ip1[4] 0.5 α′[4,0] … 0.5 α′[4,3] 0 0\n 0 0 0 0 0" }, "metadata": {} } ], "execution_count": 16, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:16.616Z", "iopub.execute_input": "2021-04-23T13:27:16.630Z", "iopub.status.idle": "2021-04-23T13:27:16.715Z" } } }, { "cell_type": "markdown", "source": [ "Term 5 is \n", "$$\n", "\\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).\n", "$$" ], "metadata": {} }, { "cell_type": "code", "source": [ "term_5 = - ( ⊙(𝐠[:,N],𝐰_0) - (1/L) .* sum(α′[N,j] .* ⊙(𝐠[:,N],𝐠[:,j]) for j in 0:N-1))" ], "outputs": [ { "output_type": "execute_result", "execution_count": 17, "data": { "text/plain": "7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7:\n 0 0 0 0 … 0 -0.5\n 0 0 0 0 0 0.5 α′[5,0]\n 0 0 0 0 0 0.5 α′[5,1]\n 0 0 0 0 0 0.5 α′[5,2]\n 0 0 0 0 0 0.5 α′[5,3]\n 0 0 0 0 … 0 0.5 α′[5,4]\n -0.5 0.5 α′[5,0] 0.5 α′[5,1] 0.5 α′[5,2] 0.5 α′[5,4] 0" }, "metadata": {} } ], "execution_count": 17, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:16.735Z", "iopub.execute_input": "2021-04-23T13:27:16.748Z", "iopub.status.idle": "2021-04-23T13:27:16.909Z" } } }, { "cell_type": "markdown", "source": [ "Okay, we are almost there. One more term. Term 6 is:\n", "$$\n", "\\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].\n", "$$" ], "metadata": {} }, { "cell_type": "code", "source": [ "term_6_part_1 = sum(λ_i_ip1[i] .* ⊙(𝐠[:,i+1],𝐰_0) for i in 0:N-1)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 18, "data": { "text/plain": "7×7 Matrix{AffExpr}:\n 0 0 0.5 λ_i_ip1[0] … 0.5 λ_i_ip1[3] 0.5 λ_i_ip1[4]\n 0 0 0 0 0\n 0.5 λ_i_ip1[0] 0 0 0 0\n 0.5 λ_i_ip1[1] 0 0 0 0\n 0.5 λ_i_ip1[2] 0 0 0 0\n 0.5 λ_i_ip1[3] 0 0 … 0 0\n 0.5 λ_i_ip1[4] 0 0 0 0" }, "metadata": {} } ], "execution_count": 18, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:16.930Z", "iopub.execute_input": "2021-04-23T13:27:16.943Z", "iopub.status.idle": "2021-04-23T13:27:16.970Z" } } }, { "cell_type": "code", "source": [ "term_6_part_2 = - (1/L) .*\n", "sum( sum( α′[i,j] .* ⊙(𝐠[:,i+1],𝐠[:,j]) for j in 0:i-1) for i in 1:N-1)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 19, "data": { "text/plain": "7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7:\n 0 0 0 0 … 0 0\n 0 0 0 -0.5 α′[1,0] -0.5 α′[3,0] -0.5 α′[4,0]\n 0 0 0 0 -0.5 α′[3,1] -0.5 α′[4,1]\n 0 -0.5 α′[1,0] 0 0 -0.5 α′[3,2] -0.5 α′[4,2]\n 0 -0.5 α′[2,0] -0.5 α′[2,1] 0 0 -0.5 α′[4,3]\n 0 -0.5 α′[3,0] -0.5 α′[3,1] -0.5 α′[3,2] … 0 0\n 0 -0.5 α′[4,0] -0.5 α′[4,1] -0.5 α′[4,2] 0 0" }, "metadata": {} } ], "execution_count": 19, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:16.992Z", "iopub.execute_input": "2021-04-23T13:27:17.005Z", "iopub.status.idle": "2021-04-23T13:27:17.091Z" } } }, { "cell_type": "code", "source": [ "term_6 = term_6_part_1 + term_6_part_2" ], "outputs": [ { "output_type": "execute_result", "execution_count": 20, "data": { "text/plain": "7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7:\n 0 0 … 0.5 λ_i_ip1[3] 0.5 λ_i_ip1[4]\n 0 0 -0.5 α′[3,0] -0.5 α′[4,0]\n 0.5 λ_i_ip1[0] 0 -0.5 α′[3,1] -0.5 α′[4,1]\n 0.5 λ_i_ip1[1] -0.5 α′[1,0] -0.5 α′[3,2] -0.5 α′[4,2]\n 0.5 λ_i_ip1[2] -0.5 α′[2,0] 0 -0.5 α′[4,3]\n 0.5 λ_i_ip1[3] -0.5 α′[3,0] … 0 0\n 0.5 λ_i_ip1[4] -0.5 α′[4,0] 0 0" }, "metadata": {} } ], "execution_count": 20, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:17.111Z", "iopub.execute_input": "2021-04-23T13:27:17.124Z", "iopub.status.idle": "2021-04-23T13:27:17.153Z" } } }, { "cell_type": "markdown", "source": [ "Time to add all these terms to construct $S\\left(\\tau,\\{\\lambda_{i,j}\\},\\{\\alpha_{i,j}^{\\prime}\\}\\right)$:\n", "$$\n", "S\\left(\\tau,\\{\\lambda_{i,j}\\},\\{\\alpha_{i,j}^{\\prime}\\}\\right) = \\sum_{i\\in [1:6]}\\texttt{term}_i\n", "$$" ], "metadata": {} }, { "cell_type": "code", "source": [ "# oof, okay constructed the terms for the LMI constraint, let us hope that there is no mistake and add them together\n", "\n", "S_mat = term_1 + term_2 + term_3 + term_4 + term_5 + term_6" ], "outputs": [ { "output_type": "execute_result", "execution_count": 21, "data": { "text/plain": "7×7 OffsetArray(::Matrix{AffExpr}, 1:7, 1:7) with eltype AffExpr with indices 1:7×1:7:\n τ … 0.5 λ_i_ip1[4] - 0.5\n -0.5 λ_star_i[0] 0.5 α′[5,0] - 0.5 α′[4,0]\n -0.5 λ_i_ip1[1] + 0.5 λ_i_ip1[0] 0.5 α′[5,1] - 0.5 α′[4,1]\n -0.5 λ_i_ip1[2] + 0.5 λ_i_ip1[1] 0.5 α′[5,2] - 0.5 α′[4,2]\n -0.5 λ_i_ip1[3] + 0.5 λ_i_ip1[2] 0.5 α′[5,3] - 0.5 α′[4,3]\n -0.5 λ_i_ip1[4] + 0.5 λ_i_ip1[3] … -0.5 λ_i_ip1[4] + 0.5 α′[5,4]\n 0.5 λ_i_ip1[4] - 0.5 0.5 λ_i_ip1[4] + 0.5 λ_star_i[5]" }, "metadata": {} } ], "execution_count": 21, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:17.179Z", "iopub.execute_input": "2021-04-23T13:27:17.194Z", "iopub.status.idle": "2021-04-23T13:27:17.221Z" } } }, { "cell_type": "markdown", "source": [ "Let's add the LMI cosntraint $S\\left(\\tau,\\{\\lambda_{i,j}\\},\\{\\alpha_{i,j}^{\\prime}\\}\\right)\\succeq0$." ], "metadata": {} }, { "cell_type": "code", "source": [ "# add the LMI constraint now\n", "\n", "@constraint(pep_model,\n", "S_mat in PSDCone()\n", ")" ], "outputs": [ { "output_type": "execute_result", "execution_count": 22, "data": { "text/plain": "[τ -0.5 λ_star_i[0] 0.5 λ_i_ip1[0] - 0.5 λ_i_ip1[1] 0.5 λ_i_ip1[1] - 0.5 λ_i_ip1[2] 0.5 λ_i_ip1[2] - 0.5 λ_i_ip1[3] 0.5 λ_i_ip1[3] - 0.5 λ_i_ip1[4] 0.5 λ_i_ip1[4] - 0.5;\n -0.5 λ_star_i[0] 0.5 λ_i_ip1[0] + 0.5 λ_star_i[0] 0.5 α′[1,0] - 0.5 λ_i_ip1[0] -0.5 α′[1,0] + 0.5 α′[2,0] -0.5 α′[2,0] + 0.5 α′[3,0] -0.5 α′[3,0] + 0.5 α′[4,0] -0.5 α′[4,0] + 0.5 α′[5,0];\n 0.5 λ_i_ip1[0] - 0.5 λ_i_ip1[1] 0.5 α′[1,0] - 0.5 λ_i_ip1[0] 0.5 λ_i_ip1[0] + 0.5 λ_i_ip1[1] + 0.5 λ_star_i[1] 0.5 α′[2,1] - 0.5 λ_i_ip1[1] -0.5 α′[2,1] + 0.5 α′[3,1] -0.5 α′[3,1] + 0.5 α′[4,1] -0.5 α′[4,1] + 0.5 α′[5,1];\n 0.5 λ_i_ip1[1] - 0.5 λ_i_ip1[2] -0.5 α′[1,0] + 0.5 α′[2,0] 0.5 α′[2,1] - 0.5 λ_i_ip1[1] 0.5 λ_i_ip1[1] + 0.5 λ_i_ip1[2] + 0.5 λ_star_i[2] 0.5 α′[3,2] - 0.5 λ_i_ip1[2] -0.5 α′[3,2] + 0.5 α′[4,2] -0.5 α′[4,2] + 0.5 α′[5,2]; in PSDCone()\n 0.5 λ_i_ip1[2] - 0.5 λ_i_ip1[3] -0.5 α′[2,0] + 0.5 α′[3,0] -0.5 α′[2,1] + 0.5 α′[3,1] 0.5 α′[3,2] - 0.5 λ_i_ip1[2] 0.5 λ_i_ip1[2] + 0.5 λ_i_ip1[3] + 0.5 λ_star_i[3] 0.5 α′[4,3] - 0.5 λ_i_ip1[3] -0.5 α′[4,3] + 0.5 α′[5,3];\n 0.5 λ_i_ip1[3] - 0.5 λ_i_ip1[4] -0.5 α′[3,0] + 0.5 α′[4,0] -0.5 α′[3,1] + 0.5 α′[4,1] -0.5 α′[3,2] + 0.5 α′[4,2] 0.5 α′[4,3] - 0.5 λ_i_ip1[3] 0.5 λ_i_ip1[3] + 0.5 λ_i_ip1[4] + 0.5 λ_star_i[4] 0.5 α′[5,4] - 0.5 λ_i_ip1[4];\n 0.5 λ_i_ip1[4] - 0.5 -0.5 α′[4,0] + 0.5 α′[5,0] -0.5 α′[4,1] + 0.5 α′[5,1] -0.5 α′[4,2] + 0.5 α′[5,2] -0.5 α′[4,3] + 0.5 α′[5,3] 0.5 α′[5,4] - 0.5 λ_i_ip1[4] 0.5 λ_i_ip1[4] + 0.5 λ_star_i[5]]", "text/latex": "$$ \\begin{bmatrix}\nτ & -0.5 λ\\_star\\_i_{0} & 0.5 λ\\_i\\_ip1_{0} - 0.5 λ\\_i\\_ip1_{1} & 0.5 λ\\_i\\_ip1_{1} - 0.5 λ\\_i\\_ip1_{2} & 0.5 λ\\_i\\_ip1_{2} - 0.5 λ\\_i\\_ip1_{3} & 0.5 λ\\_i\\_ip1_{3} - 0.5 λ\\_i\\_ip1_{4} & 0.5 λ\\_i\\_ip1_{4} - 0.5\\\\\n-0.5 λ\\_star\\_i_{0} & 0.5 λ\\_i\\_ip1_{0} + 0.5 λ\\_star\\_i_{0} & 0.5 α′_{1,0} - 0.5 λ\\_i\\_ip1_{0} & -0.5 α′_{1,0} + 0.5 α′_{2,0} & -0.5 α′_{2,0} + 0.5 α′_{3,0} & -0.5 α′_{3,0} + 0.5 α′_{4,0} & -0.5 α′_{4,0} + 0.5 α′_{5,0}\\\\\n0.5 λ\\_i\\_ip1_{0} - 0.5 λ\\_i\\_ip1_{1} & 0.5 α′_{1,0} - 0.5 λ\\_i\\_ip1_{0} & 0.5 λ\\_i\\_ip1_{0} + 0.5 λ\\_i\\_ip1_{1} + 0.5 λ\\_star\\_i_{1} & 0.5 α′_{2,1} - 0.5 λ\\_i\\_ip1_{1} & -0.5 α′_{2,1} + 0.5 α′_{3,1} & -0.5 α′_{3,1} + 0.5 α′_{4,1} & -0.5 α′_{4,1} + 0.5 α′_{5,1}\\\\\n0.5 λ\\_i\\_ip1_{1} - 0.5 λ\\_i\\_ip1_{2} & -0.5 α′_{1,0} + 0.5 α′_{2,0} & 0.5 α′_{2,1} - 0.5 λ\\_i\\_ip1_{1} & 0.5 λ\\_i\\_ip1_{1} + 0.5 λ\\_i\\_ip1_{2} + 0.5 λ\\_star\\_i_{2} & 0.5 α′_{3,2} - 0.5 λ\\_i\\_ip1_{2} & -0.5 α′_{3,2} + 0.5 α′_{4,2} & -0.5 α′_{4,2} + 0.5 α′_{5,2}\\\\\n0.5 λ\\_i\\_ip1_{2} - 0.5 λ\\_i\\_ip1_{3} & -0.5 α′_{2,0} + 0.5 α′_{3,0} & -0.5 α′_{2,1} + 0.5 α′_{3,1} & 0.5 α′_{3,2} - 0.5 λ\\_i\\_ip1_{2} & 0.5 λ\\_i\\_ip1_{2} + 0.5 λ\\_i\\_ip1_{3} + 0.5 λ\\_star\\_i_{3} & 0.5 α′_{4,3} - 0.5 λ\\_i\\_ip1_{3} & -0.5 α′_{4,3} + 0.5 α′_{5,3}\\\\\n0.5 λ\\_i\\_ip1_{3} - 0.5 λ\\_i\\_ip1_{4} & -0.5 α′_{3,0} + 0.5 α′_{4,0} & -0.5 α′_{3,1} + 0.5 α′_{4,1} & -0.5 α′_{3,2} + 0.5 α′_{4,2} & 0.5 α′_{4,3} - 0.5 λ\\_i\\_ip1_{3} & 0.5 λ\\_i\\_ip1_{3} + 0.5 λ\\_i\\_ip1_{4} + 0.5 λ\\_star\\_i_{4} & 0.5 α′_{5,4} - 0.5 λ\\_i\\_ip1_{4}\\\\\n0.5 λ\\_i\\_ip1_{4} - 0.5 & -0.5 α′_{4,0} + 0.5 α′_{5,0} & -0.5 α′_{4,1} + 0.5 α′_{5,1} & -0.5 α′_{4,2} + 0.5 α′_{5,2} & -0.5 α′_{4,3} + 0.5 α′_{5,3} & 0.5 α′_{5,4} - 0.5 λ\\_i\\_ip1_{4} & 0.5 λ\\_i\\_ip1_{4} + 0.5 λ\\_star\\_i_{5}\\\\\n\\end{bmatrix} \\in \\text{PSDCone()} $$" }, "metadata": {} } ], "execution_count": 22, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:17.241Z", "iopub.execute_input": "2021-04-23T13:27:17.254Z", "iopub.status.idle": "2021-04-23T13:27:19.595Z" } } }, { "cell_type": "markdown", "source": [ "Time to optimize the code now!" ], "metadata": {} }, { "cell_type": "code", "source": [ "# time to optimize!\n", "# -----------------\n", "optimize!(pep_model)\n", "\n", "println(\"termination status =\", termination_status(pep_model) )" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 34 \n", " Cones : 0 \n", " Scalar variables : 37 \n", " Matrix variables : 1 \n", " Integer variables : 0 \n", "\n", "Optimizer started.\n", "Presolve started.\n", "Linear dependency checker started.\n", "Linear dependency checker terminated.\n", "Eliminator started.\n", "Freed constraints in eliminator : 5\n", "Eliminator terminated.\n", "Eliminator started.\n", "Freed constraints in eliminator : 0\n", "Eliminator terminated.\n", "Eliminator - tries : 2 time : 0.00 \n", "Lin. dep. - tries : 1 time : 0.00 \n", "Lin. dep. - number : 0 \n", "Presolve terminated. Time: 0.00 \n", "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 34 \n", " Cones : 0 \n", " Scalar variables : 37 \n", " Matrix variables : 1 \n", " Integer variables : 0 \n", "\n", "Optimizer - threads : 4 \n", "Optimizer - solved problem : the primal \n", "Optimizer - Constraints : 29\n", "Optimizer - Cones : 1\n", "Optimizer - Scalar variables : 23 conic : 16 \n", "Optimizer - Semi-definite variables: 1 scalarized : 28 \n", "Factor - setup time : 0.00 dense det. time : 0.00 \n", "Factor - ML order time : 0.00 GP order time : 0.00 \n", "Factor - nonzeros before factor : 423 after factor : 435 \n", "Factor - dense dim. : 0 flops : 1.21e+04 \n", "ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n", "0 1.0e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.00 \n", "1 2.7e-01 2.7e-01 1.2e-01 2.14e+00 -2.378627944e-02 -6.117274954e-03 2.7e-01 0.00 \n", "2 1.2e-01 1.2e-01 2.9e-02 1.32e+00 7.175923608e-02 7.139302023e-02 1.2e-01 0.00 \n", "3 4.2e-02 4.2e-02 6.9e-03 1.05e+00 4.477758503e-02 4.729110878e-02 4.2e-02 0.00 \n", "4 2.3e-02 2.3e-02 2.9e-03 7.37e-01 2.880293993e-02 3.057913998e-02 2.3e-02 0.00 \n", "5 6.7e-03 6.7e-03 4.0e-04 9.85e-01 3.356646000e-02 3.302278094e-02 6.7e-03 0.00 \n", "6 2.5e-03 2.5e-03 1.1e-04 6.92e-01 2.291866424e-02 2.288414204e-02 2.5e-03 0.00 \n", "7 6.6e-04 6.6e-04 1.5e-05 8.71e-01 2.066303721e-02 2.066465445e-02 6.6e-04 0.00 \n", "8 9.8e-05 9.8e-05 8.8e-07 8.83e-01 1.883068878e-02 1.882786667e-02 9.8e-05 0.00 \n", "9 1.7e-06 1.7e-06 1.9e-09 9.93e-01 1.859217137e-02 1.859212284e-02 1.7e-06 0.00 \n", "10 4.0e-08 4.0e-08 7.3e-12 1.00e+00 1.858820828e-02 1.858820710e-02 4.0e-08 0.00 \n", "11 2.1e-09 2.1e-09 8.5e-14 1.00e+00 1.858813924e-02 1.858813917e-02 2.1e-09 0.00 \n", "12 7.9e-11 9.5e-11 6.3e-16 1.00e+00 1.858813673e-02 1.858813673e-02 7.9e-11 0.00 \n", "Optimizer terminated. Time: 0.01 \n", "\n", "termination status =OPTIMAL\n" ] } ], "execution_count": 23, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:19.623Z", "iopub.execute_input": "2021-04-23T13:27:19.642Z", "iopub.status.idle": "2021-04-23T13:27:31.156Z" } } }, { "cell_type": "markdown", "source": [ "We now collect the optimal values of our decision variables." ], "metadata": {} }, { "cell_type": "code", "source": [ "# Collect the decision variables\n", "# -----------------------------\n", "α′_opt = OffsetArray(value.(α′), 1:N, 0:N-1)\n", "λ_opt_i_ip1 = OffsetVector(value.(λ_i_ip1), 0:N-1)\n", "λ_opt_star_i = OffsetVector(value.(λ_star_i), 0:N)\n", "τ_opt = value(τ)" ], "outputs": [ { "output_type": "execute_result", "execution_count": 24, "data": { "text/plain": "0.018588136733035998" }, "metadata": {} } ], "execution_count": 24, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:31.184Z", "iopub.execute_input": "2021-04-23T13:27:31.199Z", "iopub.status.idle": "2021-04-23T13:27:32.226Z" } } }, { "cell_type": "markdown", "source": [ "Next, we recover $\\alpha_{i,j}$ from $\\alpha^\\prime_{i,j}$ using the formula\n", "$$\n", "\\alpha_{i,j}=\\begin{cases}\n", "\\frac{\\alpha_{i,j}^{\\prime}}{\\lambda_{i,i+1}}, & \\textrm{if }i\\in[1:N-1],j\\in[0:i-1]\\\\\n", "\\alpha_{N,j}^{\\prime}, & \\textrm{if }i=N.\n", "\\end{cases}\n", "$$" ], "metadata": {} }, { "cell_type": "code", "source": [ "# Recover α_{i,j} from α′_{i,j}\n", "# -----------------------------\n", "\n", "α_opt = OffsetArray(zeros(size(α′_opt)), 1:N, 0:N-1)\n", "\n", "for i in 1:N-1\n", " for j in 0:i-1\n", " α_opt[i,j] = (α′_opt[i,j] / λ_opt_i_ip1[i])\n", " end\n", "end\n", "\n", "for j in 0:N-1\n", " α_opt[N,j] = α′_opt[N,j]\n", "end" ], "outputs": [], "execution_count": 25, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:32.249Z", "iopub.execute_input": "2021-04-23T13:27:32.263Z", "iopub.status.idle": "2021-04-23T13:27:32.319Z" } } }, { "cell_type": "code", "source": [], "outputs": [], "execution_count": null, "metadata": {} }, { "cell_type": "markdown", "source": [ "We now, recover $h_{i,j}$ from $\\alpha_{i,j}$ using the formula\n", "$$\n", "\\left(\\forall i\\in[1:N]\\right)\\;\\left(\\forall j\\in[0:i-1]\\right)\\quad h_{i,j}=\\begin{cases}\n", "\\alpha_{i,j}, & \\textrm{if }j=i-1\\\\\n", "\\alpha_{i,j}-\\alpha_{i-1,j}, & \\textrm{if }j\\in[0:i-1].\n", "\\end{cases}\n", "$$" ], "metadata": {} }, { "cell_type": "code", "source": [ "# Recover h_{i,j} from α_{i,j}\n", "h_opt = OffsetArray(zeros(size(α_opt)), 1:N, 0:N-1)\n", "\n", "# set h(1,0)\n", "h_opt[1,0] = α_opt[1,0]\n", "\n", "for i in 2:N\n", " for j in 0:i-1\n", " if j == i-1\n", " h_opt[i,j] = α_opt[i,j]\n", " else\n", " h_opt[i,j] = α_opt[i,j] - α_opt[i-1,j]\n", " end\n", " end\n", "end" ], "outputs": [], "execution_count": 26, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:32.341Z", "iopub.execute_input": "2021-04-23T13:27:32.355Z", "iopub.status.idle": "2021-04-23T13:27:32.374Z" } } }, { "cell_type": "markdown", "source": [ "**Putting everything in a function.** Now, let us put everything in a function. We need to add the packages first." ], "metadata": {} }, { "cell_type": "code", "source": [ "function full_pep_solver(N, L, μ, c_w, c_f)\n", "\n", " # define all the bold vectors\n", " # --------------------------\n", "\n", " # define 𝐰_0\n", "\n", " 𝐰_0 = e_i(N+2, 1)\n", "\n", " 𝐰_star = zeros(N+2, 1)\n", "\n", " 𝐠_star = zeros(N+2,1)\n", "\n", " # define 𝐠_0, 𝐠_1, …, 𝐠_N\n", "\n", " # first we define 𝐠_Julia vectors and then 𝐠 vectors\n", "\n", " 𝐠 = OffsetArray(zeros(N+2, N+1), 1:N+2, 0:N)\n", " # 𝐠= [𝐠_0 𝐠_1 𝐠_2 ... 𝐠_N]\n", " for i in 0:N\n", " 𝐠[:,i] = e_i(N+2, i+2)\n", " end\n", "\n", "\n", " 𝐟_star = zeros(N+1,1)\n", "\n", " # time to define 𝐟_0, 𝐟_1, …, 𝐟_N\n", "\n", " 𝐟 = OffsetArray(zeros(N+1, N+1), 1:N+1, 0:N)\n", " # 𝐟 = [𝐟_0, 𝐟_1, …, 𝐟_N]\n", "\n", " for i in 0:N\n", " 𝐟[:,i] = e_i(N+1, i+1)\n", " end\n", "\n", " # Define JuMP model\n", " # ----------------\n", "\n", " pep_model = Model(optimizer_with_attributes(Mosek.Optimizer, \"MSK_DPAR_INTPNT_CO_TOL_PFEAS\" => 1e-10))\n", "\n", " #define all the variables\n", " # -------------------------------\n", "\n", " # define α′ (can be typed as \\alpha[TAB]\\prime[TAB])\n", " @variable(pep_model, α′[1:N, 0:N-1])\n", "\n", " # define λ variables\n", "\n", " @variable(pep_model, λ_i_ip1[0:N-1] >= 0)\n", " # this defines (λ_{i,i+1})_{i∈[0:N-1]} in Julia indexing\n", "\n", " @variable(pep_model, λ_star_i[0:N] >= 0)\n", " # this defines (λ_{⋆,i})_{i∈[0:N]} in Julia indexing\n", "\n", " # define τ\n", " @variable(pep_model, τ >= 0)\n", "\n", " # define objective\n", " # ------------------\n", "\n", " @objective(\n", " pep_model,\n", " Min,\n", " τ\n", " )\n", "\n", " # Add the linear equality constraint\n", " # ----------------------------------\n", " @constraint(pep_model,\n", " τ * c_f .* 𝐟[:,0] + sum(λ_i_ip1[i] .* (𝐟[:,i+1]-𝐟[:,i]) for i in 0:N-1)\n", " + sum(λ_star_i[i] .* (𝐟[:,i] - 𝐟_star) for i in 0:N)\n", " .== 𝐟[:,N] - 𝐟_star\n", " )\n", "\n", " # Add the giant LMI constraint step by step\n", " # ----------------------------------------\n", "\n", " # Define all the terms one by one\n", "\n", " term_1 = c_w * τ * ⊙(𝐰_0,𝐰_0)\n", "\n", " term_2_part_1 = sum(λ_i_ip1[i] .* ⊙(𝐠[:,i]-𝐠[:,i+1],𝐠[:,i]-𝐠[:,i+1]) for i in 0:N-1)\n", "\n", " term_2_part_2 = sum(λ_star_i[i] .* ⊙(𝐠_star - 𝐠[:,i],𝐠_star - 𝐠[:,i]) for i in 0:N)\n", "\n", " term_2 = (term_2_part_1 + term_2_part_2) ./ (2*L)\n", "\n", " term_3 = - λ_star_i[0] .* ⊙(𝐠[:,0],𝐰_0)\n", "\n", " term_4_part_1 = - sum( λ_i_ip1[i] .* ⊙(𝐠[:,i],𝐰_0) for i in 1:N-1)\n", "\n", " term_4_part_2 = (1/L) .*\n", " sum( sum(α′[i,j] .* ⊙(𝐠[:,i],𝐠[:,j]) for j in 0:i-1) for i in 1:N-1)\n", "\n", " term_4 = term_4_part_1 + term_4_part_2\n", "\n", " term_5 = - ( ⊙(𝐠[:,N],𝐰_0) - (1/L) .* sum(α′[N,j] .* ⊙(𝐠[:,N],𝐠[:,j]) for j in 0:N-1))\n", "\n", " term_6_part_1 = sum(λ_i_ip1[i] .* ⊙(𝐠[:,i+1],𝐰_0) for i in 0:N-1)\n", "\n", " term_6_part_2 = - (1/L) .*\n", " sum( sum( α′[i,j] .* ⊙(𝐠[:,i+1],𝐠[:,j]) for j in 0:i-1) for i in 1:N-1)\n", "\n", " term_6 = term_6_part_1 + term_6_part_2\n", "\n", " # oof, okay constructed the terms for the LMI constraint, let us hope that there is no mistake and add them together\n", "\n", " S_mat = term_1 + term_2 + term_3 + term_4 + term_5 + term_6\n", "\n", " # add the LMI constraint now\n", "\n", " @constraint(pep_model,\n", " S_mat in PSDCone()\n", " )\n", "\n", " # time to optimize!\n", " # -----------------\n", " optimize!(pep_model)\n", "\n", " println(\"termination status =\", termination_status(pep_model) )\n", "\n", " α′_opt = OffsetArray(value.(α′), 1:N, 0:N-1)\n", " λ_opt_i_ip1 = OffsetVector(value.(λ_i_ip1), 0:N-1)\n", " λ_opt_star_i = OffsetVector(value.(λ_star_i), 0:N)\n", " τ_opt = value(τ)\n", "\n", " # Recover α_{i,j} from α′_{i,j}\n", " # -----------------------------\n", "\n", " α_opt = OffsetArray(zeros(size(α′_opt)), 1:N, 0:N-1)\n", "\n", " for i in 1:N-1\n", " for j in 0:i-1\n", " α_opt[i,j] = (α′_opt[i,j] / λ_opt_i_ip1[i])\n", " end\n", " end\n", "\n", " for j in 0:N-1\n", " α_opt[N,j] = α′_opt[N,j]\n", " end\n", "\n", " # Recover h_{i,j} from α_{i,j}\n", " # ----------------------------\n", "\n", " h_opt = OffsetArray(zeros(size(α_opt)), 1:N, 0:N-1)\n", "\n", " # set h(1,0)\n", " h_opt[1,0] = α_opt[1,0]\n", "\n", " # set the rest of the variables\n", "\n", " for i in 2:N\n", " for j in 0:i-1\n", " if j == i-1\n", " h_opt[i,j] = α_opt[i,j]\n", " else\n", " h_opt[i,j] = α_opt[i,j] - α_opt[i-1,j]\n", " end\n", " end\n", " end\n", "\n", " # return all the outputs\n", " # ----------------------\n", "\n", " return α′_opt, λ_opt_i_ip1, λ_opt_star_i, τ_opt, α_opt, h_opt\n", "\n", "end" ], "outputs": [ { "output_type": "execute_result", "execution_count": 27, "data": { "text/plain": "full_pep_solver (generic function with 1 method)" }, "metadata": {} } ], "execution_count": 27, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:32.394Z", "iopub.execute_input": "2021-04-23T13:27:32.408Z", "iopub.status.idle": "2021-04-23T13:27:32.490Z" } } }, { "cell_type": "markdown", "source": [ "Time to run and test." ], "metadata": {} }, { "cell_type": "code", "source": [ "α′_opt, λ_opt_i_ip1, λ_opt_star_i, τ_opt, α_opt, h_opt = full_pep_solver(N, L, μ, c_w, c_f)" ], "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 34 \n", " Cones : 0 \n", " Scalar variables : 37 \n", " Matrix variables : 1 \n", " Integer variables : 0 \n", "\n", "Optimizer started.\n", "Presolve started.\n", "Linear dependency checker started.\n", "Linear dependency checker terminated.\n", "Eliminator started.\n", "Freed constraints in eliminator : 5\n", "Eliminator terminated.\n", "Eliminator started.\n", "Freed constraints in eliminator : 0\n", "Eliminator terminated.\n", "Eliminator - tries : 2 time : 0.00 \n", "Lin. dep. - tries : 1 time : 0.00 \n", "Lin. dep. - number : 0 \n", "Presolve terminated. Time: 0.00 \n", "Problem\n", " Name : \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 34 \n", " Cones : 0 \n", " Scalar variables : 37 \n", " Matrix variables : 1 \n", " Integer variables : 0 \n", "\n", "Optimizer - threads : 4 \n", "Optimizer - solved problem : the primal \n", "Optimizer - Constraints : 29\n", "Optimizer - Cones : 1\n", "Optimizer - Scalar variables : 23 conic : 16 \n", "Optimizer - Semi-definite variables: 1 scalarized : 28 \n", "Factor - setup time : 0.00 dense det. time : 0.00 \n", "Factor - ML order time : 0.00 GP order time : 0.00 \n", "Factor - nonzeros before factor : 423 after factor : 435 \n", "Factor - dense dim. : 0 flops : 1.21e+04 \n", "ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n", "0 1.0e+00 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.00 \n", "1 2.7e-01 2.7e-01 1.2e-01 2.14e+00 -2.378627944e-02 -6.117274954e-03 2.7e-01 0.00 \n", "2 1.2e-01 1.2e-01 2.9e-02 1.32e+00 7.175923608e-02 7.139302023e-02 1.2e-01 0.00 \n", "3 4.2e-02 4.2e-02 6.9e-03 1.05e+00 4.477758503e-02 4.729110878e-02 4.2e-02 0.00 \n", "4 2.3e-02 2.3e-02 2.9e-03 7.37e-01 2.880293993e-02 3.057913998e-02 2.3e-02 0.00 \n", "5 6.7e-03 6.7e-03 4.0e-04 9.85e-01 3.356646000e-02 3.302278094e-02 6.7e-03 0.00 \n", "6 2.5e-03 2.5e-03 1.1e-04 6.92e-01 2.291866424e-02 2.288414204e-02 2.5e-03 0.00 \n", "7 6.6e-04 6.6e-04 1.5e-05 8.71e-01 2.066303721e-02 2.066465445e-02 6.6e-04 0.00 \n", "8 9.8e-05 9.8e-05 8.8e-07 8.83e-01 1.883068878e-02 1.882786667e-02 9.8e-05 0.00 \n", "9 1.7e-06 1.7e-06 1.9e-09 9.93e-01 1.859217137e-02 1.859212284e-02 1.7e-06 0.02 \n", "10 4.0e-08 4.0e-08 7.3e-12 1.00e+00 1.858820828e-02 1.858820710e-02 4.0e-08 0.02 \n", "11 2.1e-09 2.1e-09 8.5e-14 1.00e+00 1.858813924e-02 1.858813917e-02 2.1e-09 0.02 \n", "12 7.9e-11 9.5e-11 6.3e-16 1.00e+00 1.858813673e-02 1.858813673e-02 7.9e-11 0.02 \n", "Optimizer terminated. Time: 0.02 \n", "\n", "termination status =OPTIMAL\n" ] }, { "output_type": "execute_result", "execution_count": 28, "data": { "text/plain": "([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])" }, "metadata": {} } ], "execution_count": 28, "metadata": { "execution": { "iopub.status.busy": "2021-04-23T13:27:32.510Z", "iopub.execute_input": "2021-04-23T13:27:32.523Z", "iopub.status.idle": "2021-04-23T13:27:36.208Z" } } } ], "metadata": { "language_info": { "file_extension": ".jl", "name": "julia", "mimetype": "application/julia", "version": "1.6.0" }, "kernelspec": { "name": "julia-1.6", "display_name": "Julia 1.6.0", "language": "julia" }, "nteract": { "version": "0.28.0" } }, "nbformat": 4, "nbformat_minor": 2 }