I aim to develop a package for solving the standard quadratic optimization problem which maximizes (minimizes) a quadratic form over the unit simplex. That is \max x^T Q x : \sum_{i \in [1:n]} x_i = 1, x\in \mathbb{R}_+^n where Q \in \mathbb{R}^{n \times n} and symmetric. I use an algorithm called infection immunization dynamics which already works quite well. Without going into details on that algorithm, I wonder you have some tips how I can improve the performance of my function solve_StQP() in the MWE below. For n=400 I obtain
0.030732 seconds (17.69 k allocations: 16.187 MiB)
and for n=4000
1.220787 seconds (189.63 k allocations: 1.221 GiB, 4.55% gc time)
Thank you very much for your inupts. A (quite long) MWE is here:
using LinearAlgebra
using Random
# set parameters
ndims = 400 # dimensions
mtype = "PSD" # matrix type "PSD" (positive semidefinite), "NSD" (negative semidefinite), "IND" (indefinite)
sense = :minimize # objective sense (:minimize or :maximize)
Random.seed!(0) # random seed
# THIS SHOULD BE IMPROVED
function solve_StQP(Q, sense::Symbol; tol::Float64=0.000001, x=ones(Float64,size(Q)[1]) ./ size(Q)[1], verbose=1, epssize=5)
if sense == :minimize
Q = -Q
end
if verbose==1
println("starting point: ", x)
end
#start vectors and matrices
Q_x::Vector{Float64} = Q * x # used to avoid excessive number of multiplications later in loop (this vector will be updated instead)
p::Float64 = x' * Q_x # current payoff x'Qx
t::Vector{Float64} = Q_x .- p # tempvector = (e^i - x | x) for all i
# define other start values
n::Int = size(Q)[1] # dimension of decision vector
epsilon::Float64 = Inf # measures violation of Nash condition. If epsilon = 0, we have a Nash equilibrium.
obj::Float64 = 0.0 # objective value
prec::Float64 = epssize * eps(Float64)# obj epsilon surrounding around zero (to avoid numerical instabilities)
i::Int = 0 # index of infective pure / or co-strategy
tempval::Float64 = 0.0 # temp value used for identifying the largest abs(t[i]), thus, best direction i
delta::Float64 = 1.0 # fraction of mutants that are infective
niterations::Int = 0 # number of iterations
num::Float64 = 0.0 # dummy numerator
denom::Float64 = 1.0 # dummy denominator
k::Float64 = 0.0 # factor needed later
#errorcount::Int = 0
#start loop
while epsilon > tol
niterations += 1
# reset values
i = 0 # index of pure / or co strategy
tempval = 0 # temp value used for identifying the largest abs(t[i]), thus, direction (or pure/co-strategy) i
# find most infective pure (co-)strategy i
for j in eachindex(x)
if (t[j] > prec) && (tempval < t[j] - prec)
i = j
tempval = t[i]
elseif (t[j] < -prec) && (x[j] > prec) && (tempval < abs(t[j] - prec))
i = j
tempval = abs(t[i])
else
# do not update i
end
end
if i == 0 # break loop, if no strategy was found (we are already optimal)
break
end
# select (infective) strategy, update delta, and inject mutants
delta = 1.0
if t[i] > prec # selection: pure strategy i is infective
num = -t[i] # numerator
denom = Q[i,i] - 2 * Q_x[i] + p # denominator
if denom < -prec
delta = min(num / denom, 1) # update delta
end
Q_x = delta * ( Q[:,i] - Q_x ) + Q_x
x = delta * ( LinearAlgebra.I[1:n,i] - x ) + x # inject mutants
elseif (t[i] < -prec) && (x[i] > prec) # selection: co-strategy of i is infective
k = x[i] / ( x[i]-1 )
num = k * (-t[i]) # numerator
denom = k^2 * ( Q[i,i] - 2 * Q_x[i] + p ) # denominator
if denom < -prec
delta = min(num / denom, 1) # update delta
end
Q_x = k * delta * ( Q[:,i] - Q_x ) + Q_x
x = delta * k * ( LinearAlgebra.I[1:n,i] -x ) +x # inject mutants
else # selection: population x is immune against all mutants (=nash equilibrium)
Q_x = Q * x
end
# update current payoff, and temp vector
p = x' * Q_x # current payoff x'Qx
t = Q_x .- p # tempvector
# update violation of Nash condition (used as stopping criteria).
epsilon=0
for i in 1:n
epsilon = epsilon + ( min( x[i], -t[i] ) )^2
end
# optional: print current payoff, and Nash condition violation
if verbose == 1
cur_obj = 0.0
if sense == :maximize
cur_obj = round(p, digits=4)
else
cur_obj = -round(p, digits=4)
end
println("obj: $cur_obj \t\tepsilon: $epsilon")
end
end# while
# set objective
if sense == :maximize
obj = p
else
obj = -p
end
# round output vector
x = round.(x, digits=6)
# @info "# errorcount / niterations: $errorcount / $niterations"
return obj, x
end
# DUMMY FUNCTIONS
function createSymmetricMatrix(ndims::Int, type::String="IND")
A = LinearAlgebra.Symmetric(rand(ndims,ndims)) * 10
if type == "PSD"
A = A*A'
end
if type == "NSD"
A = -A*A'
end
return LinearAlgebra.Symmetric(A)
end
function createRandomStartingPoint(Q)
n = size(Q)[1]
x = rand(n)
x = x ./ sum(x)
return x
end
function main(ndims, mtype, sense)
Q = createSymmetricMatrix(ndims, mtype)
xstart = createRandomStartingPoint(Q) #optional random starting point
# the performance of this function should be improved
@time obj, x = solve_StQP(Q, sense; verbose=1)
#@time obj, x = solve_StQP(Q, sense; x=xstart) #optional with different starting point
# print results
println("obj: $obj")
#println("x: $x")
return
end
main(ndims, mtype, sense)
