| Status | Coverage |
|---|---|
 |
 |
Goldfarb and Idnani quadratic programming solver in native Julia.
This is a port of Berwin A. Turlach quadprog to the Julia language.
Julia has access to several high quality convex optimizer, which are able to solve quadratic programming problems, such as, for example, COSMO. But those tend to be not optimized for pure basic quadratic problems. It supports dense as well as sparse matrices.
This solver can be 100 times faster on some problems (and without any specific code level optimization). Being native Julia, it also supports multiple precision big floats.
The solveQPcompact function implements the dual method of Goldfarb and Idnani (1982, 1983) for solving quadratic programming problems of the form min(-d^T b + 1/2 b^T D b) with the constraints A^T b >= b_0.
- Create a sparse matrix with A, then transform it to the quadprog sparse format via
amat, aind = convertSparse(A)
- and solve your problem
sol, lagr, crval, iact, nact, iter = solveQPcompact(D, d, amat, aind, b_0)
The solution is in sol. It is possible through factorized=true to skip the factorization step for further speed-up. In this case, the matrix D should contain its inverse square root. The parameter meq=0 allows to specify the number of equalities (as in quadprog R package).
The solveQP function is the same algorithm for dense matrices:
- solve your problem using
sol, lagr, crval, iact, nact, iter = solveQP(D, d, A, b_0)
A good example that shows the performance of the solver is to find the convex set closest to a set of financial option prices. This ensures the so-called arbitrage-free property of the option prices, see An arbitrage-free interpolation of class C2 for option prices.
In terms of code, we fill up the constraints in the sparse matrix G, and associated vector h. The D matrix is simply a diagonal matrix of squared weights, here we input the inverse square root directly. And then we call the GoldfarbIrdaniSolver
using GoldfarbIdnaniSolver, SparseArrays, LinearAlgebra
function filterConvexCallPrices(
strikes::Vector{T},
callPrices::Vector{T}, #undiscounted!
weights::Vector{T},
forward::T;
tol = 1e-8
)::Tuple{Vector{T},Vector{T},Vector{T}} where {T}
n = length(callPrices)
G = spzeros(T, 2 * n, n)
h = zeros(T, 2 * n)
for i = 2:n-1
dym = (strikes[i] - strikes[i-1])
dy = (strikes[i+1] - strikes[i])
G[i, i-1] = -1 / dym
G[i, i] = 1 / dym + 1 / dy
G[i, i+1] = -1 / dy
end
G[1, 1] = 1 / (strikes[2] - strikes[1])
G[1, 2] = -G[1, 1]
G[n, n] = 1 / (strikes[n] - strikes[n-1])
G[n, n-1] = -G[n, n]
for i = 1:n
h[i] = -tol
G[n+i, i] = -1
h[n+i] = -max(forward - strikes[i], 0) - tol
end
h[1] = 1 - tol
strikesf = strikes
#call solver
amat, aind = convertSparse(copy(-G'))
factorized = true
dmat = diagm(1.0 ./ weights)
dvec = @. callPrices * weights^2
nEqualities = (strikes[1] == 0) ? 1 : 0
pricesf, lagr, crval, iact, nact, iter = solveQPcompact(dmat, dvec, amat, aind, -h, meq=nEqualities, factorized = true)
return strikesf, pricesf, weights
endIn terms of Black-Scholes implied volatility, the convex filtering is very visible, and the implied volatility look smoother (linear interpolation is used in the figure):
For the COSMO solver, the code below the comment would read
W = spdiagm(weights)
z = Variable(n)
problem = minimize(square(norm(W * (z - callPrices))), G * z <= h)
#solve!(problem, () -> SCS.Optimizer(verbose = 0))
Convex.solve!(problem, () -> COSMO.Optimizer(verbose = false, eps_rel = 1e-8, eps_abs = 1e-8))
pricesf = Convex.evaluate(z)
return strikesf, pricesf, weightsWe can benchmark the solvers as follows
using BenchmarkTools, StatsBase
strikes = Float64.([20, 25, 50, 55, 75, 100, 120, 125, 140, 150, 160, 175, 180, 195, 200, 210, 230, 240, 250, 255, 260, 270, 275, 280, 285, 290, 300, 310, 315, 320, 325, 330, 335, 340, 350, 360, 370, 380, 390, 400, 410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 550, 580, 590, 600, 650, 670, 680, 690, 700])
prices = [337.9478782712897, 333.2264008151271, 310.49584376491333, 306.5789869909477, 288.3326474521462, 266.1436620149496, 249.44466994191146, 244.47045906634085, 232.02358238485036, 223.57092985061374, 215.47932505838847, 203.49665261662705, 199.45085022051444, 187.64870164975878, 183.73184487579312, 176.72338330960227, 161.59752782675685, 154.84695750485992, 147.27113520122242, 144.51478902657934, 141.11371474120148, 134.64682478802777, 131.65837649352014, 128.3346695814304, 125.60411253121669, 122.28040561912707, 115.9424612881003, 110.99712967626056, 108.21499437718802, 105.35549170482741, 102.75388027676055, 100.15226884869381, 97.60223566948576, 95.07799161470709, 90.23581650058495, 85.61989310557237, 80.97785070828226, 76.85159079958001, 72.72533089087777, 67.85118637372317, 64.96280443763155, 61.352327017517055, 58.025529966125816, 54.95662415902851, 51.810350978643, 49.07670378912775, 46.34305659961247, 43.11941604593882, 41.05628609158771, 38.735264892942695, 36.285298072150695, 34.8668962285343, 28.935397609774668, 24.370722585772782, 23.003898991015205, 21.27602765424615, 16.324515763803383, 14.545066178175517, 13.719814196435081, 13.049296961270962, 12.22404497953049]
weights = [0.18348804852853492, 0.14552754066163864, 0.06087460548064749, 0.05109612140267557, 0.03629649739986478, 0.025613847407177506, 0.020011216481149133, 0.019345036201403803, 0.01663883576809264, 0.015257762783202483, 0.014011545380885547, 0.012468184742920784, 0.012034729348791627, 0.010873040937928326, 0.010535157435308888, 0.009879102038500473, 0.008855540796164947, 0.00841084879707799, 0.008032374522333945, 0.007845822148042107, 0.00767850819280654, 0.00737017036267451, 0.007228761921893696, 0.007096391696183667, 0.00697143229125802, 0.006853569353455795, 0.006637272585167783, 0.0064509208672822085, 0.006365169495070308, 0.006284168069410054, 0.006209646991746033, 0.0061399098141406445, 0.0060750536990634395, 0.006014709985843373, 0.0059077251648834515, 0.0058176064495057645, 0.005741621638057057, 0.005682166059320373, 0.005635273659488148, 0.005598197720478721, 0.00557969512744412, 0.005570043788415912, 0.005571719928608071, 0.005583741283162423, 0.005606897345237971, 0.0056386186100743, 0.005680979257062336, 0.005742146253326264, 0.005798784362859748, 0.005870563857815834, 0.005958641403270033, 0.006029457615736634, 0.0063683433148160615, 0.0067726154457943945, 0.006927129393879053, 0.007136508419763993, 0.008060162361264996, 0.00852974277055327, 0.008785608144514201, 0.009023417426706495, 0.009332153986076881]
forward = 356.73063159822254
tte = 1.5917808219178082
strikesf, pricesf = filterConvexCallPrices(strikes, prices, weights, forward, tol = 1e-6)
rmse = rmsd(pricesf .* weights, prices .* weights)
@benchmark strikesf, pricesf = filterConvexCallPrices(strikes, prices, weights, forward, tol = 1e-6)The benchmark results are:
| Solver | RMSE | Time (ms) |
|---|---|---|
| GoldfarbIdnaniSolver | 0.0031130349998388157 | 0.205 |
| COSMO | 0.0031130309597602370 | 21.485 |
| SCS | 0.0031130429769795193 | 8.381 |
| ECOS | 0.0031131478052242090 | 2.430 |