TRDH algorithm

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
 NLPModelsModifiers = "e01155f1-5c6f-4375-a9d8-616dd036575f"
@@ -15,6 +16,7 @@ SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843"
 TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 
 [compat]
+LinearOperators = "2.5.2"
 NLPModels = "0.19"
 NLPModelsModifiers = "0.6"
 ProximalOperators = "0.15"

benchmarks/tables/bpdn-constr-table.jl

@@ -0,0 +1,29 @@
+include("regulopt-tables.jl")
+
+# model
+Random.seed!(1234)
+compound = 1
+model, nls_model, sol = bpdn_model(compound, bounds = true)
+
+f = LSR1Model(model)
+λ = norm(grad(model, zeros(model.meta.nvar)), Inf) / 10
+h = NormL1(λ)
+
+verbose = 0 # 10
+ϵ = 1.0e-6
+maxIter = 500
+maxIter_inner = 20
+options = ROSolverOptions(ν = 1.0, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options3 = ROSolverOptions(spectral = false, psb = false, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options6 = ROSolverOptions(ν = 1.0, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = false)
+options5 = ROSolverOptions(ν = 1.0, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true)
+options7 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner, reduce_TR = false)
+
+solvers = [:R2, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR]
+subsolvers = [:None, :None, :None, :None, :R2, :TRDH, :TRDH, :TRDH, :TRDH]
+solver_options = [options, options, options5, options6, options, options, options, options, options]
+subsolver_options = [options2, options2, options2, options2, options2, options7, options2, options3, options4] # n'importe lequel si subsolver = :None
+
+benchmark_table(f, 1:f.meta.nvar, sol, h, λ, solvers, subsolvers, solver_options, subsolver_options, "BPDN-cstr")

benchmarks/tables/bpdn-table.jl

@@ -0,0 +1,31 @@
+include("regulopt-tables.jl")
+
+# model
+Random.seed!(1234)
+compound = 1
+model, nls_model, sol = bpdn_model(compound, bounds = false)
+
+# parameters
+f = LSR1Model(model)
+λ = norm(grad(model, zeros(model.meta.nvar)), Inf) / 10
+h = NormL1(λ)
+
+verbose = 0 # 10
+ϵ = 1.0e-6
+maxIter = 500
+maxIter_inner = 40
+options = ROSolverOptions(ν = 1.0, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options3 = ROSolverOptions(spectral = false, psb = false, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options6 = ROSolverOptions(ν = 1.0, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = false)
+options5 = ROSolverOptions(ν = 1.0, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true)
+options7 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner, reduce_TR = false)
+options8 = ROSolverOptions(spectral = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner, reduce_TR = false)
+
+solvers = [:R2, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR, :TR]
+subsolvers = [:None, :None, :None, :None, :R2, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH]
+solver_options = [options, options, options5, options6, options, options, options, options, options, options]
+subsolver_options = [options2, options2, options2, options2, options2, options7, options2, options3, options4, options8] # n'importe lequel si subsolver = :None
+
+benchmark_table(f, 1:f.meta.nvar, sol, h, λ, solvers, subsolvers, solver_options, subsolver_options, "BPDN")

benchmarks/tables/fh-table.jl

@@ -0,0 +1,32 @@
+include("regulopt-tables.jl")
+using ADNLPModels, DifferentialEquations
+
+Random.seed!(1234)
+data, simulate, resid, misfit = RegularizedProblems.FH_smooth_term()
+model = ADNLPModel(misfit, ones(5))
+f = LBFGSModel(model)
+
+λ = 1.0
+h = NormL0(λ)
+ν = 1.0e2
+verbose = 0 #10
+maxIter = 1000
+maxIter_sub = 200 # max iter for subsolver
+ϵ = 1.0e-4
+options = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵ, ϵr = ϵ, maxIter = maxIter_sub)
+options3 = ROSolverOptions(spectral = false, psb = false, ϵa = ϵ, ϵr = ϵ, maxIter = maxIter_sub)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵ, ϵr = ϵ, maxIter = maxIter_sub)
+options5 = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true)
+options6 = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = false)
+options7 = ROSolverOptions(spectral = false, psb = true, reduce_TR = false, ϵa = ϵ, ϵr = ϵ, maxIter = maxIter_sub)
+options8 = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = false, reduce_TR = false)
+
+solvers = [:R2, :TRDH, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR]
+subsolvers = [:None, :None, :None, :None, :None, :R2, :TRDH, :TRDH, :TRDH, :TRDH]
+solver_options = [options, options, options5, options6, options8, options, options, options, options, options]
+subsolver_options = [options2, options2, options2, options2, options2, options2, options7, options2, options3, options4]
+subset = 2:10 # issues with R2 alone
+
+benchmark_table(f, 1:f.meta.nvar, [], h, λ, solvers[subset], subsolvers[subset], solver_options[subset], subsolver_options[subset],
+                "FH with ν = $ν, λ = $λ")

benchmarks/tables/nnmf-table.jl

@@ -0,0 +1,29 @@
+include("regulopt-tables.jl")
+
+Random.seed!(1234)
+m, n, k = 100, 50, 5
+model, A, selected = nnmf_model(m, n, k)
+f = LSR1Model(model)
+λ = norm(grad(model, rand(model.meta.nvar)), Inf) / 100
+h = NormL0(λ)
+ν = 1.0
+ϵ = 1.0e-5
+maxIter = 500
+maxIter_inner = 40
+verbose = 0 #10
+options = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options3 = ROSolverOptions(spectral = false, psb = false, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options5 = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true)
+options6 = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = false)
+options7 = ROSolverOptions(spectral = false, psb = true, reduce_TR = false, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options8 = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = false, reduce_TR = false)
+
+solvers = [:R2, :TRDH, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR, :TR]
+subsolvers = [:None, :None, :None, :None, :None, :R2, :TRDH, :TRDH, :TRDH, :TRDH, :TRDH]
+solver_options = [options, options, options5, options6, options8, options, options, options, options, options]
+subsolver_options = [options2, options2, options2, options2, options2, options2, options7, options2, options3, options4]
+
+benchmark_table(f, selected, [], h, λ, solvers, subsolvers, solver_options, subsolver_options,
+                "NNMF with m = $m, n = $n, k = $k, ν = $ν, λ = $λ")

benchmarks/tables/regulopt-tables.jl

@@ -0,0 +1,157 @@
+using PrettyTables
+using Random
+using LinearAlgebra
+using ProximalOperators
+using NLPModels, NLPModelsModifiers, RegularizedProblems, RegularizedOptimization, ShiftedProximalOperators, SolverBenchmark
+using Printf
+
+# utils for extracting stats / display table
+modelname(nlp::LSR1Model) = "LSR1"
+modelname(nlp::LBFGSModel) = "LBFGS"
+modelname(nlp::SpectralGradientModel) = "SpectralGradient"
+modelname(nlp::DiagonalQNModel) = "DiagonalQN"
+subsolvername(subsolver::Symbol) = subsolver == :None ? "" : string("-", subsolver)
+function options_str(options::ROSolverOptions, solver::Symbol, subsolver_options::ROSolverOptions, subsolver::Symbol)
+  if solver == :TRDH
+    out_str = !options.spectral ? (options.psb ? "-DiagQN-PSB" : "-DiagQN-Andrei") : "-Spectral"
+    out_str = (options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+  elseif solver == :TR && subsolver == :TRDH
+    out_str = !subsolver_options.spectral ? (subsolver_options.psb ? "-DiagQN-PSB" : "-DiagQN-Andrei") : "-Spectral"
+    out_str = (subsolver_options.reduce_TR) ? out_str : string(out_str, "-noredTR")
+  else
+    out_str = ""
+  end
+  return out_str
+end
+grad_evals(nlp::AbstractNLPModel) = neval_grad(nlp)
+grad_evals(nls::AbstractNLSModel) = neval_jtprod_residual(nls) + neval_jprod_residual(nls)
+function nb_prox_evals(stats, solver::Symbol)
+  if solver ∈ [:TR, :R2, :TRDH]
+    prox_evals = sum(stats.solver_specific[:SubsolverCounter])
+  else
+    error("not implemented")
+  end
+  return prox_evals
+end
+
+function benchmark_table(
+  f::AbstractNLPModel,
+  selected,
+  sol,
+  h,
+  λ,
+  solvers,
+  subsolvers,
+  solver_options,
+  subsolver_options,
+  pb_name::String,
+)
+
+  row_names = ["$(solver)$(subsolvername(subsolver))$(options_str(opt, solver, subsolver_opt, subsolver))" 
+                for (solver, opt, subsolver, subsolver_opt) in zip(solvers, solver_options, subsolvers, subsolver_options)]
+
+  n∇f_evals = []
+  nprox_evals = []
+  solver_stats = []
+
+  for (solver, subsolver, opt, sub_opt) in zip(solvers, subsolvers, solver_options, subsolver_options)
+    @info " using $solver with subsolver = $subsolver"
+    args = solver == :R2 ? () : (NormLinf(1.0),)
+    if subsolver == :None
+      solver_out = eval(solver)(f, h, args..., opt, x0 = f.meta.x0, selected = selected)
+    else
+      solver_out = eval(solver)(f, h, args..., opt, x0 = f.meta.x0, subsolver = eval(subsolver),
+                                subsolver_options = sub_opt, selected = selected)
+    end
+    push!(n∇f_evals, grad_evals(f))
+    push!(nprox_evals, nb_prox_evals(solver_out, solver))
+    push!(solver_stats, solver_out)
+    reset!(f)
+
+  end
+
+  if length(sol) == 0
+    header = [
+      "f(x)",
+      "h(x)/λ",
+      "ξ",
+      "∇f evals",
+      "prox calls",
+    ]
+  else
+    header = [
+      "f(x) (true = $(round(obj(model, sol); sigdigits = 4)))",
+      "h(x)/λ",
+      "ξ",
+      "||x-x*||/||x*||",
+      "∇f evals",
+      "prox calls",
+    ]
+  end
+
+  n_solvers = length(row_names)
+  data = Matrix{Any}(undef, n_solvers, length(header))
+  for i=1:n_solvers
+    solver_out = solver_stats[i]
+    x = solver_out.solution
+    fx = solver_out.solver_specific[:Fhist][end]
+    hx = solver_out.solver_specific[:Hhist][end]
+    ξ = solver_out.dual_feas
+    n∇f = n∇f_evals[i]
+    nprox = nprox_evals[i]
+    if length(sol) == 0
+      data[i, :] .= [fx, hx / λ, ξ, n∇f, nprox]
+    else
+      err = norm(x - sol) / norm(sol)
+      data[i, :] .= [fx, hx / λ, ξ, err, n∇f, nprox]
+    end
+  end
+
+  if length(sol) == 0
+    print_formats = ft_printf(["%7.3e", "%7.1e", "%7.1e", "%i","%i"], 1:length(header))
+  else
+    print_formats = ft_printf(["%7.3e", "%7.1e", "%7.1e", "%7.3e", "%i","%i"], 1:length(header))
+  end
+
+  return pretty_table(data; 
+      header = header,
+      row_names= row_names,
+      title = "$pb_name $(modelname(f)) $(typeof(h).name.name)",
+      # backend = Val(:latex),
+      formatters = (print_formats,
+        # (v, i, j) -> (SolverBenchmark.safe_latex_AbstractFloat(v)),
+        ),
+    )
+end
+
+
+
+# λ = norm(grad(model, rand(model.meta.nvar)), Inf) / 100000
+# h = NormL1(λ)
+# benchmark_table(f, selected, [], h, λ, solvers, subsolvers, solver_options, subsolver_options,
+#                 "NNMF with m = $m, n = $n, k = $k, ν = 1.0e-3,")
+
+
+
+
+# header = ["TR LSR1 L0Box", "R2 LSR1 L0Box", "LM L0Box", "LMTR L0Box"]
+# TR_out = TR(f, h, χ, options, x0 = f.meta.x0)
+# n∇f_TR = neval_grad(f)
+# prox_evals_TR = sum(TR_out.solver_specific[:SubsolverCounter])
+# reset!(f)
+# R2_out = R2(f, h, options, x0 = f.meta.x0)
+# n∇f_R2 = neval_grad(f)
+# prox_evals_R2 = R2_out.iter
+# reset!(f)
+# LM_out = LM(nls_model, h, options, x0 = nls_model.meta.x0)
+# n∇f_LM = neval_jtprod_residual(nls_model) + neval_jprod_residual(nls_model)
+# prox_evals_LM = sum(LM_out.solver_specific[:SubsolverCounter])
+# reset!(nls_model)
+# LMTR_out = LMTR(nls_model, h, χ, options, x0 = nls_model.meta.x0)
+# n∇f_LMTR = neval_jtprod_residual(nls_model) + neval_jprod_residual(nls_model)
+# prox_evals_LMTR = sum(LMTR_out.solver_specific[:SubsolverCounter])
+# reset!(nls_model)
+# n∇f_evals = [n∇f_TR, n∇f_R2, n∇f_LM, n∇f_LMTR]
+# nprox_evals = [prox_evals_TR, prox_evals_R2, prox_evals_LM, prox_evals_LMTR]
+
+# solver_stats = [TR_out, R2_out, LM_out, LMTR_out]

benchmarks/tables/svm-table.jl

@@ -0,0 +1,32 @@
+include("regulopt-tables.jl")
+using MLDatasets
+
+Random.seed!(1234)
+nlp_train, nls_train, sol_train  = RegularizedProblems.svm_train_model()
+nlp_test, nls_test, sol_test = RegularizedProblems.svm_test_model()
+f = LSR1Model(nlp_train)
+f_test = LSR1Model(nlp_test)
+λ = 10.0 #norm(grad(model, rand(model.meta.nvar)), Inf) / 10
+h = NormL1(λ)
+
+ν = 1.0e0
+verbose = 0 #10
+ϵ = 1.0e-6
+maxIter = 500
+maxIter_inner = 40
+options = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = true)
+options2 = ROSolverOptions(spectral = false, psb = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options3 = ROSolverOptions(spectral = false, psb = false, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options4 = ROSolverOptions(spectral = true, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+options5 = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = true)
+options6 = ROSolverOptions(ν = ν, ϵa = ϵ, ϵr = ϵ, verbose = verbose, maxIter = maxIter, spectral = false, psb = false)
+options7 = ROSolverOptions(spectral = false, psb = true, reduce_TR = false, ϵa = ϵ, ϵr = ϵ, maxIter=maxIter_inner)
+
+solvers = [:R2, :TRDH, :TRDH, :TRDH, :TR, :TR, :TR, :TR, :TR]
+subsolvers = [:None, :None, :None, :None, :R2, :TRDH, :TRDH, :TRDH, :TRDH]
+solver_options = [options, options, options5, options6, options, options, options, options, options]
+subsolver_options = [options2, options2, options2, options2, options2, options2, options7, options3, options4]
+subset = 1:9
+
+benchmark_table(f, 1:f.meta.nvar, [], h, λ, solvers[subset], subsolvers[subset], solver_options[subset], subsolver_options[subset],
+                "SVM with ν = $ν, λ = $λ")

examples/demo-bpdn-constr.jl

@@ -9,7 +9,7 @@ include("plot-utils-bpdn.jl")
 Random.seed!(1234)
 
 function demo_solver(f, nls, sol, h, χ, suffix = "l0-linf")
-  options = ROSolverOptions(ν = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10)
+  options = ROSolverOptions(ν = 1.0, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10, spectral = false, psb = true)
 
   @info " using TR to solve with" h χ
   reset!(f)
@@ -22,6 +22,12 @@ function demo_solver(f, nls, sol, h, χ, suffix = "l0-linf")
   R2_out = R2(f, h, options, x0 = f.meta.x0)
   @info "R2 relative error" norm(R2_out.solution - sol) / norm(sol)
   plot_bpdn(R2_out, sol, "constr-r2-$(suffix)")
+
+  @info " using TRDH to solve with" h χ
+  reset!(f)
+  TRDH_out = TRDH(f, h, χ, options, x0 = f.meta.x0)
+  @info "TRDH relative error" norm(TRDH_out.solution - sol) / norm(sol)
+  # plot_bpdn(TRDH_out, sol, "constr-trdh-$(suffix)")
 
   @info " using LMTR to solve with" h χ
   reset!(nls)

examples/demo-nnmf-constr.jl

@@ -9,16 +9,32 @@ include("plot-utils-nnmf.jl")
 Random.seed!(1234)
 
 function demo_solver(f, h, χ, selected, Avec, m, n, k, suffix = "l0-linf")
-  options = ROSolverOptions(ν = 1.0e-3, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10, maxIter = 500)
+  options = ROSolverOptions(ν = 1.0e-3, β = 1e16, ϵa = 1e-6, ϵr = 1e-6, verbose = 10, maxIter = 500, spectral = true)
   @info " using TR to solve with" h χ
   reset!(f)
   TR_out = TR(f, h, χ, options, selected = selected)
-  plot_nnmf(TR_out, Avec, m, n, k, "tr-r2-$suffix")
+  # plot_nnmf(TR_out, Avec, m, n, k, "tr-r2-$suffix")
 
   @info " using R2 to solve with" h
   reset!(f)
   R2_out = R2(f, h, options, selected = selected)
-  plot_nnmf(R2_out, Avec, m, n, k, "r2-$suffix")
+  # plot_nnmf(R2_out, Avec, m, n, k, "r2-$suffix")
+
+  @info " using TR with R2 as subproblem to solve with" h χ
+  reset!(f)
+  TR_out = TR(f, h, χ, options, selected = selected)
+  # plot_nnmf(TR_out, Avec, m, n, k, "tr-r2-$suffix")
+
+  subsolver_options = ROSolverOptions(spectral = false, psb = true)
+  @info " using TR with TRDH as subproblem to solve with" h χ
+  reset!(f)
+  TR2_out = TR(f, h, χ, options, selected = selected, subsolver = TRDH, subsolver_options = subsolver_options)
+  # plot_nnmf(TR2_out, Avec, m, n, k, "tr-trdh-$suffix")
+
+  @info " using TRDH to solve with" h χ
+  reset!(f)
+  TRDH_out = TRDH(f, h, χ, options, selected = selected)
+  # plot_nnmf(TRDH_out, Avec, m, n, k, "trdh-$suffix")
 end
 
 function demo_nnmf()