plot_sparsity (generic function with 1 method)Accelerations for quadratic functions and symmetric problems
This example illustrates how to exploit symmetry to reduce the dimension of the problem via SymmetricLMO. Moreover, active set based algorithms can be accelerated by using the specialized structure ActiveSetQuadratic.
The specific problem we consider here comes from quantum information and some context can be found here. Formally, we want to find the distance between a tensor of size m^N and the N-partite local polytope which is defined by its vertices
\[d^{\vec{a}^{(1)}\ldots \vec{a}^{(N)}}_{x_1\ldots x_N}\coloneqq\prod_{n=1}^Na^{(n)}_{x_n}\]
labeled by $\vec{a}^{(n)}=a^{(n)}_1\ldots a^{(n)}_m$ for $n\in[1,N]$, where $a^{(n)}_x=\pm1$. In the bipartite case (N=2), this polytope is affinely equivalent to the cut polytope.
Import and setup
We first import the necessary packages.
import Combinatorics
import FrankWolfe
import LinearAlgebra
import TullioThen we can define our custom LMO, together with the method compute_extreme_point, which simply enumerates the vertices $d^{\vec{a}^{(1)}}$ defined above. This structure is specialized for the case N=5 and contains pre-allocated fields used to accelerate the enumeration. Note that the output type (full tensor) is quite naive, but this is enough to illustrate the syntax in this toy example.
struct BellCorrelationsLMO{T} <: FrankWolfe.LinearMinimizationOracle
m::Int # size of the tensor
tmp1::Array{T, 1}
tmp2::Array{T, 2}
tmp3::Array{T, 3}
tmp4::Array{T, 4}
end
function FrankWolfe.compute_extreme_point(lmo::BellCorrelationsLMO{T}, A::Array{T, 5}; kwargs...) where {T <: Number}
ax = [ones(T, lmo.m) for n in 1:5]
sc1 = zero(T)
sc2 = one(T)
axm = [zeros(T, lmo.m) for n in 1:5]
scm = typemax(T)
L = 2^lmo.m
aux = zeros(Int, lmo.m)
for λa5 in 0:(L÷2)-1
digits!(aux, λa5, base=2)
ax[5] .= 2aux .- 1
Tullio.@tullio lmo.tmp4[x1, x2, x3, x4] = A[x1, x2, x3, x4, x5] * ax[5][x5]
for λa4 in 0:L-1
digits!(aux, λa4, base=2)
ax[4] .= 2aux .- 1
Tullio.@tullio lmo.tmp3[x1, x2, x3] = lmo.tmp4[x1, x2, x3, x4] * ax[4][x4]
for λa3 in 0:L-1
digits!(aux, λa3, base=2)
ax[3] .= 2aux .- 1
Tullio.@tullio lmo.tmp2[x1, x2] = lmo.tmp3[x1, x2, x3] * ax[3][x3]
for λa2 in 0:L-1
digits!(aux, λa2, base=2)
ax[2] .= 2aux .- 1
LinearAlgebra.mul!(lmo.tmp1, lmo.tmp2, ax[2])
for x1 in 1:lmo.m
ax[1][x1] = lmo.tmp1[x1] > zero(T) ? -one(T) : one(T)
end
sc = LinearAlgebra.dot(ax[1], lmo.tmp1)
if sc < scm
scm = sc
for n in 1:5
axm[n] .= ax[n]
end
end
end
end
end
end
return [axm[1][x1]*axm[2][x2]*axm[3][x3]*axm[4][x4]*axm[5][x5] for x1 in 1:lmo.m, x2 in 1:lmo.m, x3 in 1:lmo.m, x4 in 1:lmo.m, x5 in 1:lmo.m]
endThen we define our specific instance, coming from a GHZ state measured with measurements forming a regular polygon on the equator of the Bloch sphere. See this article for definitions and references.
function correlation_tensor_GHZ_polygon(::Type{T}, N::Int, m::Int) where {T <: Number}
res = zeros(T, m*ones(Int, N)...)
tab_cos = [cos(x*T(pi)/m) for x in 0:N*m]
tab_cos[abs.(tab_cos) .< Base.rtoldefault(T)] .= zero(T)
for ci in CartesianIndices(res)
res[ci] = tab_cos[sum(ci.I)-N+1]
end
return res
end
T = Float64
verbose = true
max_iteration = 10^4
m = 5
p = 0.23correlation_tensor_GHZ_polygon(T, 5, m)
x0 = zeros(T, size(p))The objective function is simply $\frac12\|x-p\|_2^2$, which we decompose in different terms for speed.
normp2 = LinearAlgebra.dot(p, p) / 2
f = let p = p, normp2 = normp2
x -> LinearAlgebra.dot(x, x) / 2 - LinearAlgebra.dot(p, x) + normp2
end
grad! = let p = p
(storage, x) -> begin
@inbounds for i in eachindex(x)
storage[i] = x[i] - p[i]
end
end
endNaive run
If we run the blended pairwise conditional gradient algorithm without modifications, convergence is not reached in 10000 iterations.
lmo_naive = BellCorrelationsLMO{T}(m, zeros(T, m), zeros(T, m, m), zeros(T, m, m, m), zeros(T, m, m, m, m))
as_naive = FrankWolfe.ActiveSet([(one(T), x0)])
@time FrankWolfe.blended_pairwise_conditional_gradient(f, grad!, lmo_naive, as_naive; verbose, lazy=true, line_search=FrankWolfe.Shortstep(one(T)), max_iteration)
Blended Pairwise Conditional Gradient Algorithm.
MEMORY_MODE: FrankWolfe.InplaceEmphasis() STEPSIZE: Shortstep EPSILON: 1.0e-7 MAXITERATION: 10000 TYPE: Float64
GRADIENTTYPE: Array{Float64, 5} LAZY: true lazy_tolerance: 2.0
LMO: Main.BellCorrelationsLMO{Float64}
----------------------------------------------------------------------------------------------------------------
Type Iteration Primal Dual Dual Gap Time It/sec #ActiveSet
----------------------------------------------------------------------------------------------------------------
I 1 4.132813e+01 -4.029553e+01 8.162365e+01 0.000000e+00 Inf 1
LD 46 1.102218e+01 -2.888659e+01 3.990877e+01 1.078171e+01 4.266485e+00 46
LD 112 4.160807e+00 -1.560197e+01 1.976278e+01 2.124813e+01 5.271053e+00 88
LD 215 1.198957e+00 -4.345338e+00 5.544295e+00 3.092935e+01 6.951326e+00 126
LD 369 7.209992e-01 -2.020307e+00 2.741306e+00 3.967060e+01 9.301599e+00 160
P 1000 4.811978e-01 -2.260109e+00 2.741306e+00 7.233756e+01 1.382408e+01 289
P 2000 2.153509e-01 -2.525956e+00 2.741306e+00 1.164784e+02 1.717056e+01 451
LD 2126 1.869199e-01 -1.162610e+00 1.349530e+00 1.240793e+02 1.713420e+01 479
P 3000 1.233569e-01 -1.226173e+00 1.349530e+00 1.310373e+02 2.289425e+01 498
P 4000 5.988011e-02 -1.289650e+00 1.349530e+00 1.591784e+02 2.512903e+01 607
LD 4561 2.768136e-02 -3.409421e-01 3.686235e-01 1.728020e+02 2.639437e+01 653
P 5000 1.825337e-02 -3.503701e-01 3.686235e-01 1.733525e+02 2.884296e+01 625
LD 5916 1.173454e-02 -8.843410e-02 1.001686e-01 1.740029e+02 3.399943e+01 625
P 6000 1.147967e-02 -8.868896e-02 1.001686e-01 1.742948e+02 3.442443e+01 625
P 7000 9.885996e-03 -9.028264e-02 1.001686e-01 1.750002e+02 3.999995e+01 625
LD 7903 9.475200e-03 -1.638953e-02 2.586473e-02 1.756650e+02 4.498904e+01 625
P 8000 9.453302e-03 -1.641143e-02 2.586473e-02 1.759657e+02 4.546340e+01 625
P 9000 9.327751e-03 -1.653698e-02 2.586473e-02 1.766691e+02 5.094269e+01 625
P 10000 9.290467e-03 -1.657427e-02 2.586473e-02 1.773698e+02 5.637938e+01 625
Last 10001 9.290427e-03 1.176967e-03 8.113461e-03 1.776883e+02 5.628395e+01 625
----------------------------------------------------------------------------------------------------------------
PP 10001 9.290427e-03 1.176967e-03 8.113461e-03 1.779855e+02 5.618998e+01 625
----------------------------------------------------------------------------------------------------------------
178.220923 seconds (3.29 G allocations: 76.985 GiB, 4.89% gc time, 0.03% compilation time)Faster active set for quadratic functions
A first acceleration can be obtained by using the active set specialized for the quadratic objective function, whose gradient is here $x-p$, explaining the hessian and linear part provided as arguments. The speedup is obtained by pre-computing some scalar products to quickly obtained, in each iteration, the best and worst atoms currently in the active set.
asq_naive = FrankWolfe.ActiveSetQuadratic([(one(T), x0)], LinearAlgebra.I, -p)
@time FrankWolfe.blended_pairwise_conditional_gradient(f, grad!, lmo_naive, asq_naive; verbose, lazy=true, line_search=FrankWolfe.Shortstep(one(T)), max_iteration)
Blended Pairwise Conditional Gradient Algorithm.
MEMORY_MODE: FrankWolfe.InplaceEmphasis() STEPSIZE: Shortstep EPSILON: 1.0e-7 MAXITERATION: 10000 TYPE: Float64
GRADIENTTYPE: Array{Float64, 5} LAZY: true lazy_tolerance: 2.0
LMO: Main.BellCorrelationsLMO{Float64}
----------------------------------------------------------------------------------------------------------------
Type Iteration Primal Dual Dual Gap Time It/sec #ActiveSet
----------------------------------------------------------------------------------------------------------------
I 1 4.132813e+01 -4.029553e+01 8.162365e+01 0.000000e+00 Inf 1
LD 46 1.102218e+01 -2.888659e+01 3.990877e+01 1.082342e+01 4.250043e+00 46
LD 106 4.630016e+00 -1.445646e+01 1.908648e+01 1.899251e+01 5.581146e+00 79
LD 225 1.132318e+00 -4.781394e+00 5.913712e+00 3.044171e+01 7.391174e+00 126
LD 363 7.274406e-01 -2.207406e+00 2.934847e+00 3.870982e+01 9.377465e+00 157
P 1000 4.624255e-01 -2.472421e+00 2.934847e+00 7.539692e+01 1.326314e+01 306
LD 1836 2.172535e-01 -1.224204e+00 1.441458e+00 1.205135e+02 1.523481e+01 486
P 2000 1.982804e-01 -1.243177e+00 1.441458e+00 1.207776e+02 1.655936e+01 484
P 3000 1.217963e-01 -1.319661e+00 1.441458e+00 1.333047e+02 2.250484e+01 520
P 4000 4.862297e-02 -1.392835e+00 1.441458e+00 1.612399e+02 2.480776e+01 625
LD 4311 2.750372e-02 -3.639247e-01 3.914285e-01 1.682647e+02 2.562034e+01 648
P 5000 1.488184e-02 -3.765466e-01 3.914285e-01 1.686278e+02 2.965111e+01 626
LD 5521 1.184656e-02 -8.883781e-02 1.006844e-01 1.687472e+02 3.271759e+01 625
P 6000 1.057020e-02 -9.011417e-02 1.006844e-01 1.690627e+02 3.548979e+01 625
P 7000 9.610191e-03 -9.107418e-02 1.006844e-01 1.692628e+02 4.135581e+01 625
LD 7455 9.459506e-03 -1.662258e-02 2.608208e-02 1.693719e+02 4.401557e+01 625
P 8000 9.366945e-03 -1.671514e-02 2.608208e-02 1.696977e+02 4.714266e+01 625
P 9000 9.300757e-03 -1.678132e-02 2.608208e-02 1.698969e+02 5.297331e+01 625
LD 9578 9.287270e-03 2.510164e-03 6.777106e-03 1.700244e+02 5.633310e+01 625
P 10000 9.282033e-03 2.504927e-03 6.777106e-03 1.703558e+02 5.870067e+01 625
Last 10001 9.282014e-03 3.724653e-03 5.557361e-03 1.706567e+02 5.860302e+01 625
----------------------------------------------------------------------------------------------------------------
PP 10001 9.282014e-03 3.724653e-03 5.557361e-03 1.709555e+02 5.850059e+01 625
----------------------------------------------------------------------------------------------------------------
171.202053 seconds (3.30 G allocations: 77.213 GiB, 5.04% gc time, 0.03% compilation time)In this small example, the acceleration is quite minimal, but as soon as one of the following conditions is met, significant speedups (factor ten at least) can be expected:
- quite expensive scalar product between atoms, for instance, due to a high dimension (say, more than 10000),
- high number of atoms in the active set (say, more than 1000),
- high number of iterations (say, more than 100000), spending most of the time redistributing the weights in the active set.
Dimension reduction via symmetrization
Permutation of the tensor axes
It is easy to see that our specific instance remains invariant under permutation of the dimensions of the tensor. This means that all computations can be performed in the symmetric subspace, which leads to an important speedup, owing to the reduced dimension (hence reduced size of the final active set and reduced number of iterations).
The way to operate this in the FrankWolfe package is to use a symmetrized LMO, which basically does the following:
- symmetrize the gradient, which is not necessary here as the gradient remains symmetric throughout the algorithm,
- call the standard LMO,
- symmetrize its output, which amounts to averaging over its orbit with respect to the group considered (here the symmetric group permuting the dimensions of the tensor).
function reynolds_permutedims(atom::Array{T, N}, lmo::BellCorrelationsLMO{T}) where {T <: Number, N}
res = zeros(T, size(atom))
for per in Combinatorics.permutations(1:N)
res .+= permutedims(atom, per)
end
res ./= factorial(N)
return res
endNote that the second argument lmo is not used here but could in principle be exploited to obtain a very small speedup by precomputing and storing Combinatorics.permutations(1:N) in a dedicated field of our custom LMO.
lmo_permutedims = FrankWolfe.SymmetricLMO(lmo_naive, reynolds_permutedims)
asq_permutedims = FrankWolfe.ActiveSetQuadratic([(one(T), x0)], LinearAlgebra.I, -p)
@time FrankWolfe.blended_pairwise_conditional_gradient(f, grad!, lmo_permutedims, asq_permutedims; verbose, lazy=true, line_search=FrankWolfe.Shortstep(one(T)), max_iteration)
Blended Pairwise Conditional Gradient Algorithm.
MEMORY_MODE: FrankWolfe.InplaceEmphasis() STEPSIZE: Shortstep EPSILON: 1.0e-7 MAXITERATION: 10000 TYPE: Float64
GRADIENTTYPE: Array{Float64, 5} LAZY: true lazy_tolerance: 2.0
LMO: FrankWolfe.SymmetricLMO{Main.BellCorrelationsLMO{Float64}, typeof(Main.reynolds_permutedims), FrankWolfe.var"#10#11"}
----------------------------------------------------------------------------------------------------------------
Type Iteration Primal Dual Dual Gap Time It/sec #ActiveSet
----------------------------------------------------------------------------------------------------------------
I 1 4.132813e+01 -4.029553e+01 8.162365e+01 0.000000e+00 Inf 1
LD 11 1.153403e+01 -2.734161e+01 3.887563e+01 1.798675e+00 6.115614e+00 9
LD 31 3.285725e+00 -1.376230e+01 1.704803e+01 3.073525e+00 1.008614e+01 13
LD 50 1.646307e+00 -6.869795e+00 8.516102e+00 3.823962e+00 1.307544e+01 15
LD 93 7.245404e-01 -2.896319e+00 3.620859e+00 5.806974e+00 1.601523e+01 20
LD 172 2.996542e-01 -1.378052e+00 1.677706e+00 7.568174e+00 2.272675e+01 24
LD 329 7.534658e-02 -6.003412e-01 6.756878e-01 9.619555e+00 3.420116e+01 29
LD 654 2.017129e-02 -2.325514e-01 2.527227e-01 1.092043e+01 5.988776e+01 32
LD 918 1.120386e-02 -5.331182e-02 6.451568e-02 1.121052e+01 8.188734e+01 26
P 1000 1.035776e-02 -5.415791e-02 6.451568e-02 1.150581e+01 8.691258e+01 26
LD 1324 9.405322e-03 -6.061341e-03 1.546666e-02 1.156681e+01 1.144655e+02 26
LD 1763 9.281060e-03 6.068677e-03 3.212383e-03 1.185873e+01 1.486669e+02 26
P 2000 9.275508e-03 6.063125e-03 3.212383e-03 1.214353e+01 1.646968e+02 26
LD 2213 9.274505e-03 8.598448e-03 6.760579e-04 1.218141e+01 1.816703e+02 26
LD 2646 9.274227e-03 9.121418e-03 1.528092e-04 1.249676e+01 2.117348e+02 26
P 3000 9.274215e-03 9.121406e-03 1.528092e-04 1.281493e+01 2.341020e+02 26
LD 3077 9.274215e-03 9.237360e-03 3.685452e-05 1.282990e+01 2.398305e+02 26
LD 3514 9.274214e-03 9.265220e-03 8.994042e-06 1.314787e+01 2.672676e+02 26
LD 3915 9.274214e-03 9.272267e-03 1.947318e-06 1.345815e+01 2.909017e+02 26
P 4000 9.274214e-03 9.272267e-03 1.947318e-06 1.373295e+01 2.912703e+02 26
LD 4376 9.274214e-03 9.273758e-03 4.562791e-07 1.379790e+01 3.171498e+02 26
LD 4777 9.274214e-03 9.274102e-03 1.122895e-07 1.411101e+01 3.385301e+02 26
P 5000 9.274214e-03 9.274102e-03 1.122895e-07 1.439181e+01 3.474198e+02 26
LD 5151 9.274214e-03 9.274188e-03 2.584730e-08 1.441883e+01 3.572412e+02 26
Last 5151 9.274214e-03 9.274188e-03 2.584730e-08 1.492404e+01 3.451478e+02 26
----------------------------------------------------------------------------------------------------------------
PP 5151 9.274214e-03 9.274188e-03 2.584730e-08 1.516964e+01 3.395599e+02 26
----------------------------------------------------------------------------------------------------------------
15.418026 seconds (274.86 M allocations: 6.612 GiB, 4.65% gc time, 0.38% compilation time)Now, convergence is reached within 10000 iterations, and the size of the final active set is considerably smaller than before, thanks to the reduced dimension.
Uniqueness pattern
In this specific case, there is a bigger symmetry group that we can exploit. Its action roughly allows us to work in the subspace respecting the structure of the objective point p, that is, to average over tensor entries that have the same value in p. Although quite general, this kind of symmetry is not always applicable, and great care has to be taken when using it, in particular, to ensure that there exists a suitable group action whose Reynolds operator corresponds to this averaging procedure. In our current case, the theoretical study enabling this further symmetrization can be found here.
function build_reynolds_unique(p::Array{T, N}) where {T <: Number, N}
ptol = round.(p; digits=8)
ptol[ptol .== zero(T)] .= zero(T) # transform -0.0 into 0.0 as isequal(0.0, -0.0) is false
uniquetol = unique(ptol[:])
indices = [ptol .== u for u in uniquetol]
return function(A::Array{T, N}, lmo)
res = zeros(T, size(A))
for ind in indices
@view(res[ind]) .= sum(A[ind]) / sum(ind) # average over ind
end
return res
end
end
lmo_unique = FrankWolfe.SymmetricLMO(lmo_naive, build_reynolds_unique(p))
asq_unique = FrankWolfe.ActiveSetQuadratic([(one(T), x0)], LinearAlgebra.I, -p)
@time FrankWolfe.blended_pairwise_conditional_gradient(f, grad!, lmo_unique, asq_unique; verbose, lazy=true, line_search=FrankWolfe.Shortstep(one(T)), max_iteration)
Blended Pairwise Conditional Gradient Algorithm.
MEMORY_MODE: FrankWolfe.InplaceEmphasis() STEPSIZE: Shortstep EPSILON: 1.0e-7 MAXITERATION: 10000 TYPE: Float64
GRADIENTTYPE: Array{Float64, 5} LAZY: true lazy_tolerance: 2.0
LMO: FrankWolfe.SymmetricLMO{Main.BellCorrelationsLMO{Float64}, Main.var"#43#45"{5, Float64, Vector{BitArray{5}}}, FrankWolfe.var"#10#11"}
----------------------------------------------------------------------------------------------------------------
Type Iteration Primal Dual Dual Gap Time It/sec #ActiveSet
----------------------------------------------------------------------------------------------------------------
I 1 4.132813e+01 -4.029553e+01 8.162365e+01 0.000000e+00 Inf 1
LD 4 2.991558e+00 -4.896575e+00 7.888134e+00 2.710608e-01 1.475684e+01 3
LD 7 4.010409e-01 -3.282121e+00 3.683161e+00 5.235709e-01 1.336973e+01 3
LD 15 1.802237e-01 -4.675134e-01 6.477370e-01 1.540416e+00 9.737627e+00 5
LD 32 2.002807e-02 -1.556256e-01 1.756537e-01 2.267057e+00 1.411522e+01 5
LD 72 1.368120e-02 -4.123301e-02 5.491421e-02 2.497637e+00 2.882725e+01 3
LD 130 9.506383e-03 -9.057423e-03 1.856381e-02 3.008074e+00 4.321702e+01 4
LD 145 9.308486e-03 1.142443e-03 8.166042e-03 3.226554e+00 4.493959e+01 4
LD 160 9.278954e-03 6.339760e-03 2.939194e-03 3.464956e+00 4.617663e+01 4
LD 178 9.275038e-03 8.020945e-03 1.254093e-03 3.743968e+00 4.754314e+01 4
LD 193 9.274328e-03 8.822599e-03 4.517287e-04 3.955830e+00 4.878875e+01 4
LD 211 9.274234e-03 9.080960e-03 1.932742e-04 4.197688e+00 5.026577e+01 4
LD 226 9.274217e-03 9.204592e-03 6.962503e-05 4.439777e+00 5.090346e+01 4
LD 244 9.274214e-03 9.244275e-03 2.993943e-05 4.691582e+00 5.200804e+01 4
LD 259 9.274214e-03 9.263418e-03 1.079580e-05 4.947783e+00 5.234667e+01 4
LD 277 9.274214e-03 9.269571e-03 4.642816e-06 5.186889e+00 5.340388e+01 4
LD 292 9.274214e-03 9.272535e-03 1.678918e-06 5.444781e+00 5.362934e+01 4
LD 307 9.274214e-03 9.273611e-03 6.027276e-07 5.679246e+00 5.405647e+01 4
LD 323 9.274214e-03 9.273974e-03 2.398615e-07 5.920923e+00 5.455230e+01 4
LD 338 9.274214e-03 9.274110e-03 1.036546e-07 6.186729e+00 5.463307e+01 4
LD 351 9.274214e-03 9.274168e-03 4.575080e-08 6.412647e+00 5.473558e+01 4
Last 351 9.274214e-03 9.274168e-03 4.575080e-08 6.923227e+00 5.069890e+01 4
----------------------------------------------------------------------------------------------------------------
PP 351 9.274214e-03 9.274168e-03 4.575080e-08 7.140386e+00 4.915700e+01 4
----------------------------------------------------------------------------------------------------------------
7.351973 seconds (142.16 M allocations: 3.339 GiB, 4.95% gc time, 0.72% compilation time)Reduction of the memory footprint of the iterate
In the previous run, the dimension reduction is mathematically exploited to accelerate the algorithm, but it is not used to effectively work in a subspace of reduced dimension. Indeed, the iterate, although symmetric, was still a full tensor. As a last example of the speedup obtainable through symmetry reduction, we show how to map the computations into a space whose physical dimension is also reduced during the algorithm. This makes all in-place operations marginally faster, which can lead, in bigger instances, to significant accelerations, especially for active set based algorithms in the regime where many lazy iterations are performed. We refer to the example symmetric.jl for a small benchmark with symmetric matrices.
function build_reduce_inflate(p::Array{T, N}) where {T <: Number, N}
ptol = round.(p; digits=8)
ptol[ptol .== zero(T)] .= zero(T) # transform -0.0 into 0.0 as isequal(0.0, -0.0) is false
uniquetol = unique(ptol[:])
dim = length(uniquetol) # reduced dimension
indices = [ptol .== u for u in uniquetol]
mul = [sum(ind) for ind in indices] # multiplicities, used to have matching scalar products
sqmul = sqrt.(mul) # precomputed for speed
return function(A::Array{T, N}, lmo)
vec = zeros(T, dim)
for (i, ind) in enumerate(indices)
vec[i] = sum(A[ind]) / sqmul[i]
end
return FrankWolfe.SymmetricArray(A, vec)
end, function(x::FrankWolfe.SymmetricArray, lmo)
for (i, ind) in enumerate(indices)
@view(x.data[ind]) .= x.vec[i] / sqmul[i]
end
return x.data
end
end
reduce, inflate = build_reduce_inflate(p)
p_reduce = reduce(p, nothing)
x0_reduce = reduce(x0, nothing)
f_reduce = let p_reduce = p_reduce, normp2 = normp2
x -> LinearAlgebra.dot(x, x) / 2 - LinearAlgebra.dot(p_reduce, x) + normp2
end
grad_reduce! = let p_reduce = p_reduce
(storage, x) -> begin
@inbounds for i in eachindex(x)
storage[i] = x[i] - p_reduce[i]
end
end
endNote that the objective function and its gradient have to be explicitly rewritten. In this simple example, their shape remains unchanged, but in general this may need some reformulation, which falls to the user.
lmo_reduce = FrankWolfe.SymmetricLMO(lmo_naive, reduce, inflate)
asq_reduce = FrankWolfe.ActiveSetQuadratic([(one(T), x0_reduce)], LinearAlgebra.I, -p_reduce)
@time FrankWolfe.blended_pairwise_conditional_gradient(f_reduce, grad_reduce!, lmo_reduce, asq_reduce; verbose, lazy=true, line_search=FrankWolfe.Shortstep(one(T)), max_iteration)
Blended Pairwise Conditional Gradient Algorithm.
MEMORY_MODE: FrankWolfe.InplaceEmphasis() STEPSIZE: Shortstep EPSILON: 1.0e-7 MAXITERATION: 10000 TYPE: Float64
GRADIENTTYPE: FrankWolfe.SymmetricArray{false, Float64, Array{Float64, 5}, Vector{Float64}} LAZY: true lazy_tolerance: 2.0
LMO: FrankWolfe.SymmetricLMO{Main.BellCorrelationsLMO{Float64}, Main.var"#48#52"{5, Float64, Vector{Float64}, Vector{BitArray{5}}, Int64}, Main.var"#49#53"{Vector{Float64}, Vector{BitArray{5}}}}
----------------------------------------------------------------------------------------------------------------
Type Iteration Primal Dual Dual Gap Time It/sec #ActiveSet
----------------------------------------------------------------------------------------------------------------
I 1 4.132813e+01 -4.029553e+01 8.162365e+01 0.000000e+00 Inf 1
LD 4 2.991558e+00 -4.896575e+00 7.888134e+00 3.038019e-01 1.316648e+01 3
LD 7 4.010409e-01 -3.282121e+00 3.683161e+00 5.546945e-01 1.261956e+01 3
LD 16 1.663895e-01 -3.019507e-01 4.683402e-01 1.524055e+00 1.049831e+01 5
LD 30 1.909390e-02 -6.502115e-02 8.411505e-02 2.024959e+00 1.481511e+01 4
LD 131 9.738911e-03 -1.106636e-02 2.080527e-02 2.578751e+00 5.079980e+01 4
LD 155 9.318517e-03 5.789313e-03 3.529204e-03 2.830931e+00 5.475231e+01 4
LD 207 9.275282e-03 8.691635e-03 5.836471e-04 3.073480e+00 6.735036e+01 4
LD 263 9.274251e-03 9.174862e-03 9.938911e-05 3.309358e+00 7.947161e+01 4
LD 320 9.274215e-03 9.255732e-03 1.848329e-05 3.554988e+00 9.001438e+01 4
LD 376 9.274214e-03 9.271193e-03 3.021432e-06 3.796650e+00 9.903469e+01 4
LD 427 9.274214e-03 9.273665e-03 5.489460e-07 4.031752e+00 1.059093e+02 4
LD 483 9.274214e-03 9.274124e-03 9.004990e-08 4.284591e+00 1.127295e+02 4
Last 483 9.274214e-03 9.274124e-03 9.004990e-08 4.768426e+00 1.012913e+02 4
----------------------------------------------------------------------------------------------------------------
PP 483 9.274214e-03 9.274124e-03 9.004990e-08 5.028572e+00 9.605113e+01 4
----------------------------------------------------------------------------------------------------------------
5.270069 seconds (99.55 M allocations: 2.339 GiB, 5.57% gc time, 1.14% compilation time)This page was generated using Literate.jl.