diff --git a/Project.toml b/Project.toml index f0bdd4fe8..ca298fa1a 100644 --- a/Project.toml +++ b/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"] diff --git a/examples/blended_pairwise_with_direct.jl b/examples/blended_pairwise_with_direct.jl index 576361db1..1304749d2 100644 --- a/examples/blended_pairwise_with_direct.jl +++ b/examples/blended_pairwise_with_direct.jl @@ -16,23 +16,14 @@ 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) @@ -40,9 +31,6 @@ 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) diff --git a/examples/blended_pairwise_with_direct_non-standard-quadratic.jl b/examples/blended_pairwise_with_direct_non-standard-quadratic.jl index ff7736be1..5710e92d9 100644 --- a/examples/blended_pairwise_with_direct_non-standard-quadratic.jl +++ b/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) diff --git a/examples/blended_pairwise_with_sparsify.jl b/examples/blended_pairwise_with_sparsify.jl index db2a98998..941df79cd 100644 --- a/examples/blended_pairwise_with_sparsify.jl +++ b/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] diff --git a/examples/simple_stateless.jl b/examples/simple_stateless.jl index 6a550a963..6bf230e26 100644 --- a/examples/simple_stateless.jl +++ b/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, diff --git a/src/FrankWolfe.jl b/src/FrankWolfe.jl index e17e40038..cc2fef145 100644 --- a/src/FrankWolfe.jl +++ b/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") diff --git a/src/active_set.jl b/src/active_set.jl index 62393b47f..c02145666 100644 --- a/src/active_set.jl +++ b/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 diff --git a/src/active_set_quadratic.jl b/src/active_set_quadratic.jl index d252d186d..5739198f2 100644 --- a/src/active_set_quadratic.jl +++ b/src/active_set_quadratic.jl @@ -49,7 +49,8 @@ function detect_quadratic_function(grad!, x0; test=true) end function ActiveSetQuadratic(tuple_values::AbstractVector{Tuple{R,AT}}, grad!::Function) where {AT,R} - return ActiveSetQuadratic(tuple_values, detect_quadratic_function(grad!, tuple_values[1][2])...) + A, b = detect_quadratic_function(grad!, tuple_values[1][2]) + return ActiveSetQuadratic(tuple_values, A, b) end function ActiveSetQuadratic(tuple_values::AbstractVector{Tuple{R,AT}}, A::H, b) where {AT,R,H} @@ -71,13 +72,15 @@ function ActiveSetQuadratic(tuple_values::AbstractVector{Tuple{R,AT}}, A::H, b) dots_b[idx] = fast_dot(b, atoms[idx]) end x = similar(b) - as = ActiveSetQuadratic{AT,R,typeof(x),H}(weights, atoms, x, A, b, dots_x, dots_A, dots_b, weights_prev, modified) + as = ActiveSetQuadratic(weights, atoms, x, A, b, dots_x, dots_A, dots_b, weights_prev, modified) + compute_active_set_iterate!(as) return as end -function ActiveSetQuadratic{AT,R}(tuple_values::AbstractVector{<:Tuple{<:Number,<:Any}}, grad!) where {AT,R} - return ActiveSetQuadratic{AT,R}(tuple_values, detect_quadratic_function(grad!, tuple_values[1][2])...) +function ActiveSetQuadratic{AT,R}(tuple_values::AbstractVector{<:Tuple{<:Number,<:Any}}, grad!::Function) where {AT,R} + A, b = detect_quadratic_function(grad!, tuple_values[1][2]) + return ActiveSetQuadratic{AT,R}(tuple_values, A, b) end function ActiveSetQuadratic{AT,R}(tuple_values::AbstractVector{<:Tuple{<:Number,<:Any}}, A::H, b) where {AT,R,H} @@ -116,6 +119,7 @@ function Base.:*(a::Identity, b) return a.λ * b end end + function ActiveSetQuadratic(tuple_values::AbstractVector{Tuple{R,AT}}, A::UniformScaling, b) where {AT,R} return ActiveSetQuadratic(tuple_values, Identity(A.λ), b) end @@ -147,6 +151,7 @@ function Base.push!(as::ActiveSetQuadratic{AT,R}, (λ, a)) where {AT,R} return as end +# TODO multi-indices version function Base.deleteat!(as::ActiveSetQuadratic, idx::Int) @inbounds for i in 1:idx-1 as.dots_x[i] -= as.weights_prev[idx] * as.dots_A[idx][i] @@ -310,3 +315,9 @@ function update_active_set_quadratic!(warm_as::ActiveSetQuadratic{AT,R,IT,H}, b) compute_active_set_iterate!(warm_as) return warm_as end + +function update_weights!(as::ActiveSetQuadratic, new_weights) + as.weights_prev .= as.weights + as.weights .= new_weights + as.modified .= true +end diff --git a/src/active_set_quadratic_direct_solve.jl b/src/active_set_quadratic_direct_solve.jl new file mode 100644 index 000000000..0d71db82f --- /dev/null +++ b/src/active_set_quadratic_direct_solve.jl @@ -0,0 +1,248 @@ + +""" + ActiveSetQuadraticLinearSolve{AT, R, IT} + +Represents an active set of extreme vertices collected in a FW algorithm, +along with their coefficients `(λ_i, a_i)`. +`R` is the type of the `λ_i`, `AT` is the type of the atoms `a_i`. +The iterate `x = ∑λ_i a_i` is stored in x with type `IT`. +The objective function is assumed to be of the form `f(x)=½⟨x,Ax⟩+⟨b,x⟩+c` +so that the gradient is `∇f(x)=Ax+b`. + +This active set stores an inner `active_set` that keeps track of the current set of vertices and convex decomposition. +It therefore delegates all update, deletion, and addition operations to this inner active set. +The `weight`, `atoms`, and `x` fields should only be accessed to read and are effectively the same objects as those in the inner active set. + +The structure also contains a scheduler struct which is called with the `should_solve_lp` function. +To define a new frequency at which the LP should be solved, one can define another scheduler struct and implement the corresponding method. +""" +struct ActiveSetQuadraticLinearSolve{AT, R <: Real, IT, H, BT, OT <: MOI.AbstractOptimizer, AS <: AbstractActiveSet, SF} <: AbstractActiveSet{AT,R,IT} + weights::Vector{R} + atoms::Vector{AT} + x::IT + A::H # Hessian matrix + b::BT # linear term + active_set::AS + lp_optimizer::OT + scheduler::SF + counter::Base.RefValue{Int} +end + +""" + ActiveSetQuadraticLinearSolve(tuple_values::Vector{Tuple{R,AT}}, grad!::Function, lp_optimizer) + +Creates an ActiveSetQuadraticLinearSolve by computing the Hessian and linear term from `grad!`. +""" +function ActiveSetQuadraticLinearSolve(tuple_values::Vector{Tuple{R,AT}}, grad!::Function, lp_optimizer; scheduler=LogScheduler()) where {AT,R} + A, b = detect_quadratic_function(grad!, tuple_values[1][2]) + return ActiveSetQuadraticLinearSolve(tuple_values, A, b, lp_optimizer, scheduler=scheduler) +end + +""" + ActiveSetQuadraticLinearSolve(tuple_values::Vector{Tuple{R,AT}}, A, b, lp_optimizer) + +Creates an `ActiveSetQuadraticLinearSolve` from the given Hessian `A`, linear term `b` and `lp_optimizer` by creating an inner `ActiveSetQuadratic` active set. +""" +function ActiveSetQuadraticLinearSolve(tuple_values::Vector{Tuple{R,AT}}, A::H, b, lp_optimizer; scheduler=LogScheduler()) where {AT,R,H} + inner_as = ActiveSetQuadratic(tuple_values, A, b) + return ActiveSetQuadraticLinearSolve(inner_as.weights, inner_as.atoms, inner_as.x, inner_as.A, inner_as.b, inner_as, lp_optimizer, scheduler, Ref(0)) +end + +function ActiveSetQuadraticLinearSolve(inner_as::AbstractActiveSet, A, b, lp_optimizer; scheduler=LogScheduler()) + as = ActiveSetQuadraticLinearSolve(inner_as.weights, inner_as.atoms, inner_as.x, A, b, inner_as, lp_optimizer, scheduler, Ref(0)) + compute_active_set_iterate!(as) + return as +end + +function ActiveSetQuadraticLinearSolve(inner_as::AbstractActiveSet, A::LinearAlgebra.UniformScaling, b, lp_optimizer; scheduler=LogScheduler()) + as = ActiveSetQuadraticLinearSolve(inner_as.weights, inner_as.atoms, inner_as.x, A, b, inner_as, lp_optimizer, scheduler, Ref(0)) + compute_active_set_iterate!(as) + return as +end + +function ActiveSetQuadraticLinearSolve(inner_as::AbstractActiveSet, grad!::Function, lp_optimizer; scheduler=LogScheduler()) + A, b = detect_quadratic_function(grad!, inner_as.atoms[1]) + return ActiveSetQuadraticLinearSolve(inner_as, A, b, lp_optimizer; scheduler=scheduler) +end + +function ActiveSetQuadraticLinearSolve{AT,R}(tuple_values::Vector{<:Tuple{<:Number,<:Any}}, grad!::Function, lp_optimizer; scheduler=LogScheduler()) where {AT,R} + A, b = detect_quadratic_function(grad!, tuple_values[1][2]) + return ActiveSetQuadraticLinearSolve{AT,R}(tuple_values, A, b, lp_optimizer; scheduler=scheduler) +end + +function ActiveSetQuadraticLinearSolve{AT,R}(tuple_values::Vector{<:Tuple{<:Number,<:Any}}, A::H, b, lp_optimizer; scheduler=LogScheduler()) where {AT,R,H} + inner_as = ActiveSetQuadratic{AT,R}(tuple_values, A, b) + as = ActiveSetQuadraticLinearSolve{AT,R,typeof(x),H}(inner_as.weights, inner_as.atoms, inner_as.x, A, b, inner_as, lp_optimizer, scheduler, Ref(0)) + compute_active_set_iterate!(as) + return as +end + +function ActiveSetQuadraticLinearSolve(tuple_values::Vector{Tuple{R,AT}}, A::UniformScaling, b, lp_optimizer; scheduler=LogScheduler()) where {AT,R} + return ActiveSetQuadraticLinearSolve(tuple_values, Identity(A.λ), b, lp_optimizer; scheduler=scheduler) +end +function ActiveSetQuadraticLinearSolve{AT,R}(tuple_values::Vector{<:Tuple{<:Number,<:Any}}, A::UniformScaling, b, lp_optimizer; scheduler=LogScheduler()) where {AT,R} + return ActiveSetQuadraticLinearSolve{AT,R}(tuple_values, Identity(A.λ), b, lp_optimizer; scheduler=scheduler) +end + +# all mutating functions are delegated to the inner active set + +Base.push!(as::ActiveSetQuadraticLinearSolve, tuple) = push!(as.active_set, tuple) + +Base.deleteat!(as::ActiveSetQuadraticLinearSolve, idx) = deleteat!(as.active_set, idx) + +Base.empty!(as::ActiveSetQuadraticLinearSolve) = empty!(as.active_set) + +function active_set_update!( + as::ActiveSetQuadraticLinearSolve{AT,R}, + lambda, atom, renorm=true, idx=nothing; + weight_purge_threshold=weight_purge_threshold_default(R), + add_dropped_vertices=false, + vertex_storage=nothing, +) where {AT,R} + if idx === nothing + idx = find_atom(as, atom) + end + active_set_update!(as.active_set, lambda, atom, renorm, idx; weight_purge_threshold=weight_purge_threshold, add_dropped_vertices=add_dropped_vertices, vertex_storage=vertex_storage) + # new atom introduced, we can solve the auxiliary LP + if idx < 0 + as.counter[] += 1 + if should_solve_lp(as, as.scheduler) + solve_quadratic_activeset_lp!(as) + end + end + return as +end + +active_set_renormalize!(as::ActiveSetQuadraticLinearSolve) = active_set_renormalize!(as.active_set) + +function active_set_argmin(as::ActiveSetQuadraticLinearSolve, direction) + return active_set_argmin(as.active_set, direction) +end + +function active_set_argminmax(as::ActiveSetQuadraticLinearSolve, direction; Φ=0.5) + return active_set_argminmax(as.active_set, direction; Φ=Φ) +end + +# generic quadratic with quadratic information provided + +""" + solve_quadratic_activeset_lp!(as::ActiveSetQuadraticLinearSolve{AT, R, IT, H})) + +Solves the auxiliary LP over the current active set. +The method is specialized by type `H` of the Hessian matrix `A`. +""" +function solve_quadratic_activeset_lp!(as::ActiveSetQuadraticLinearSolve{AT, R, IT, <: AbstractMatrix}) where {AT, R, IT} + nv = length(as) + o = as.lp_optimizer + MOI.empty!(o) + λ = MOI.add_variables(o, nv) + # λ ≥ 0, ∑ λ == 1 + MOI.add_constraint.(o, λ, MOI.GreaterThan(0.0)) + sum_of_variables = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, λ), 0.0) + MOI.add_constraint(o, sum_of_variables, MOI.EqualTo(1.0)) + # Vᵗ A V λ == -Vᵗ b + for atom in as.atoms + lhs = MOI.ScalarAffineFunction{Float64}([], 0.0) + Base.sizehint!(lhs.terms, nv) + # replaces direct sum because of MOI and MutableArithmetic slow sums + for j in 1:nv + push!(lhs.terms, MOI.ScalarAffineTerm(fast_dot(atom, as.A, as.atoms[j]), λ[j])) + end + rhs = -dot(atom, as.b) + MOI.add_constraint(o, lhs, MOI.EqualTo(rhs)) + end + MOI.set(o, MOI.ObjectiveFunction{typeof(sum_of_variables)}(), sum_of_variables) + MOI.set(o, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.optimize!(o) + if MOI.get(o, MOI.TerminationStatus()) ∉ (MOI.OPTIMAL, MOI.FEASIBLE_POINT, MOI.ALMOST_OPTIMAL) + return as + end + indices_to_remove = Int[] + new_weights = R[] + for idx in eachindex(λ) + weight_value = MOI.get(o, MOI.VariablePrimal(), λ[idx]) + if weight_value <= 1e-5 + push!(indices_to_remove, idx) + else + push!(new_weights, weight_value) + end + end + deleteat!(as.active_set, indices_to_remove) + @assert length(as) == length(new_weights) + update_weights!(as.active_set, new_weights) + active_set_cleanup!(as) + active_set_renormalize!(as) + compute_active_set_iterate!(as) + return as +end + +# special case of scaled identity Hessian + +function solve_quadratic_activeset_lp!(as::ActiveSetQuadraticLinearSolve{AT, R, IT, <: Union{Identity, LinearAlgebra.UniformScaling}}) where {AT, R, IT} + hessian_scaling = as.A.λ + nv = length(as) + o = as.lp_optimizer + MOI.empty!(o) + λ = MOI.add_variables(o, nv) + # λ ≥ 0, ∑ λ == 1 + MOI.add_constraint.(o, λ, MOI.GreaterThan(0.0)) + sum_of_variables = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, λ), 0.0) + MOI.add_constraint(o, sum_of_variables, MOI.EqualTo(1.0)) + # a Aᵗ A λ == -Aᵗ b + for atom in as.atoms + lhs = MOI.ScalarAffineFunction{Float64}([], 0.0) + Base.sizehint!(lhs.terms, nv) + # replaces direct sum because of MOI and MutableArithmetic slow sums + for j in 1:nv + push!(lhs.terms, MOI.ScalarAffineTerm(hessian_scaling * dot(atom, as.atoms[j]), λ[j])) + end + rhs = -dot(atom, as.b) + MOI.add_constraint(o, lhs, MOI.EqualTo(rhs)) + end + MOI.set(o, MOI.ObjectiveFunction{typeof(sum_of_variables)}(), sum_of_variables) + MOI.set(o, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.optimize!(o) + if MOI.get(o, MOI.TerminationStatus()) ∉ (MOI.OPTIMAL, MOI.FEASIBLE_POINT, MOI.ALMOST_OPTIMAL) + return as + end + indices_to_remove = Int[] + new_weights = R[] + for idx in eachindex(λ) + weight_value = MOI.get(o, MOI.VariablePrimal(), λ[idx]) + if weight_value <= 2 * weight_purge_threshold_default(typeof(weight_value)) + push!(indices_to_remove, idx) + else + push!(new_weights, weight_value) + end + end + deleteat!(as.active_set, indices_to_remove) + @assert length(as) == length(new_weights) + update_weights!(as.active_set, new_weights) + active_set_cleanup!(as) + @assert all(>=(0), new_weights) + active_set_renormalize!(as) + compute_active_set_iterate!(as) + return as +end + +struct LogScheduler{T} + start_time::Int + scaling_factor::T + max_interval::Int + current_interval::Base.RefValue{Int} + last_solve_counter::Base.RefValue{Int} +end + +LogScheduler(; start_time=20, scaling_factor=1.5, max_interval=1000) = LogScheduler(start_time, scaling_factor, max_interval, Ref(start_time), Ref(0)) + +function should_solve_lp(as::ActiveSetQuadraticLinearSolve, scheduler::LogScheduler) + if as.counter[] - scheduler.last_solve_counter[] >= scheduler.current_interval[] + scheduler.last_solve_counter[] = as.counter[] + scheduler.current_interval[] = min( + round(Int, scheduler.scaling_factor * scheduler.current_interval[]), + scheduler.max_interval, + ) + return true + end + return false +end diff --git a/src/active_set_sparsifier.jl b/src/active_set_sparsifier.jl index 6bb523217..3759adf88 100644 --- a/src/active_set_sparsifier.jl +++ b/src/active_set_sparsifier.jl @@ -1,35 +1,35 @@ struct ActiveSetSparsifier{AT,R,IT,AS<:AbstractActiveSet{AT,R,IT},OT<:MOI.AbstractOptimizer} <: - AbstractActiveSet{AT,R,IT} - active_set::AS - weights::Vector{R} - atoms::Vector{AT} - x::IT - optimizer::OT - minimum_vertices::Int - solve_frequency::Int - counter::Base.RefValue{Int} + AbstractActiveSet{AT,R,IT} + active_set::AS + weights::Vector{R} + atoms::Vector{AT} + x::IT + optimizer::OT + minimum_vertices::Int + solve_frequency::Int + counter::Base.RefValue{Int} end function ActiveSetSparsifier( - active_set::AbstractActiveSet, - optimizer::MOI.AbstractOptimizer; - minimum_vertices=50, - solve_frequency=50, + active_set::AbstractActiveSet, + optimizer::MOI.AbstractOptimizer; + minimum_vertices=50, + solve_frequency=50, ) - return ActiveSetSparsifier( - active_set, - active_set.weights, - active_set.atoms, - active_set.x, - optimizer, - minimum_vertices, - solve_frequency, - Ref(0), - ) + return ActiveSetSparsifier( + active_set, + active_set.weights, + active_set.atoms, + active_set.x, + optimizer, + minimum_vertices, + solve_frequency, + Ref(0), + ) end function Base.push!(as::ActiveSetSparsifier, (λ, a)) - return push!(as.active_set, (λ, a)) + return push!(as.active_set, (λ, a)) end Base.deleteat!(as::ActiveSetSparsifier, idx::Int) = deleteat!(as.active_set, idx) @@ -37,60 +37,60 @@ Base.deleteat!(as::ActiveSetSparsifier, idx::Int) = deleteat!(as.active_set, idx Base.empty!(as::ActiveSetSparsifier) = empty!(as.active_set) function active_set_update!( - as::ActiveSetSparsifier{AS,OT,AT,R,IT}, - lambda, - atom, - renorm=true, - idx=nothing; - kwargs..., + as::ActiveSetSparsifier{AS,OT,AT,R,IT}, + lambda, + atom, + renorm=true, + idx=nothing; + kwargs..., ) where {AS,OT,AT,R,IT} - active_set_update!(as.active_set, lambda, atom, renorm, idx; kwargs...) - n = length(as) - as.counter[] += 1 - if n > as.minimum_vertices && mod(as.counter[], as.solve_frequency) == 0 - # sparsifying active set - MOI.empty!(as.optimizer) - x0 = as.active_set.x - λ = MOI.add_variables(as.optimizer, n) - # λ ∈ Δ_n ⇔ λ ≥ 0, ∑ λ == 1 - MOI.add_constraint.(as.optimizer, λ, MOI.GreaterThan(0.0)) - MOI.add_constraint(as.optimizer, sum(λ; init=0.0), MOI.EqualTo(1.0)) - x_sum = 0 * as.active_set.atoms[1] - for (idx, atom) in enumerate(as.active_set.atoms) - x_sum += λ[idx] * atom - end - for idx in eachindex(x_sum) - MOI.add_constraint(as.optimizer, x_sum[idx], MOI.EqualTo(x0[idx])) - end - # Set a dummy objective (minimize ∑λ) - dummy_objective = sum(λ; init=0.0) - MOI.set(as.optimizer, MOI.ObjectiveFunction{typeof(dummy_objective)}(), dummy_objective) - MOI.set(as.optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE) - MOI.optimize!(as.optimizer) - if MOI.get(as.optimizer, MOI.TerminationStatus()) in - (MOI.OPTIMAL, MOI.FEASIBLE_POINT, MOI.ALMOST_OPTIMAL) - indices_to_remove = Int[] - new_weights = R[] - for idx in eachindex(λ) - weight_value = MOI.get(as.optimizer, MOI.VariablePrimal(), λ[idx]) - if weight_value <= sqrt(weight_purge_threshold_default(typeof(weight_value))) - push!(indices_to_remove, idx) - else - push!(new_weights, weight_value) - end - end - deleteat!(as.active_set.atoms, indices_to_remove) - deleteat!(as.active_set.weights, indices_to_remove) - @assert length(as) == length(new_weights) - as.active_set.weights .= new_weights - active_set_renormalize!(as) - end - end - return as + active_set_update!(as.active_set, lambda, atom, renorm, idx; kwargs...) + n = length(as) + as.counter[] += 1 + if n > as.minimum_vertices && mod(as.counter[], as.solve_frequency) == 0 + # sparsifying active set + MOI.empty!(as.optimizer) + x0 = as.active_set.x + λ = MOI.add_variables(as.optimizer, n) + # λ ∈ Δ_n ⇔ λ ≥ 0, ∑ λ == 1 + MOI.add_constraint.(as.optimizer, λ, MOI.GreaterThan(0.0)) + MOI.add_constraint(as.optimizer, sum(λ; init=0.0), MOI.EqualTo(1.0)) + x_sum = 0 * as.active_set.atoms[1] + for (idx, atom) in enumerate(as.active_set.atoms) + x_sum += λ[idx] * atom + end + for idx in eachindex(x_sum) + MOI.add_constraint(as.optimizer, x_sum[idx], MOI.EqualTo(x0[idx])) + end + # Set a dummy objective (minimize ∑λ) + dummy_objective = sum(λ; init=0.0) + MOI.set(as.optimizer, MOI.ObjectiveFunction{typeof(dummy_objective)}(), dummy_objective) + MOI.set(as.optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE) + MOI.optimize!(as.optimizer) + if MOI.get(as.optimizer, MOI.TerminationStatus()) in + (MOI.OPTIMAL, MOI.FEASIBLE_POINT, MOI.ALMOST_OPTIMAL) + indices_to_remove = Int[] + new_weights = R[] + for idx in eachindex(λ) + weight_value = MOI.get(as.optimizer, MOI.VariablePrimal(), λ[idx]) + if weight_value <= sqrt(weight_purge_threshold_default(typeof(weight_value))) + push!(indices_to_remove, idx) + else + push!(new_weights, weight_value) + end + end + deleteat!(as.active_set.atoms, indices_to_remove) + deleteat!(as.active_set.weights, indices_to_remove) + @assert length(as) == length(new_weights) + as.active_set.weights .= new_weights + active_set_renormalize!(as) + end + end + return as end active_set_renormalize!(as::ActiveSetSparsifier) = active_set_renormalize!(as.active_set) active_set_argmin(as::ActiveSetSparsifier, direction) = active_set_argmin(as.active_set, direction) active_set_argminmax(as::ActiveSetSparsifier, direction; Φ=0.5) = - active_set_argminmax(as.active_set, direction; Φ=Φ) + active_set_argminmax(as.active_set, direction; Φ=Φ) \ No newline at end of file diff --git a/src/blended_pairwise.jl b/src/blended_pairwise.jl index a1341135a..40435d946 100644 --- a/src/blended_pairwise.jl +++ b/src/blended_pairwise.jl @@ -30,10 +30,6 @@ function blended_pairwise_conditional_gradient( add_dropped_vertices=false, use_extra_vertex_storage=false, recompute_last_vertex=true, - squadratic=false, # standard quadratic objective - quadratic=false, # quadratic objective but not standard quadratic - lp_solver=nothing, # lp_solver used for direct solve or sparsification - ds_scale_up=2, # direct solve scale up parameter ) # add the first vertex to active set from initialization active_set = ActiveSet([(1.0, x0)]) @@ -63,10 +59,6 @@ function blended_pairwise_conditional_gradient( add_dropped_vertices=add_dropped_vertices, use_extra_vertex_storage=use_extra_vertex_storage, recompute_last_vertex=recompute_last_vertex, - squadratic=squadratic, - quadratic=quadratic, # Pass the new parameter - lp_solver=lp_solver, - ds_scale_up=ds_scale_up, # Pass the new parameter ) end @@ -100,10 +92,6 @@ function blended_pairwise_conditional_gradient( add_dropped_vertices=false, use_extra_vertex_storage=false, recompute_last_vertex=true, - squadratic=false, # non-standard quadratic objective - quadratic=false, # quadratic objective but not standard quadratic - lp_solver=nothing, # lp_solver used for direct solve or sparsification - ds_scale_up=2, # direct solve scale up parameter ) where {AT,R} # format string for output of the algorithm @@ -132,8 +120,6 @@ function blended_pairwise_conditional_gradient( end t = 0 - next_direct_solve_interval = 10 - last_direct_solve = 0 compute_active_set_iterate!(active_set) x = get_active_set_iterate(active_set) primal = convert(eltype(x), Inf) @@ -146,10 +132,6 @@ function blended_pairwise_conditional_gradient( gradient = collect(x) end - if squadratic || quadratic # when using LPs then lp_solver must be provided - @assert lp_solver !== nothing "When squadratic, sparsify, or quadratic is true, lp_solver must be provided" - end - if verbose println("\nBlended Pairwise Conditional Gradient Algorithm.") NumType = eltype(x) @@ -159,9 +141,6 @@ function blended_pairwise_conditional_gradient( grad_type = typeof(gradient) println("GRADIENTTYPE: $grad_type LAZY: $lazy lazy_tolerance: $lazy_tolerance") println("LMO: $(typeof(lmo))") - if lp_solver !== nothing - println("LP SOLVER: $(lp_solver)") - end if use_extra_vertex_storage && !lazy @info("vertex storage only used in lazy mode") end @@ -312,20 +291,6 @@ function blended_pairwise_conditional_gradient( end vertex_taken = v # short circuit for quadratic case - if (squadratic || quadratic) && ((t > last_direct_solve && t - last_direct_solve >= next_direct_solve_interval)) # Updated condition - next_direct_solve_interval *= ds_scale_up - last_direct_solve = t - if quadratic - active_set = direct_solve_gq(active_set, grad!, lp_solver) # direct solve for non-standard quadratic problem - else - active_set = direct_solve(active_set, grad!, lp_solver) # direct solve for standard quadratic problem - end - active_set_cleanup!(active_set; weight_purge_threshold=weight_purge_threshold) - compute_active_set_iterate!(active_set) - x = get_active_set_iterate(active_set) - grad!(gradient, x) - step_type = ST_DIRECT - end dual_gap = fast_dot(gradient, x) - fast_dot(gradient, v) # if we are about to exit, compute dual_gap with the cleaned-up x if dual_gap ≤ epsilon @@ -519,4 +484,4 @@ function blended_pairwise_conditional_gradient( end return (x=x, v=v, primal=primal, dual_gap=dual_gap, traj_data=traj_data, active_set=active_set) -end \ No newline at end of file +end diff --git a/src/utils.jl b/src/utils.jl index bf707b5ec..a4b18c4a6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -163,6 +163,77 @@ function fast_dot(A::Matrix{T1}, B::SparseArrays.SparseMatrixCSC{T2}) where {T1, return s end +fast_dot(a, Q, b) = dot(a, Q, b) + +function fast_dot(a::SparseArrays.AbstractSparseVector{<:Real}, Q::Diagonal{<:Real}, b::AbstractVector{<:Real}) + if a === b + return _fast_quadratic_form_symmetric(a, Q) + end + d = Q.diag + nzvals = SparseArrays.nonzeros(a) + nzinds = SparseArrays.nonzeroinds(a) + return sum(eachindex(nzvals); init=zero(eltype(a))) do nzidx + nzvals[nzidx] * d[nzinds[nzidx]] * b[nzinds[nzidx]] + end +end + +function fast_dot(a::SparseArrays.AbstractSparseVector{<:Real}, Q::Diagonal{<:Real}, b::SparseArrays.AbstractSparseVector{<:Real}) + if a === b + return _fast_quadratic_form_symmetric(a, Q) + end + n = length(a) + if length(b) != n + throw( + DimensionMismatch("Vector a has a length $n but b has a length $(length(b))") + ) + end + anzind = SparseArrays.nonzeroinds(a) + bnzind = SparseArrays.nonzeroinds(b) + anzval = SparseArrays.nonzeros(a) + bnzval = SparseArrays.nonzeros(b) + s = zero(Base.promote_eltype(a, Q, b)) + + if isempty(anzind) || isempty(bnzind) + return s + end + + a_idx = 1 + b_idx = 1 + a_idx_last = length(anzind) + b_idx_last = length(bnzind) + + # go through the nonzero indices of a and b simultaneously + @inbounds while a_idx <= a_idx_last && b_idx <= b_idx_last + ia = anzind[a_idx] + ib = bnzind[b_idx] + if ia == ib + s += dot(anzval[a_idx], Q.diag[ia], bnzval[b_idx]) + a_idx += 1 + b_idx += 1 + elseif ia < ib + a_idx += 1 + else + b_idx += 1 + end + end + return s +end + + +function _fast_quadratic_form_symmetric(a, Q) + d = Q.diag + if length(d) != length(a) + throw(DimensionMismatch()) + end + nzvals = SparseArrays.nonzeros(a) + nzinds = SparseArrays.nonzeroinds(a) + s = zero(Base.promote_eltype(a, Q)) + @inbounds for nzidx in eachindex(nzvals) + s += nzvals[nzidx]^2 * d[nzinds[nzidx]] + end + return s +end + """ trajectory_callback(storage) diff --git a/test/quadratic_lp_active_set.jl b/test/quadratic_lp_active_set.jl new file mode 100644 index 000000000..a86c0b713 --- /dev/null +++ b/test/quadratic_lp_active_set.jl @@ -0,0 +1,110 @@ +using FrankWolfe +using LinearAlgebra +using Random +using Test +import HiGHS +import MathOptInterface as MOI + +n = Int(1e4) +k = 5000 + +s = 10 +@info "Seed $s" +Random.seed!(s) + +xpi = rand(n); +total = sum(xpi); + +const xp = xpi ./ total; + +f(x) = norm(x - xp)^2 +function grad!(storage, x) + @. storage = 2 * (x - xp) +end + +lmo = FrankWolfe.KSparseLMO(5, 1.0) + +const x00 = FrankWolfe.compute_extreme_point(lmo, rand(n)) + +function build_callback(trajectory_arr) + return function callback(state, active_set, args...) + return push!(trajectory_arr, (FrankWolfe.callback_state(state)..., length(active_set))) + end +end + +trajectoryBPCG_standard = [] +x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient( + f, + grad!, + lmo, + copy(x00), + max_iteration=k, + line_search=FrankWolfe.Shortstep(2.0), + verbose=false, + callback=build_callback(trajectoryBPCG_standard), +); + +as_quad_direct = FrankWolfe.ActiveSetQuadraticLinearSolve( + [(1.0, copy(x00))], + 2 * LinearAlgebra.I, -2xp, + MOI.instantiate(MOI.OptimizerWithAttributes(HiGHS.Optimizer, MOI.Silent() => true)), +) + +trajectoryBPCG_quadratic_direct_specialized = [] +x, v, primal, dual_gap, _ = FrankWolfe.blended_pairwise_conditional_gradient( + f, + grad!, + lmo, + as_quad_direct, + max_iteration=k, + line_search=FrankWolfe.Shortstep(2.0), + verbose=false, + callback=build_callback(trajectoryBPCG_quadratic_direct_specialized), +); + +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)), +) + +trajectoryBPCG_quadratic_direct_generic = [] +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)), +) + +trajectoryBPCG_quadratic_noqas = [] +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), +); + + +dual_gaps_quadratic_specialized = getindex.(trajectoryBPCG_quadratic_direct_specialized, 4) +dual_gaps_quadratic_generic = getindex.(trajectoryBPCG_quadratic_direct_generic, 4) +dual_gaps_quadratic_noqas = getindex.(trajectoryBPCG_quadratic_noqas, 4) +dual_gaps_bpcg = getindex.(trajectoryBPCG_standard, 4) + + +@test dual_gaps_quadratic_specialized[end] < dual_gaps_bpcg[end] +@test dual_gaps_quadratic_noqas[end] < dual_gaps_bpcg[end] +@test norm(dual_gaps_quadratic_noqas - dual_gaps_quadratic_noqas) ≤ k * 1e-5 diff --git a/test/runtests.jl b/test/runtests.jl index 360b7cd19..40f5ac366 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -798,6 +798,13 @@ using Test end end +module LpDirectSolveTest +using Test +@testset "LP solving for quadratic functions and active set" begin + include("quadratic_lp_active_set.jl") +end +end + module SparsifyingActiveSetTest using Test @testset "Sparsifying active set" begin