define RegularizedOptimizationModel type

A type representing the problem as a whole is necessary in order
to run benchmarks with SolverBenchmark.jl.

Author: dpo

Project.toml

@@ -19,14 +19,15 @@ MLDatasets = "^0.7.4"
 NLPModels = "0.16, 0.17, 0.18, 0.19, 0.20"
 Noise = "0.2"
 Requires = "1"
-julia = "^1.3.0"
+julia = "^1.6.0"
 
 [extras]
 ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 MLDatasets = "eb30cadb-4394-5ae3-aed4-317e484a6458"
+ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537"
 QuadraticModels = "f468eda6-eac5-11e8-05a5-ff9e497bcd19"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ADNLPModels", "DifferentialEquations", "MLDatasets", "QuadraticModels", "Test"]
+test = ["ADNLPModels", "DifferentialEquations", "MLDatasets", "ProximalOperators", "QuadraticModels", "Test"]

src/RegularizedProblems.jl

@@ -14,6 +14,12 @@ include("group_lasso_model.jl")
 include("nnmf.jl")
 
 function __init__()
+  @require ProximalOperators = "a725b495-10eb-56fe-b38b-717eba820537" begin
+    include("testset_bpdn.jl")
+    include("testset_lrcomp.jl")
+    include("testset_matrand.jl")
+    include("testset_group_lasso.jl")
+  end
   @require ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a" begin
     @require DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" begin
       include("fh_model.jl")

src/group_lasso_model.jl

@@ -1,10 +1,22 @@
 export group_lasso_model
 
-function group_lasso_data(m::Int, n::Int, g::Int, ag::Int, noise::Float64 = 0.01)
-  (m ≤ n) || error("number of rows ($m) should be ≤ number of columns ($n)")
-  (mod(n, g) == 0) || error("number of groups ($g) must divide evenly into number of rows ($n)")
-  (ag ≤ g) || error("number of active groups ($ag) must be smaller than the number of groups ($g)")
-
+function group_lasso_data(;
+  m::Int = 200,
+  n::Int = 512,
+  g::Int = 16,
+  ag::Int = 5,
+  noise::Float64 = 0.01,
+  compound::Int = 1,
+)
+  m ≤ n || error("number of rows ($m) should be ≤ number of columns ($n)")
+  mod(n, g) == 0 || error("number of groups ($g) must divide evenly into number of rows ($n)")
+  ag ≤ g || error("number of active groups ($ag) must be smaller than the number of groups ($g)")
+  compound > 0 || error("compound factor must be positive")
+
+  m = compound * m
+  n = compound * n
+  g = compound * g
+  ag = compound * ag
   x0 = zeros(n)
   active_groups = sort(randperm(g)[1:ag]) # pick out active groups
   group_eles = Int(n / g) # get number of elements in a group
@@ -25,12 +37,8 @@ function group_lasso_data(m::Int, n::Int, g::Int, ag::Int, noise::Float64 = 0.01
   A, b, b0, x0, g, active_groups, indset
 end
 
-group_lasso_data(compound::Int = 1, args...) =
-  group_lasso_data(200 * compound, 512 * compound, 16 * compound, 5 * compound, args...)
-
 """
-    model, nls_model, sol = group_lasso_model(args...)
-    model, nls_model, sol = group_lasso_model(compound = 1, args...)
+    model, nls_model, sol = group_lasso_model(; kwargs...)
 
 Return an instance of an `NLPModel` and `NLSModel` representing the group-lasso
 problem, i.e., the under-determined linear least-squares objective
@@ -42,28 +50,23 @@ vector following a normal distribution with mean zero and standard deviation σ.
 Note that with this format, all groups have a the same number of elements and the number of
 groups divides evenly into the total number of elements.
 
-## Arguments
-
-* `m :: Int`: the number of rows of A
-* `n :: Int`: the number of columns of A (with `n` ≥ `m`)
-* `g :: Int : the number of groups`
-* `ag :: Array{Int}`:  group-index denoting which groups are active (with `max(ag) ≤ g`), i.e. `[1, 4, 5]` when there are 7 groups
-* `noise :: Float64`: noise amount ϵ (default: 0.01).
-
-The second form calls the first form with arguments
+## Keyword Arguments
 
-    m = 200 * compound
-    n = 512 * compound
-    k =  10 * compound
+* `m :: Int`: the number of rows of A (default: 200)
+* `n :: Int`: the number of columns of A, with `n` ≥ `m` (default: 512)
+* `g :: Int`: the number of groups (default: 16)
+* `ag :: Int`: the number of active groups (default: 5)
+* `noise :: Float64`: noise amount (default: 0.01)
+* `compound :: Int`: multiplier for `m`, `n`, `g`, and `ag` (default: 1).
 
 ## Return Value
 
 An instance of a `FirstOrderModel` that represents the group-lasso problem.
 An instance of a `FirstOrderNLSModel` that represents the group-lasso problem.
-Also returns true x, number of groups g, group-index denoting which groups are active, and a Matrix where rows are group indices of x
+Also returns true x, number of groups g, group-index denoting which groups are active, and a Matrix where rows are group indices of x.
 """
-function group_lasso_model(args...)
-  A, b, b0, x0, g, active_groups, indset = group_lasso_data(args...)
+function group_lasso_model(args...; kwargs...)
+  A, b, b0, x0, g, active_groups, indset = group_lasso_data(args...; kwargs...)
   r = similar(b)
 
   function resid!(r, x)

src/lrcomp_model.jl

@@ -5,7 +5,7 @@ function lrcomp_data(m::Int, n::Int; T::DataType = Float64)
   A
 end
 
-function lrcomp_model(m::Int, n::Int; T::DataType = Float64)
+function lrcomp_model(; m::Int = 100, n::Int = 100, T::DataType = Float64)
   A = lrcomp_data(m, n, T = T)
   r = vec(similar(A))
 

src/matrand_model.jl

@@ -7,7 +7,7 @@ function mat_rand(m::Int, n::Int, r::Int, sr::Float64, va::Float64, vb::Float64,
   Ω = findall(<(sr), rand(m, n))
   B = xs[Ω]
   B = (1 - c) * add_gauss(B, va, 0; clip = true) + c * add_gauss(B, vb, 0; clip = true)
-  ω = zeros(Int64, size(Ω, 1))   # Vectorize Omega 
+  ω = zeros(Int64, size(Ω, 1))   # Vectorize Omega
   for i = 1:size(Ω, 1)
     ω[i] = Ω[i][1] + size(Ω, 2) * (Ω[i][2] - 1)
   end
@@ -44,7 +44,7 @@ function matrix_completion_model(xs, B, ω)
 end
 
 """
-    model, nls_model, sol = random_matrix_completion_model(args...)
+    model, nls_model, sol = random_matrix_completion_model(; kwargs...)
 
 Return an instance of an `NLPModel` and an instance of an `NLSModel` representing
 the same matrix completion problem, i.e., the square linear least-squares objective
@@ -55,38 +55,38 @@ in the Frobenius norm, where X is the unknown image represented as an m x n matr
 A is a fixed image, and the operator P only retains a certain subset of pixels of
 X and A.
 
-## Arguments
+## Keyword Arguments
 
-* `m :: Int`: the number of rows of X and A
-* `n :: Int`: the number of columns of X and A
-* `r :: Int`: the desired rank of A
+* `m :: Int`: the number of rows of X and A (default: 100)
+* `n :: Int`: the number of columns of X and A (default: 100)
+* `r :: Int`: the desired rank of A (default: 5)
 * `sr :: AbstractFloat`: a threshold between 0 and 1 used to determine the set of pixels
-  retained by the operator P
-* `va :: AbstractFloat`: the variance of a first Gaussian perturbation to be applied to A
-* `vb :: AbstractFloat`: the variance of a second Gaussian perturbation to be applied to A
-* `c :: AbstractFloat`: the coefficient of the convex combination of the two Gaussian perturbations.
+retained by the operator P (default: 0.8)
+* `va :: AbstractFloat`: the variance of a first Gaussian perturbation to be applied to A (default: 1.0e-4)
+* `vb :: AbstractFloat`: the variance of a second Gaussian perturbation to be applied to A (default: 1.0e-2)
+* `c :: AbstractFloat`: the coefficient of the convex combination of the two Gaussian perturbations (default: 0.2).
 
 ## Return Value
 
 An instance of a `FirstOrderModel` and of a `FirstOrderNLSModel` that represent the same
 matrix completion problem, and the exact solution.
 """
-function random_matrix_completion_model(
-  m::Int,
-  n::Int,
-  r::Int,
-  sr::R,
-  va::R,
-  vb::R,
-  c::R,
-) where {R <: AbstractFloat}
+function random_matrix_completion_model(;
+  m::Int = 100,
+  n::Int = 100,
+  r::Int = 5,
+  sr::Float64 = 0.8,
+  va::Float64 = 1.0e-4,
+  vb::Float64 = 1.0e-2,
+  c::Float64 = 0.2,
+)
   xs, B, ω = mat_rand(m, n, r, sr, va, vb, c)
   matrix_completion_model(xs, B, ω)
 end
 
 function perturb(I, c = 0.8, p = 0.8)
   Ω = findall(<(p), rand(256, 256))
-  ω = zeros(Int, size(Ω, 1))   # Vectorize Omega 
+  ω = zeros(Int, size(Ω, 1))   # Vectorize Omega
   for i = 1:size(Ω, 1)
     ω[i] = Ω[i][1] + 256 * (Ω[i][2] - 1)
   end
@@ -98,7 +98,7 @@ function perturb(I, c = 0.8, p = 0.8)
 end
 
 """
-    model, nls_model, sol = MIT_matrix_completion_model(args...)
+    model, nls_model, sol = MIT_matrix_completion_model()
 
 A special case of matrix completion problem in which the exact image is a noisy
 MIT logo.

src/testset_bpdn.jl

@@ -0,0 +1,22 @@
+# Predefine a set of common problem instances.
+export setup_bpdn_l0, setup_bpdn_l1, setup_bpdn_B0
+
+function setup_bpdn_l0(args...; kwargs...)
+  model, nls_model, _ = bpdn_model(args...)
+  λ = norm(grad(model, zeros(model.meta.nvar)), Inf) / 10
+  h = NormL0(λ)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end
+
+function setup_bpdn_l1(args...; kwargs...)
+  model, nls_model, _ = bpdn_model(args...)
+  λ = norm(grad(model, zeros(model.meta.nvar)), Inf) / 10
+  h = NormL1(λ)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end
+
+function setup_bpdn_B0(compound = 1, args...; kwargs...)
+  model, nls_model, _ = bpdn_model(compound, args...)
+  h = IndBallL0(10 * compound)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end

src/testset_group_lasso.jl

@@ -0,0 +1,11 @@
+# Predefine a set of common problem instances.
+export setup_group_lasso_l12
+
+function setup_group_lasso_l12(args...; kwargs...)
+  model, nls_model, ng, _, idx = group_lasso_model(; kwargs...)
+  idx = [idx[i, :] for i = 1:ng]
+  λ = 0.2 * ones(ng)
+  h = GroupNormL2(λ, idx)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end
+

src/testset_lrcomp.jl

@@ -0,0 +1,16 @@
+# Predefine a set of common problem instances.
+export setup_lrcomp_rank, setup_lrcomp_nuclear
+
+function setup_lrcomp_rank(args...; kwargs...)
+  model, nls_model, _ = lrcomp_model(args...; kwargs...)
+  λ = 0.1
+  h = Rank(λ)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end
+
+function setup_lrcomp_nuclear(args...; kwargs...)
+  model, nls_model, _ = lrcomp_model(args...; kwargs...)
+  λ = 0.1
+  h = NuclearNorm(λ)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end

src/testset_matrand.jl

@@ -0,0 +1,31 @@
+# Predefine a set of common problem instances.
+export setup_random_completion_rank, setup_random_completion_nuclear
+export setup_mit_completion_rank, setup_mit_completion_nuclear
+
+function setup_random_completion_rank(args...; kwargs...)
+  model, nls_model, _ = random_matrix_completion_model(; kwargs...)
+  λ = 0.1
+  h = Rank(λ)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end
+
+function setup_random_completion_nuclear(args...; kwargs...)
+  model, nls_model, _ = random_matrix_completion_model(; kwargs...)
+  λ = 0.1
+  h = NuclearNorm(λ)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end
+
+function setup_mit_completion_rank(args...; kwargs...)
+  model, nls_model, _ = MIT_matrix_completion_model()
+  λ = 0.1
+  h = Rank(λ)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end
+
+function setup_mit_completion_nuclear(args...; kwargs...)
+  model, nls_model, _ = MIT_matrix_completion_model()
+  λ = 0.1
+  h = NuclearNorm(λ)
+  return RegularizedNLPModel(model, h), RegularizedNLSModel(nls_model, h)
+end

src/types.jl

@@ -1,4 +1,5 @@
-export FirstOrderModel, FirstOrderNLSModel
+export FirstOrderModel,
+  FirstOrderNLSModel, AbstractRegularizedNLPModel, RegularizedNLPModel, RegularizedNLSModel
 
 """
     model = FirstOrderModel(f, ∇f!; name = "first-order model")
@@ -131,3 +132,98 @@ function NLPModels.jtprod_residual!(
   nls.jtprod_resid!(Jtv, x, v)
   Jtv
 end
+
+abstract type AbstractRegularizedNLPModel{T, S} <: AbstractNLPModel{T, S} end
+
+"""
+    rmodel = RegularizedNLPModel(model, regularizer)
+    rmodel = RegularizedNLSModel(model, regularizer)
+
+An aggregate type to represent a regularized optimization model, .i.e.,
+of the form
+
+    minimize f(x) + h(x),
+
+where f is smooth (and is usually assumed to have Lipschitz-continuous gradient),
+and h is lower semi-continuous (and may have to be prox-bounded).
+
+The regularized model is made of
+
+- `model <: AbstractNLPModel`: the smooth part of the model, for example a `FirstOrderModel`
+- `h`: the nonsmooth part of the model; typically a regularizer defined in `ProximalOperators.jl`
+- `selected`: the subset of variables to which the regularizer h should be applied (default: all).
+
+This aggregate type can be used to call solvers with a single object representing the
+model, but is especially useful for use with SolverBenchmark.jl, which expects problems
+to be defined by a single object.
+"""
+mutable struct RegularizedNLPModel{T, S, M <: AbstractNLPModel{T, S}, H, I} <:
+               AbstractRegularizedNLPModel{T, S}
+  model::M     # smooth  model
+  h::H         # regularizer
+  selected::I  # set of variables to which the regularizer should be applied
+end
+
+function RegularizedNLPModel(model::AbstractNLPModel{T, S}, h::H) where {T, S, H}
+  selected = 1:get_nvar(model)
+  RegularizedNLPModel{T, S, typeof(model), typeof(h), typeof(selected)}(model, h, selected)
+end
+
+mutable struct RegularizedNLSModel{T, S, M <: AbstractNLSModel{T, S}, H, I} <:
+               AbstractRegularizedNLPModel{T, S}
+  model::M     # smooth  model
+  h::H         # regularizer
+  selected::I  # set of variables to which the regularizer should be applied
+end
+
+function RegularizedNLSModel(model::AbstractNLSModel{T, S}, h::H) where {T, S, H}
+  selected = 1:get_nvar(model)
+  RegularizedNLSModel{T, S, typeof(model), typeof(h), typeof(selected)}(model, h, selected)
+end
+
+function NLPModels.obj(rnlp::AbstractRegularizedNLPModel, x::AbstractVector)
+  # The size check on x will be performed when evaluating the smooth model.
+  # We intentionally do not increment an objective evaluation counter here
+  # because the relevant counters are inside the smooth term.
+  obj(rnlp.model, x) + rnlp.h(x)
+end
+
+# Forward meta getters so they grab info from the smooth model
+for field ∈ fieldnames(NLPModels.NLPModelMeta)
+  meth = Symbol("get_", field)
+  if field == :name
+    @eval NLPModels.$meth(rnlp::RegularizedNLPModel) =
+      NLPModels.$meth(rnlp.model) * "/" * string(typeof(rnlp.h).name.wrapper)
+    @eval NLPModels.$meth(rnls::RegularizedNLSModel) =
+      NLPModels.$meth(rnls.model) * "/" * string(typeof(rnls.h).name.wrapper)
+  else
+    @eval NLPModels.$meth(rnlp::RegularizedNLPModel) = NLPModels.$meth(rnlp.model)
+  end
+end
+
+for field in fieldnames(NLPModels.NLSMeta)
+  meth = Symbol("get_", field)
+  @eval NLPModels.$meth(rnls::RegularizedNLSModel) = NLPModels.$meth(rnls.model)
+end
+
+# Forward counter getters so they grab info from the smooth model
+for model_type ∈ (RegularizedNLPModel, RegularizedNLSModel)
+  for counter in fieldnames(Counters)
+    @eval NLPModels.$counter(rnlp::$model_type) = NLPModels.$counter(rnlp.model)
+  end
+end
+
+for counter in fieldnames(NLSCounters)
+  counter == :counters && continue
+  @eval NLPModels.$counter(rnls::RegularizedNLSModel) = NLPModels.$counter(rnls.model)
+end
+
+# simple show method for now
+function Base.show(io::IO, rnlp::AbstractRegularizedNLPModel)
+  print(io, "Smooth model: ")
+  show(io, rnlp.model)
+  print(io, "\nRegularizer: ")
+  show(io, rnlp.h)
+  print(io, "\n\nSelected variables: ")
+  show(io, rnlp.selected)
+end