The interior point solver can sometimes be sped up dramatically by exploiting the simultanious structure of $Q$,$A$ and $G$. Here we will consider the simple problem:
$$\mbox{minimize}\quad\frac{1}{2}x^{T}Qx-b^{T}x,\qquad\mbox{s.t.}\quad x\geq0.$$The solver can be called with no special parameters,
using ConicIP
n = 1000
Q = sparse(randn(n,n));
Q = Q'*Q;
c = ones(n,1);
A = speye(n);
b = zeros(n,1);
𝐾 = [("R",n)];
@time conicIP( Q , c , A , b , 𝐾 , verbose = true);
The speed of the solver is reasonable, as the deault solver exploits the sparsity of the constraint matrix. We can do better, however.
To speed up the solver, we require a function solve3x3gen which returns a function solving the KKT system
Where $F$ is a Block matrix with blocks corrosponding to the cones specified in 𝐾,
The operation $vec$ is the vectorization operator on symmetric matrices $$ \mbox{vec}(U)=(U_{11},\sqrt{2}U_{21},\dots,\sqrt{2}U_{p1},U_{22},\sqrt{2}U_{32},\dots,\sqrt{2}U_{p2},\dots,U_{p-1,p-1},\sqrt{2}U_{p,p-1},U_{pp}) $$
and $\mbox{mat}$ is the inverse transformation back to a symmetric matrix.
The matrix $\mbox{diag}(u)$ is represented by type Diag, $αI+uu$ is represented by type SymWoodbury, and $A_{U_i}$ is represented by the type VecCongurance.
In this example, since we have no linear constraints, $G$ is empty, and our KKT system is $$ \left(\begin{array}{cc} Q & I\\ I & -F^{T}F \end{array}\right)\left(\begin{array}{c} a\\ b \end{array}\right)=\left(\begin{array}{c} x\\ y \end{array}\right) $$
The system can be solved by pivoting on the second block, as follows: $$ \big(Q+(F^TF)^{-1} \big)\,a=x+(F^TF)^{-1}y,\qquad b=(F^TF)^{-1}(a-y) $$
Because we only have polyhedral constraints, $F^{-2}$ is a diagonal matrix, thus the first equation is a diagonal perturbation to $Q$ which can be solved via a Cholesky Factorization. Pivoting allows us to solve a 1000x1000 system rather than a 2000x2000 system (albeit with a some sparsity structure).
function kktsolver(Q, A, G)
function solve3x3gen(F, F⁻¹)
invFᵀF = inv(F'F)
QpD⁻¹ = cholfact(Q + spdiagm( (F[1].diag).^(-2) ))
function solve3x3(x, z, y)
a = QpD⁻¹\(x + A'*(invFᵀF*y))
b = invFᵀF*(y - A*a)
c = zeros(0,1) # empty vector.
return(a, c, b)
end
end
end
@time sol = conicIP( Q , c , A , b , 𝐾 ,
kktsolver = kktsolver;
verbose = true);
This results in a 5-fold improvement in speed, and a dramatic drop in memory usage!
This pattern of pivoting on the third block happens often enough that we have encapsulated it in the convenience function pivot, which transforms a $2x2$ solver of the system
into a $3x3$ solver. This is illustrated below
function kktsolver2x2(Q, A, G)
function solve3x3gen(F, F⁻¹)
QpD⁻¹ = cholfact(Q + spdiagm( (F[1].diag).^(-2) ))
return (y, x) -> (QpD⁻¹\y, zeros(0,1))
end
end
@time sol = conicIP( Q , c , A , b , 𝐾 ,
kktsolver = pivot(kktsolver2x2);
verbose = true);
And as a bonus we even get an extra boost in speed!
HTML(readall(open("style.css")))