Quadratic solve structure

Project.toml

@@ -59,4 +59,4 @@ Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea"
 
 [targets]
-test = ["CSV", "Combinatorics", "DataFrames", "Distributions", "DoubleFloats", "DynamicPolynomials", "FiniteDifferences", "ForwardDiff", "GLPK", "JSON", "JuMP", "LaTeXStrings", "MAT", "MultivariatePolynomials", "Plots", "PlotThemes", "Polyhedra", "ReverseDiff", "ZipFile", "Test", "Tullio", "Clp", "HiGHS", "Hypatia"]
+test = ["CSV", "Combinatorics", "DataFrames", "Distributions", "DoubleFloats", "DynamicPolynomials", "FiniteDifferences", "ForwardDiff", "GLPK", "HiGHS", "JSON", "JuMP", "LaTeXStrings", "MAT", "MultivariatePolynomials", "Plots", "PlotThemes", "Polyhedra", "ReverseDiff", "ZipFile", "Test", "Tullio", "Clp", "Hypatia"]

examples/blended_pairwise_with_direct.jl

@@ -16,33 +16,21 @@ using FrankWolfe
 using LinearAlgebra
 using Random
 
-import Pkg
-Pkg.add("GLPK")
-Pkg.add("HiGHS")
-import GLPK
 import HiGHS
-
-# lp_solver = GLPK.Optimizer
-lp_solver = HiGHS.Optimizer
-# lp_solver = Clp.Optimizer # buggy / does not work properly
-
+import MathOptInterface as MOI
 
 include("../examples/plot_utils.jl")
 
 n = Int(1e4)
-k = 10000
+k = 10_000
 
-# s = rand(1:100)
 s = 10
 @info "Seed $s"
 Random.seed!(s)
 
 xpi = rand(n);
 total = sum(xpi);
 
-# here the optimal solution lies in the interior if you want an optimal solution on a face and not the interior use:
-# const xp = xpi;
-
 const xp = xpi ./ total;
 
 f(x) = norm(x - xp)^2
@@ -52,13 +40,7 @@ end
 
 lmo = FrankWolfe.KSparseLMO(5, 1.0)
 
-## other LMOs to try
-# lmo_big = FrankWolfe.KSparseLMO(100, big"1.0")
-# lmo = FrankWolfe.LpNormLMO{Float64,5}(1.0)
-# lmo = FrankWolfe.ProbabilitySimplexOracle(1.0);
-# lmo = FrankWolfe.UnitSimplexOracle(1.0);
-
-x00 = FrankWolfe.compute_extreme_point(lmo, rand(n))
+const x00 = FrankWolfe.compute_extreme_point(lmo, rand(n))
 
 function build_callback(trajectory_arr)
     return function callback(state, active_set, args...)
@@ -67,75 +49,112 @@ function build_callback(trajectory_arr)
 end
 
 
-FrankWolfe.benchmark_oracles(f, grad!, () -> randn(n), lmo; k=100)
-
 trajectoryBPCG_standard = []
-callback = build_callback(trajectoryBPCG_standard)
-
-x0 = deepcopy(x00)
 @time x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
     f,
     grad!,
     lmo,
-    x0,
+    copy(x00),
     max_iteration=k,
     line_search=FrankWolfe.Shortstep(2.0),
-    print_iter=k / 10,
-    memory_mode=FrankWolfe.InplaceEmphasis(),
     verbose=true,
-    trajectory=true,
-    callback=callback,
-    squadratic=false,
+    callback=build_callback(trajectoryBPCG_standard),
 );
 
+# Just projection quadratic
 trajectoryBPCG_quadratic = []
-callback = build_callback(trajectoryBPCG_quadratic)
+as_quad = FrankWolfe.ActiveSetQuadratic([(1.0, copy(x00))], 2 * LinearAlgebra.I, -2xp)
+@time x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+    f,
+    grad!,
+    lmo,
+    as_quad,
+    max_iteration=k,
+    line_search=FrankWolfe.Shortstep(2.0),
+    verbose=true,
+    callback=build_callback(trajectoryBPCG_quadratic),
+);
+
+as_quad = FrankWolfe.ActiveSetQuadratic(
+    [(1.0, copy(x00))],
+    2 * LinearAlgebra.I, -2xp,
+)
 
-x0 = deepcopy(x00)
+# with quadratic active set
+trajectoryBPCG_quadratic_as = []
 @time x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
     f,
     grad!,
     lmo,
-    x0,
+    as_quad,
     max_iteration=k,
     line_search=FrankWolfe.Shortstep(2.0),
-    print_iter=k / 10,
-    memory_mode=FrankWolfe.InplaceEmphasis(),
     verbose=true,
-    trajectory=true,
-    callback=callback,
-    squadratic=true,
-    lp_solver=lp_solver,
+    callback=build_callback(trajectoryBPCG_quadratic_as),
 );
 
-# Just using that qudratic (more expensive)
-trajectoryBPCG_quadratic_nosquad = []
-callback = build_callback(trajectoryBPCG_quadratic_nosquad)
+as_quad_direct = FrankWolfe.ActiveSetQuadraticLinearSolve(
+    [(1.0, copy(x00))],
+    2 * LinearAlgebra.I, -2xp,
+    MOI.instantiate(MOI.OptimizerWithAttributes(HiGHS.Optimizer, MOI.Silent() => true)),
+)
 
-x0 = deepcopy(x00)
+# with LP acceleration
+trajectoryBPCG_quadratic_direct = []
 @time x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
     f,
     grad!,
     lmo,
-    x0,
+    as_quad_direct,
     max_iteration=k,
     line_search=FrankWolfe.Shortstep(2.0),
-    print_iter=k / 10,
-    memory_mode=FrankWolfe.InplaceEmphasis(),
     verbose=true,
-    trajectory=true,
-    callback=callback,
-    quadratic=true,
-    squadratic=false,
-    lp_solver=lp_solver,
+    callback=build_callback(trajectoryBPCG_quadratic_direct),
 );
 
-# Update the data and labels for plotting
-dataSparsity = [trajectoryBPCG_standard, trajectoryBPCG_quadratic, trajectoryBPCG_quadratic_nosquad]
-labelSparsity = ["BPCG (Standard)", "BPCG (Specific Direct)", "BPCG (Generic Direct)"]
+as_quad_direct_generic = FrankWolfe.ActiveSetQuadraticLinearSolve(
+    [(1.0, copy(x00))],
+    2 * Diagonal(ones(length(xp))), -2xp,
+    MOI.instantiate(MOI.OptimizerWithAttributes(HiGHS.Optimizer, MOI.Silent() => true)),
+)
+
+# with LP acceleration
+trajectoryBPCG_quadratic_direct_generic = []
+@time x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+    f,
+    grad!,
+    lmo,
+    as_quad_direct_generic,
+    max_iteration=k,
+    line_search=FrankWolfe.Shortstep(2.0),
+    verbose=true,
+    callback=build_callback(trajectoryBPCG_quadratic_direct_generic),
+);
+
+as_quad_direct_basic_as = FrankWolfe.ActiveSetQuadraticLinearSolve(
+    FrankWolfe.ActiveSet([1.0], [copy(x00)], collect(x00)),
+    2 * LinearAlgebra.I, -2xp,
+    MOI.instantiate(MOI.OptimizerWithAttributes(HiGHS.Optimizer, MOI.Silent() => true)),
+)
+
+# with LP acceleration
+trajectoryBPCG_quadratic_noqas = []
 
-# Plot sparsity
-# plot_sparsity(dataSparsity, labelSparsity, legend_position=:topright)
+@time x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+    f,
+    grad!,
+    lmo,
+    as_quad_direct_basic_as,
+    max_iteration=k,
+    line_search=FrankWolfe.Shortstep(2.0),
+    verbose=true,
+    callback=build_callback(trajectoryBPCG_quadratic_noqas),
+);
+
+
+# Update the data and labels for plotting
+data_trajectories = [trajectoryBPCG_standard, trajectoryBPCG_quadratic, trajectoryBPCG_quadratic_as, trajectoryBPCG_quadratic_direct, trajectoryBPCG_quadratic_direct_generic, trajectoryBPCG_quadratic_noqas]
+labels_trajectories = ["BPCG (Standard)", "BPCG (Specific Direct)", "AS_Quad", "Reloaded", "Reloaded_generic", "Reloaded_noqas"]
 
 # Plot trajectories
-plot_trajectories(dataSparsity, labelSparsity, xscalelog=false)
+plot_trajectories(data_trajectories, labels_trajectories, xscalelog=false)

examples/blended_pairwise_with_direct_non-standard-quadratic.jl

@@ -14,16 +14,8 @@ using FrankWolfe
 using LinearAlgebra
 using Random
 
-import Pkg
-Pkg.add("GLPK")
-Pkg.add("HiGHS")
-import GLPK
 import HiGHS
-
-# lp_solver = GLPK.Optimizer
-lp_solver = HiGHS.Optimizer
-# lp_solver = Clp.Optimizer # buggy / does not work properly
-
+import MathOptInterface as MOI
 
 include("../examples/plot_utils.jl")
 
@@ -39,7 +31,6 @@ A = let
     A = randn(n, n)
     A' * A
 end
-
 @assert isposdef(A) == true
 
 const y = Random.rand(Bool, n) * 0.6 .+ 0.3
@@ -73,17 +64,14 @@ function build_callback(trajectory_arr)
 end
 
 
-FrankWolfe.benchmark_oracles(f, grad!, () -> randn(n), lmo; k=100)
-
 trajectoryBPCG_standard = []
 callback = build_callback(trajectoryBPCG_standard)
 
-x0 = deepcopy(x00)
-@time x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
     f,
     grad!,
     lmo,
-    x0,
+    copy(x00),
     max_iteration=k,
     line_search=FrankWolfe.Adaptive(),
     print_iter=k / 10,
@@ -94,33 +82,64 @@ x0 = deepcopy(x00)
     squadratic=false,
 );
 
-trajectoryBPCG_quadratic = []
-callback = build_callback(trajectoryBPCG_quadratic)
+active_set_quadratic_automatic = FrankWolfe.ActiveSetQuadraticLinearSolve(
+    [(1.0, copy(x00))],
+    grad!,
+    MOI.instantiate(MOI.OptimizerWithAttributes(HiGHS.Optimizer, MOI.Silent() => true)),
+    scheduler=FrankWolfe.LogScheduler(start_time=100, scaling_factor=1.2, max_interval=100),
+)
+trajectoryBPCG_quadratic_automatic = []
+x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+    f,
+    grad!,
+    lmo,
+    active_set_quadratic_automatic,
+    # print_iter=1,
+    max_iteration=k,
+    verbose=true,
+    callback=build_callback(trajectoryBPCG_quadratic_automatic),
+);
 
-x0 = deepcopy(x00)
-@time x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+active_set_quadratic_automatic2 = FrankWolfe.ActiveSetQuadraticLinearSolve(
+    [(1.0, copy(x00))],
+    grad!,
+    MOI.instantiate(MOI.OptimizerWithAttributes(HiGHS.Optimizer, MOI.Silent() => true)),
+    scheduler=FrankWolfe.LogScheduler(start_time=10, scaling_factor=2),
+)
+trajectoryBPCG_quadratic_automatic2 = []
+x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
     f,
     grad!,
     lmo,
-    x0,
+    active_set_quadratic_automatic2,
+    # print_iter=1,
     max_iteration=k,
-    line_search=FrankWolfe.Adaptive(),
-    print_iter=k / 10,
-    memory_mode=FrankWolfe.InplaceEmphasis(),
     verbose=true,
-    trajectory=true,
-    callback=callback,
-    quadratic=true,
-    lp_solver=lp_solver,
+    callback=build_callback(trajectoryBPCG_quadratic_automatic2),
+);
+
+
+active_set_quadratic_automatic_standard = FrankWolfe.ActiveSetQuadraticLinearSolve(
+    FrankWolfe.ActiveSet([(1.0, copy(x00))]),
+    grad!,
+    MOI.instantiate(MOI.OptimizerWithAttributes(HiGHS.Optimizer, MOI.Silent() => true)),
+    scheduler=FrankWolfe.LogScheduler(start_time=10, scaling_factor=2),
+)
+trajectoryBPCG_quadratic_automatic_standard = []
+x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+    f,
+    grad!,
+    lmo,
+    active_set_quadratic_automatic_standard,
+    max_iteration=k,
+    verbose=true,
+    callback=build_callback(trajectoryBPCG_quadratic_automatic_standard),
 );
 
-# Reduction primal/dual error vs. sparsity of solution
 
-dataSparsity = [trajectoryBPCG_standard, trajectoryBPCG_quadratic]
-labelSparsity = ["BPCG (Standard)", "BPCG (Direct)"]
+dataSparsity = [trajectoryBPCG_standard, trajectoryBPCG_quadratic_automatic, trajectoryBPCG_quadratic_automatic_standard]
+labelSparsity = ["BPCG (Standard)", "AS_Quad", "AS_Standard"]
 
-# Plot sparsity
-# plot_sparsity(dataSparsity, labelSparsity, legend_position=:topright)
 
 # Plot trajectories
 plot_trajectories(dataSparsity, labelSparsity,xscalelog=false)

examples/blended_pairwise_with_sparsify.jl

@@ -14,7 +14,6 @@ using FrankWolfe
 using LinearAlgebra
 using Random
 
-import Pkg
 import HiGHS
 
 # lp_solver = GLPK.Optimizer
@@ -110,6 +109,27 @@ active_set_sparse = FrankWolfe.ActiveSetSparsifier(FrankWolfe.ActiveSet([1.0], [
     callback=callback,
 );
 
+trajectoryBPCG_as_sparse = []
+callback = build_callback(trajectoryBPCG_as_sparse)
+
+active_set_sparse = FrankWolfe.ActiveSetSparsifier(FrankWolfe.ActiveSet([1.0], [x00], similar(x00)), HiGHS.Optimizer())
+
+@time x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient(
+    f,
+    grad!,
+    lmo,
+    copy(x00),
+    max_iteration=k,
+    line_search=FrankWolfe.Shortstep(2.0),
+    print_iter=k / 10,
+    memory_mode=FrankWolfe.InplaceEmphasis(),
+    verbose=true,
+    trajectory=true,
+    callback=callback,
+    sparsify=false,
+#    squadratic=true, # activate to see the effect of the numerical precision of the LP solver
+);
+
 # Reduction primal/dual error vs. sparsity of solution
 
 dataSparsity = [trajectoryBPCG_standard, trajectoryBPCG_quadratic]

examples/simple_stateless.jl

@@ -71,7 +71,7 @@ plot_trajectories([trajectory_simple[1:end], trajectory_restart[1:end]], ["simpl
 
 # simple step iterations about 33% faster
 
-@test line_search.factor <= 44
+@test line_search.factor <= 51
 
 x, v, primal, dual_gap, trajectory_restart_highpres = FrankWolfe.frank_wolfe(
     f,

src/FrankWolfe.jl

@@ -34,6 +34,7 @@ include("moi_oracle.jl")
 include("function_gradient.jl")
 include("active_set.jl")
 include("active_set_quadratic.jl")
+include("active_set_quadratic_direct_solve.jl")
 include("active_set_sparsifier.jl")
 
 include("blended_cg.jl")

src/active_set.jl

@@ -326,3 +326,7 @@ function compute_active_set_iterate!(active_set::AbstractActiveSet{<:ScaledHotVe
     end
     return active_set.x
 end
+
+function update_weights!(as::AbstractActiveSet, new_weights)
+    as.weights .= new_weights
+end