# KKT Callbacks with IntPoint¶

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 [4]:
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);

 > INTERIOR POINT SOLVER v0.7 (July 2016)

Iter  Mu          prFeas      duFeas      muFeas      refine
1  1.7653e+00  1.3071e+00  4.1667e+01  1.5633e+00  1
2  3.8817e-01  2.7260e-01  8.6896e+00  4.3973e-01  1
3  1.3674e-01  7.3342e-02  2.3379e+00  2.9856e-01  1
4  4.0684e-02  1.6830e-02  5.3649e-01  1.4986e-01  1
5  1.0436e-02  2.5143e-03  8.0146e-02  6.5562e-02  1
6  1.3460e-03  1.1343e-14  6.1191e-17  8.5276e-03  1
7  1.8750e-04  8.4381e-15  3.9569e-17  1.1522e-03  1
8  2.3185e-05  7.7170e-15  4.1380e-17  1.4949e-04  1
9  2.1053e-06  8.3087e-15  4.1707e-17  1.5936e-05  1
10  9.6541e-08  8.0378e-15  4.3236e-17  1.3858e-06  1

> EXIT -- Below Tolerance!

4.550797 seconds (29.30 k allocations: 1.328 GB, 10.45% gc time)


The speed of the solver is reasonable, as the deault solver exploits the sparsity of the constraint matrix. We can do better, however.

## KKT Callbacks¶

To speed up the solver, we require a function solve3x3gen which 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"}\\ α_iI+u_i u_i^T & \mbox{if }K_{i}\mbox{ is "Q"} \quad \mbox{ for some } u_i, \alpha_i\\ 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 [12]:
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);

 > INTERIOR POINT SOLVER v0.7 (July 2016)

Iter  Mu          prFeas      duFeas      muFeas      refine
1  1.7653e+00  1.3071e+00  4.1667e+01  1.5633e+00  1
2  3.8817e-01  2.7260e-01  8.6896e+00  4.3973e-01  1
3  1.3674e-01  7.3342e-02  2.3379e+00  2.9856e-01  1
4  4.0684e-02  1.6830e-02  5.3649e-01  1.4986e-01  1
5  1.0436e-02  2.5143e-03  8.0146e-02  6.5562e-02  1
6  1.3460e-03  1.2157e-14  5.1492e-17  8.5276e-03  1
7  1.8750e-04  9.1014e-15  3.9022e-17  1.1522e-03  1
8  2.3185e-05  8.4269e-15  3.5807e-17  1.4949e-04  1
9  2.1053e-06  8.0958e-15  3.7206e-17  1.5936e-05  1
10  9.6541e-08  8.1248e-15  3.6720e-17  1.3858e-06  1

> EXIT -- Below Tolerance!

3.213126 seconds (32.19 k allocations: 1.035 GB, 29.09% gc time)


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

$$\left(\begin{array}{cc} Q+A^{T}(F^{T}F)^{-1}A & G\\ G & 0 \end{array}\right)\left(\begin{array}{c} a\\ b \end{array}\right) = \left(\begin{array}{c} x\\ y \end{array}\right)$$

into a $3x3$ solver. This is illustrated below

In [8]:
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);

 > INTERIOR POINT SOLVER v0.7 (July 2016)

Iter  Mu          prFeas      duFeas      muFeas      refine
1  1.7653e+00  1.3071e+00  4.1667e+01  1.5633e+00  1
2  3.8817e-01  2.7260e-01  8.6896e+00  4.3973e-01  1
3  1.3674e-01  7.3342e-02  2.3379e+00  2.9856e-01  1
4  4.0684e-02  1.6830e-02  5.3649e-01  1.4986e-01  1
5  1.0436e-02  2.5143e-03  8.0146e-02  6.5562e-02  1
6  1.3460e-03  1.2747e-14  6.1372e-17  8.5276e-03  1
7  1.8750e-04  8.0473e-15  4.5367e-17  1.1522e-03  1
8  2.3185e-05  8.2632e-15  4.5512e-17  1.4949e-04  1
9  2.1053e-06  8.6647e-15  4.3328e-17  1.5936e-05  1
10  9.6541e-08  8.3346e-15  4.2522e-17  1.3858e-06  1

> EXIT -- Below Tolerance!

2.077634 seconds (44.26 k allocations: 1.036 GB, 12.61% gc time)


And as a bonus we even get an extra boost in speed!

In [9]:
HTML(readall(open("style.css")))

Out[9]: