A Quantum Interior Point Method for
LPs and SDPs
Iordanis Kerenidis 1

Anupam Prakash 1

1 CNRS, IRIF, Université Paris Diderot, Paris, France.

September 26, 2018

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Semi Definite Programs

A Semidefinite Program (SDP) is an optimization problem
with inputs a vector c ∈ Rm and matrices A(1) , . . . , A(m) , B in
Rn×n .
The primal and dual SDP are given by,X
Opt(P) = minm {c t x |xk A(k)  B}.
x∈Rk∈[m]Opt(D) = max{Tr (BY ) | Y  0, Tr (YA(j) ) = cj }.
Y 0

We will be working in the case where P, D are strictly feasible,
in this case strong duality holds and Opt(P) = Opt(D).

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Semi Definite Programs

A linear programs (LP) is a special cases of an SDP where(A(i) , B) are diagonal matrices,Opt(P) = minm {c t x |x∈RXxi ai ≥ b, ai ∈ Rn }i∈[m]Opt(D) = max{b t y | y t aj = cj }y ≥0SDPs are one of the most general class of optimization
problems for which we have efficient algorithms.
SDPs capture a large class of convex optimization problems,
they are also used for approximate solutions to NP hardproblems.

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

SDP algorithms

The running time for SDP algorithms will be given in terms ofm, n and , here we consider m = O(n2 ).
The first polynomial time algorithms for LPs and SDPs wereobtained using the ellipsoid method and the interior pointmethod.
The best known LP and SDP algorithms have complexity
O(n2.5 log(1/)) and O((m3 + mnω + mn2 s) log(mn/)).
In addition there is the Arora Kale method, whose complexity
is upper bounded by,Õ nmsIordanis Kerenidis , Anupam Prakash

Rr

4


+ nsRr

7



QuICS Workshop, Maryland, 2018

Quantum SDP algorithms

Quantum SDP algorithms based on Arora-Kale framework
were proposed by [Brandao-Svore 17] and subsequently
improved by [van Appeldoorn-Gribling-Gilyen-de Wolf 17].
These algorithms were recently improved even further in[BKLLSW18] and [AG18].
The best known running time for a quantum SDP algorithm
using the Arora-Kale framework is,!
   4√√√
RrRrÕm+ nn .
For Max-Cut and scheduling LPs , the complexity is at leastO(n6 ) [AGGW17, Theorem 20].
Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Our Results

We provide a quantum interior point method with complexity
e n2.5e n1.5O(µκ3 log(1/)) for SDPs and O(µκ3 log(1/)) forξ2ξ2LPs .
The output of our algorithm is a pair of matrices (S, Y ) thatare -optimal ξ-approximate SDP solutions.
√√
The parameter µ is at most 2n for SDPs and 2n for LPs .
The parameter κ is an upper bound on the condition numberof the intermediate solution matrices.
If the intermediate matrices are ’well conditioned’, the runninge 3.5 ) and O(ne 2 ).
time scales as O(nIordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Input models
Sparse oracle model [BS16, AGGW17]: The input matrices
A(i) are assumed to be s-sparse and we have access toOA : |i, k, l, 0i → |i, k, l, index(i, k, l)i.
(i)

(i)

Quantum state model [BKLLSW17]: A(i) = A+ + A− andwe have access to the purifications of the corresponding
density matrices for all i ∈ [m].
Operator model [AG18]: Access to unitary block encodings ofthe A(i) , that is implementations of: (j)A /αj .
Uj =.
.
.
How to construct block encodings for A and what α can beachieved?
Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

QRAM data structure model
QRAM data structure model: Access to efficient data
structure storing A(i) , i ∈ [m] in a QRAM (Quantum Random
Access Memory).
Given di , i ∈ [N] stored in the QRAM, the following queriesrequire time polylog(N),|i, 0i → |i, xi i
Definition
A QRAM data structure for storing a dataset D of size N isefficient if it can be constructed in a single pass over the entries(i, di ) for i ∈ [N] and the update time per entry is O(polylog(N)).
Introduced to address the state preparation problem in
Quantum Machine Learning, to prepare arbitrary vector states√|xi without incurring an O( n) overhead.
Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Block encodings using QRAM
The optimal value of α = kAk, any minor of a unitary matrix
has spectral norm at most 1.
The quantum linear system problem with A ∈ Rn×n is to
produce states |Abi , |A−1 bi, it is scale invariant so we assume
kAk = 1.
PDefine sp (A) = maxi∈[n] j∈[n] Apij and letqµ(A) = minp∈[0,1] (kAkF , s2p (A)s(1−2p) (AT )).
Theorem (KP16, KP17, CGJ18)
There are efficient QRAM data structures, that allow a block
encodings for A ∈ Rn×n with α = µ(A) to be implemented in timeO(polylog(n)).
√√
Notice that µ(A) < n is sub-linear, it can be O( n) in theworst case.
Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Quantum linear system solvers
Given an efficient block encoding for A, there is a quantum
2 (A)/)elinear system solver with running time O(µ(A)κ
[KP16, KP17].
Given an efficient block encoding for A, there is a quantumelinear system solver with running time O(µ(A)κ(A)log(1/)).
[CGJ18, GSLW18].
Composing block encodings: Given block encodings for
M1 , M2 with paramaters µ1 , µ2 , the linear system in
M = M1 M2 can be solved in time O((µ1 + µ2 )κ(M) log(1/)).
Can quantum linear systems be leveraged for optimization
using iterative methods? Gradient descent with affine update
rules. [KP17].
This work: Quantum LP and SDP solvers are not muchharder than quantum linear systems!
Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Interior Point Method overview

The classical IPM starts with feasible solutions (S, Y ) to theSDP and updates them (S, Y ) → (S + dS, Y + dY )iteratively.
The updates (dS, dY ) are obtained by solving a n2 × n2 linearsystem called the Newton linear system.
The matrix for the Newton linear system is not explicit and isexpensive to compute from the data.
√After O( n log(1/)) iterations, the method converges tofeasible solutions (S, Y ) with duality gap at most Iordanis Kerenidis , Anupam PrakashQuICS Workshop, Maryland, 2018

Classical Interior Point MethodAlgorithm 1 Classical interior point method.
Require: Matrices A(k) with k ∈ [m], B ∈ Rn×n , c ∈ Rm in memory, precision  > 0.
1

23Find feasible initial point (S0 , Y0 , ν0 ) close to the analyticcenter.
Startingwith (S0 , Y0 , ν0 ) repeat the following steps√
O( n log(1/)) times.
1

Solve the Newton linear system to get (dS, dY ).

2Update S ← S + dS, Y ← Y + dY , ν ← Tr (SY )/n.

Output (S, Y ).

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Quantum Interior Point Method overview

We construct block encodings for the Newton linear system
matrix which allows us to solve this linear system with low
cost in the quantum setting.
We need to find (dS, dY ) classically to write the Newton
linear system for the next iteration.
We give a linear time tomography algorithm that reconstructs
dd-dimensional state to error δ with complexity O( d log).
δ2We show that with tomography precision δ = O( κ1 ) themethod converges at the same rate as the classical IPM.
The solutions output by the QIPM are ξ-approximately
feasible as (dS, dY ) are reconstructed by tomography.

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Quantum Interior Point Method
Algorithm 2 Quantum interior point method.
Require: Same as classical algorithm with inputs stored in QRAM.
1

2

Find feasible initial point (S, Y , ν) = (S0 , Y0 , ν0 ) andstore in QRAM.
Repeat the for T iterations and output the final (S, Y ).
1

Solve Newton linear system to obtain state close to|dS ◦ dY i to error δ 2 /n3 .

2Estimate kdS ◦ dY k, perform tomography and use thenorm estimate to obtain,dS ◦ dY − dS ◦ dY32

≤ 2δ kdS ◦ dY k2 .

Update Y ← Y + dY and S ← S + dS and store in
QRAM. Update ν ← Tr (SY )/n.

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Running time overview

2.5e n 2 µκ3 log(1/)).
Running time for SDPs is O(ξThere are O(n0.5 log(1/)) iterations.
In each iteration, we solve Newton linear system having sizeO(n2 ) in time O(µκ log(1/)).
2 2

We then perform tomography in time O( nξκ2 ).
e n1.5
Running time for LPs is O(µκ3 log(1/)), the linearξ2system has size O(n).

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Talk Overview

Tomography for efficient vector states.
Classical interior point method.
Analysis of the approximate interior point method.
Quantum interior point method.

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Vector state tomography

The vector state for x ∈ Rd is defined as |xi =P

i xi |ii.

We assume that |xi can be prepared efficiently, that is wehave have access to a unitary U that prepares copies of |xi.
We need to learn sgn(xi ) for vector state tomography, wewould need to learn a phase e −2πθi for pure state tomography.
Theoremd)
There is an algorithm with time and query complexity O( d logδ2dthat produces an estimate xe ∈ R with ke
x k2 = 1 such thatkex − xk2 ≤ δ with probability at least (1 − 1/poly (d)).

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Vector state tomography algorithm
Measure N = 36dδ2ln d copies of |xi in the standard basis toobtain estimates pi = nNi where ni is the number of times
outcome i is observed.
√Store Ppi for i ∈ [d] in QRAM data structure so that√|pi = i∈[d] pi |ii can be prepared efficiently.
PP√
Create N copies of √12 |0i i∈[d] xi |ii + √12 |1i i∈[d] pi |ii
using a control qubit.
Apply a Hadamard
gate on the first qubit of each copy of the√√1P
state to obtain 2 i∈[d] [(xi + pi ) |0, ii + (xi − pi ) |1, ii].
Measure each copy in the standard basis and maintain countsn(b, i) of the number of times outcome |b, ii is observed forb ∈ 0, 1. Set σi = 1 if n(0, i) > 0.4pi N and −1 otherwise.
√Output the unit vector xe with xei = σi pi .
Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Analysis of tomography algorithm
Let S = {i ∈ [d] | xi2 ≥ δ 2 /d}, note that S has at most a δ 2
fraction of the `2 norm.
Claim1: The sign is estimated correctly for all i ∈ S withprobability at least 1 − 1/poly (d).
√√Claim2: For all i ∈ S we have that |xi − pi | ≤ δ/ d withprobability at least 1 − 1/poly (d).
PClaim3: i∈S pi ≤ 2δ 2 with probability at least
1 − 1/poly (d).
Error analysis using the three claims:
XXX
√
√
√
(xi − σ(i) pi )2 =(|xi | − pi )2 +(|xi | + pi )2i∈[d]i∈S2i∈S≤δ +2Xi∈SIordanis Kerenidis , Anupam Prakash(xi2 + pi ) ≤ 3δ 2 + 2Xi∈S

QuICS Workshop, Maryland, 2018pi

Classical Interior Point Method
Recall the primal and dual SDP ,Opt(P) = minm {c t x |x∈RXxk A(k)  B}.

k∈[m]

Opt(D) = max{Tr (BY ) | Y  0, Tr (YA(j) ) = cj }.
Y 0Define L = Spank∈[m] (A(k) ), let L⊥ be the orthogonal
complement.
Let CPbe an arbitrary dual feasible solution andS = k∈[m] xk A(k) − B.
The SDP pair can be written in a more symmetric form,Opt(P 0 ) = min{Tr (CS) + Tr (BC ) | S ∈ (L − B)}
S0Opt(D) = max{Tr (BY ) | Y ∈ (L⊥ + C )}Y 0Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Classical Interior Point Method
As Tr ((S + B)(Y − C )) = 0 we have that the duality gap isTr (SY ).
The logarithmic barrier is defined as K (X ) = − log(det(X )),it is defined on the interior of the psd cone.
The central path is parametrized by ν and is given by thesolutions to,Opt(Pν ) = min{Tr (CS) + νK (S) | S ∈ (L − B)}
S0Opt(Dν ) = max{Tr (BY ) − νK (Y ) | Y ∈ (L⊥ + C )}Y 0As ν → 0 we recover solutions to the original SDP .
Theorem: The optimal solutions (Sν , Yν ) on the central pathsatisfy Sν Yν = Yν Sν = νI .
Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Classical Interior Point Method

An ideal interior point method would follow the central pathin the direction ν → 0. The actual method stays close to thepath.
2Define the distance d(S, Y , ν) = I − ν −1 S 1/2 YS 1/2 F .
Theorem: The duality gap and distance from central path are
related as,√√
ν(n − nd(S, Y , ν)) ≤ Tr (SY ) ≤ ν(n + nd(S, Y , ν))It suffices to stay close to the central path, if d(S, Y , ν)) ≤ ηfor η ∈ [0, 1] then Tr (SY ) ≤ 2νn.

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Classical Interior Point Method
The interior point method starts with a pair of feasible
solutions (S, Y ) with duality gap Tr (SY ) = νn andd(S, Y , ν) ≤ η for a constant η ≤ 0.1.
A single step of the method updates the solution to(S 0 = S + dS, Y 0 = Y + dY ) such that Tr (S 0 Y 0 ) = ν 0 n for√ν 0 = (1 − χ/ n)ν where χ ≤ η is a positive constant.
The updates are found by solving the Newton linear system,dS ∈ L, dY ∈ L⊥
dSY + SdY = ν 0 I − SY
The classical analysis shows that: (i) The Newton linear
system has a unique solution. (ii) The updated solutions(S 0 , Y 0 ) are positive definite. (iii) The distanced(S, Y , ν 0 ) ≤ η.
Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

The approximate interior point method

Approximate interior point method analysis:Theorem
ξand dS − dS F ≤ kYξk beLet dY − dY F ≤ kY −1k2
2
approximate solutions to the Newton linear system and let(S = S + dS, Y = Y + dY ) be the updated solution. Then,• The updated solution is positive definite, that is S  0 andY  0.
• The updated solution satisfies d(S, Y , ν) ≤ η and Tr (S Y ) = νnfor ν = (1 − √αn )ν for a constant 0 < α ≤ χ.

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

The quantum interior point method
How to solve the Newton linear system? If we use thevariables (dS, dY ) then dS ∈ L is hard to express.
With variables (dx, dY ) the Newton linear system is given byM(dx, dY ) = (ν 0 I − SY , 0m ).
(A(1) Y )11..
 (1).
 (A Y )ij


..
.
M=
(A(1) Y )nn0

..
.
0

. . . (A(m) Y )11..
..
.
.
(m). . . (A Y )ij
..
..
.
.
. . . (A(m) Y )nn...
0..
..
.
.
...
0Iordanis Kerenidis , Anupam Prakash
(1 ⊗ S1 )T..
.
T
(j ⊗ Si )


..
.
T(n ⊗ Sn ) (vec(A(1) ))T ..
.
(m)T(vec(A ))QuICS Workshop, Maryland, 2018

The quantum interior point method
Let Z ∈ Rn×n , then define matrix Ze with rows Zeij = (i ⊗ Zj )Tand Zb with rows Zbij = (j ⊗ Zi )T .
Ze is a block diagonal matrix with n copies of Z on diagonal
blocks while Zb is obtained by permuting the rows of Ze.
The Newton linear system matrix can be factorized as aproduct of matrices,
Ye 0AT
M = M1 M2 =.
0 Im0!
g−1bY SA

Block encoding for M1 : Is the same as constructing a blockencoding for Y .
Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

The quantum interior point method−1 Sb are tensor products ofBlock encoding for M2 : Rows of Yg−1gthe rows of Yand S, that is Y −1 Sb = (Yj−1 ⊗ Si )T .

It suffices to prepare the rows and columns of M2 efficiently, if
we pre-compute Y −1 the rows and columns can be preparedefficiently.
In addition we provide a procedure for preparing |a ◦ bi given
unitaries for preparing |ai , |bi in time O(T (Ua ) + T (Ub )).
Further technical details: Recovery of dS, the precision O( κ1 )for tomography, linear programs.
Guarantees: Tr (SY ) ≤  and (S, Y ) ∈ (L − B 0 , L⊥ + C 0 ) withkB ⊕ C − B 0 ⊕ C 0 k ≤ ξ kB ⊕ C kF .

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018

Open questions

For what optimization problems can one get a polynomial
quantum speedup?
Find quantum analogs of interior point methods for the caseof sparse SDPs .
Improve the classical step in the IPM, find better quantum
algorithms for LPs .
Quantum algorithms for convex optimization with provable
polynomial speedups?

Iordanis Kerenidis , Anupam Prakash

QuICS Workshop, Maryland, 2018