NExOS.jl
¶Shuvomoy Das Gupta
The problem in consideration can be written as:
$$ \begin{equation} \begin{array}{ll} & \textrm{minimize} & \left\Vert \mathcal{A}(X)-b\right\Vert _{2}^{2}\\ & \textrm{subject to} & \mathop{{\bf rank}}(X)\leq r\\ & & \|X\|_{2}\leq M\\ & & X\in\mathbf{R}^{m\times n}, \end{array} \end{equation} $$where $X\in\mathbf{R}^{m\times n}$ is the decision variable, $b\in\mathbf{R}^{k}$ is a noisy measurement data, and $\mathcal{A}:\mathbf{R}^{m\times n}\to\mathbf{R}^{k}$ is a linear map. The affine map $\mathcal{A}$ can be determined by $k$ matrices $A_{1},\ldots,A_{k}$ in $\mathbf{R}^{m\times n}$ such that
$$ \mathcal{A}(X)=(\mathbf{tr}(A_{1}^{T}X),\ldots,\mathbf{tr}(A_{k}^{T}X)). $$First we generate our data for this example.
using Random, NExOS, ProximalOperators, LinearAlgebra
Random.seed!(1234)
m = 10
n = 2*m
M = 1.0
k = convert(Int64, m*n/20)
10
Here, $k$ is the number of components in output of the affine operator $\mathcal{A}$, i.e., for any matrix $X \in \mathbf{X}$, we have $\mathcal{A}(X) \in \mathbf{R}^k$.
r = convert(Int64,round(m*.35))
4
Here, $r$ corresponds to the rank of the matrix $\mathbf{rank}(X) <= r$.
The affine operator $\mathcal{A}$ is passed as a $k \times mn$ matrix $\bar{A}$, so that when it acts on $\mathbf{vec}(X)$ we have $\bar{A} (\mathbf{vec}(X)) = \mathcal{A}(X)$.
barA = randn(k, m*n)
10×200 Array{Float64,2}: 0.867347 -0.560501 1.56417 … 0.384485 0.685573 -0.858928 -0.901744 -0.0192918 -1.39674 0.315489 1.21533 1.12188 -0.494479 0.128064 1.1055 -0.382068 1.29772 -2.45236 -0.902914 1.85278 -1.10673 0.691941 -1.71088 -2.30555 0.864401 -0.827763 -3.21136 0.0473293 -0.747106 1.54823 2.21188 0.110096 -0.0740145 … -0.455901 0.0330671 -0.297855 0.532813 -0.251176 0.150976 0.100961 2.05177 1.58331 -0.271735 0.369714 0.769278 -1.12375 1.05237 0.562472 0.502334 0.0721164 -0.310153 -0.579068 0.430384 0.85342 -0.516984 -1.50343 -0.602707 -0.493044 0.211279 -0.321671
The vector $b$ is the observed output.
b = randn(k)
10-element Array{Float64,1}: -0.4667377936599136 -0.17047970356351902 0.8292034377961193 -0.4500585937344793 -1.3038856208344294 0.5869555339344937 0.17548586215288262 -0.2760027307659979 -0.2631151380019278 -1.1348769546238908
Let us chose an initial point.
Z0 = zeros(m,n);
Create the objective function.
f = LeastSquaresOverMatrix(barA, b, 1.0, iterative = true);
Let us create the constraint set.
C = RankSet(M, r)
RankSet{Float64,Int64}(1.0, 4)
Create the settings
file, note that we are working with the :safe
rule for updating the proximal parameter $\gamma$. Also, the solver output is turned off here, we could turn it on by settings verbose=true
.
settings = Settings(μ_max = 2, μ_mult_fact = 0.5, verbose = false, freq = 500, γ_updt_rule = :safe)
problem = Problem(f, C, settings.β, Z0)
Problem{ProximalOperators.LeastSquaresIterative{1,Float64,Float64,Array{Float64,2},Array{Float64,1},ProximalOperators.AAc},RankSet{Float64,Int64},Float64,Array{Float64,2}}(description : Least squares penalty domain : n/a expression : n/a parameters : n/a, RankSet{Float64,Int64}(1.0, 4), 1.0e-10, [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])
Solve the problem!
state_final = solve!(problem, settings)
State{Array{Float64,2},Int64,Float64}([0.009640555740470996 -0.0039049066980724172 … 0.0032744246264381385 -0.011401924226827814; 0.017039021787684776 -0.006486467044692794 … -0.006723985185855889 0.011346398802052203; … ; -0.0031155465177579336 0.012103228512837969 … -0.0030847566736254646 0.006789864137943311; 0.0037643322883237993 0.011064929502908967 … 0.003141509734574783 -0.006391901052935991], [0.009640555310037102 -0.003904906555286084 … 0.0032744247618116023 -0.011401924380723396; 0.017039022426592925 -0.006486466874719959 … -0.006723985197637846 0.011346398967171602; … ; -0.0031155462203965234 0.012103228542114463 … -0.0030847567074370006 0.006789863932342989; 0.0037643323913579916 0.011064929307709616 … 0.0031415097913574316 -0.00639190102063273], [0.009640555709940084 -0.003904906817854433 … 0.003274424252976659 -0.011401923965769611; 0.017039021755258916 -0.006486467029139976 … -0.006723985147301422 0.011346398638919817; … ; -0.0031155461667072546 0.012103228769906214 … -0.003084756747081626 0.006789864408033303; 0.0037643322644109536 0.011064929211436649 … 0.003141509432661576 -0.006391901114060031], 1, 9.267131556491698e-10, 2.192690473634684e-15, 7.450580596923828e-9, 0.002154434690031884)
Observe the output. Let us check the objective value first, we want it to be small.
f(state_final.x)
2.1913058045348726e-14
Finally, we want to check what if the best state found by our algorithm has converged to a locally optimal solution, which is checked by testing whether the best state has reached the desired fixed point gap and feasibility gap.
log10(state_final.fxd_pnt_gap) <= -4
log10(state_final.fsblt_gap) <= -4
true