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,
InΒ [17]:
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 kktsolver which returns a function solve3x3gen which, in turn returns a function solving the KKT system
$$\left(\begin{array}{ccc}
Q & G^{T} & A^{T}\\
G\\
A & & -F^{T}F
\end{array}\right)
\left(\begin{array}{c}
a\\
c\\
b
\end{array}\right) = \left(\begin{array}{c}
x\\
z\\
y
\end{array}\right)$$
Where $F$ is a Block matrix with blocks corrosponding to the cones specified in πΎ,
$$F[i]=\begin{cases}
\mbox{diag}(u_i) & \mbox{if }K_{i}\mbox{ is "R"}\\
Ξ±_iJ+u_i u_i^T & \mbox{if }K_{i}\mbox{ is "Q"} \quad
\mbox{ for some } u_i, \alpha_i \mbox{ $J$ is the hyperbolic identity }\\
A_{U_i} & \mbox{if }K_i\mbox{ is "S" }\, \quad \mbox{where } A_{U_i}x=\mbox{vec}(U_i^{T}\mbox{mat}(x)U_i) \mbox{ for all $x$}.
\end{cases}
$$
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).
InΒ [11]:
function kktsolver(Q, A, G, cone_dims)
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)
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
InΒ [14]:
function kktsolver2x2(Q, A, G, cone_dims)
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!
InΒ [16]:
HTML(readall(open("style.css")))
Out[16]: