Skip to content

Conic programming: Solving LPs, SOCPs, and SDPs

Overview

Conic programs are a specific subclass of convex optimization problems, where the objective function is linear and the convex constraints are restrictions to the intersection of affine spaces and certain cones within \(\mathbb{R}^n\). Commonly considered cones are the positive orthant, the second-order cone ("ice-cream cone"), and the semidefinite cone, which give rise to linear programs (LPs), second-order cone programs (SOCPs), and semidefinite programs (SDPs), respectively. This framework remains quite general and many real-world problems can be reduced to a conic program. However, the additional structure of the program allows for more efficient classical and quantum algorithms, compared to completely general convex problems.

Algorithms for LPs, SOCPs, and SDPs have long been a topic of study. Today, the best classical algorithms are based on interior point methods (IPMs), but other algorithms based on the multiplicative weights update (MWU) method exist and can be superior in a regime where high precision is not required. Both of these approaches can be turned into quantum algorithms with potential to deliver asymptotic quantum speedup for general LPs, SOCPs, and SDPs. However, the runtime of the quantum algorithm typically depends on additional instance-specific parameters, which makes it difficult to produce a general apples-to-apples comparison with classical algorithms.

Actual end-to-end problem(s) solved

  • Linear programs (LPs) are the simplest convex program. An LP instance is specified by an \(m \times n\) matrix \(A\), an \(n\)-dimensional vector \(c\) and an \(m\)-dimensional vector \(b\). The problem can then be written as
    \[\begin{equation} \label{eq:LP} \begin{align} & \min_{x \in \mathbb{R}^n} \langle c, x \rangle \\ \text{subject to } & Ax=b \\ & x_i \geq 0 \text{ for } i=1,\ldots,n \end{align} \end{equation}\]

    where notation \(\langle u, v \rangle\) denotes the standard dot product of vectors \(u\) and \(v\). The function \(\langle c, x\rangle\), which is linear in \(x\), is called the objective function, and a point \(x\) is called feasible if it satisfies the linear equality1 constraints \(Ax=b\) as well as the positivity constraints \(x_i \geq 0\) for all \(i\). We denote the feasible point that optimizes the objective function by \(x^*\). Let \(\epsilon\) be a precision parameter. The actual end-to-end problem solved is to take as input a classical description of the problem instance \((c,A,b,\epsilon)\) and output a classical description of a point \(x\) for which \(\langle c, x \rangle \leq \langle c, x^*\rangle + \epsilon\). The set of points that obey the positivity constraints \(x_i \geq 0\) forms the positive orthant of the vector space \(\mathbb{R}^n\). This set meets the mathematical definition of a convex cone: for any points \(u\) and \(v\) in the set and any non-negative scalars \(\alpha,\beta\geq 0\), the point \(\alpha u+\beta v\) is also in the set. - Second-order cone programs (SOCPs) are formed by replacing the positivity constraints in the definition of LPs with one or more second-order cone constraints, where the second-order cone of dimension \(k\) is defined to include points \((x_0;x_1;\ldots;x_{k-1}) \in \mathbb{R}^k\) for which \(x_0^2 \geq x_1^2+\ldots + x_{k-1}^2\). - Semidefinite programs (SDPs) are formed by replacing the \(n\)-dimensional vector \(x\) in the definition of LPs with a \(n \times n\) symmetric matrix \(X\) and replacing the positive orthant constraint with the conic constraint that \(X\) is a positive semidefinite matrix. Denote the set of \(n \times n\) symmetric matrices by \(\mathbb{S}^n\), and for any pair of matrices \(U,V \in \mathbb{S}^n\), define the notation \(\langle U, V\rangle = \text{tr}(UV)\) (which generalizes the standard dot product). Then, an SDP instance is specified by matrices \(C, A^{(1)},A^{(2)},\ldots, A^{(m)} \in \mathbb{S}^n\), as well as \(b \in \mathbb{R}^m\), and can be written as

    \[\begin{equation} \label{eq:SDP} \begin{align} & \min_{X \in \mathbb{S}^n} \langle C, X \rangle \\ \text{subject to } & \langle A^{(j)}, X\rangle = b_j \text{ for } j = 1,\ldots,m \\ & X \succeq 0 \end{align} \end{equation}\]

    where \(X\succeq 0\) denotes the constraint that \(X\) is positive semidefinite.

In the LP or SDP case, we might also require as input parameters \(R\) and \(r\), where \(R\) is a known upper bound on the size of the solution in the sense that \(\sum_i |x_i| \leq R\) (LP) or \(\text{tr}(X) \leq R\) (SDP), and where \(r\) is an analogous upper bound on the size of the solution to the dual program (not written explicitly here, see [1]).

Dominant resource cost/complexity

Two separate approaches to solving conic programs with quantum algorithms have been proposed in the literature. Both methods start with classical algorithms and replace some of the subroutines with quantum algorithms.

  1. Quantum interior point methods (QIPMs) for LPs [2], SOCPs [3, 4], and SDPs [2, 5, 6] have been proposed. These methods start with classical interior point methods, for which the core step is solving a linear system, and simply replace the classical linear system solver with a quantum linear system solver (QLSS), combined with pure state quantum tomography. Given a linear system \(Gu = v\), the QLSS produces a quantum state \(\ket{u}\), and quantum tomography is subsequently used to gain a classical estimate of the amplitudes of \(\ket{u}\) in the computational basis. The QLSS ingredient introduces complexity dependence on a parameter \(\kappa=\lVert G \rVert \lVert G^{-1}\rVert\), the condition number of \(G\), where \(\lVert \cdot \rVert\) denotes the spectral norm. Additionally, the QLSS requires that the classical data defining \(G\) be loaded in the form of a block-encoding, for which the standard construction introduces a dependence on the factor \(\zeta = \lVert G \rVert_F \lVert G \rVert^{-1}\), where \(\lVert \cdot \rVert_F\) denotes the Frobenius norm. Finally, the tomography ingredient introduces a complexity dependence on a parameter \(\xi\), defined as the precision to which the vector \(u\) must be classically learned, measured in \(\ell_2\) norm. Assuming \(m\) is on the order of the number of degrees of freedom (i.e. \(n\) in the case of LP and SOCP, and \(n^2\) in the case of SDP), the number of queries the QIPM makes to block-encodings of the input matrices is
    \[\begin{align} \text{LP, SOCP \cite{augustino2022inexact}:} & \qquad \tilde{\mathcal{O}}\left(\frac{n^{1.5}\zeta \kappa}{\xi}\log(1/\epsilon)\right) \\ \text{SDP \cite{kerenidis2018QIntPoint,augustino2021quantum}:} & \qquad \tilde{\mathcal{O}}\left(\frac{n^{2.5}\zeta \kappa}{\xi}\log(1/\epsilon)\right) \end{align}\]

    where the \(\tilde{\mathcal{O}}\) notation hides logarithmic factors. Note that depending on how \(\xi\) is defined, extra factors of \(\kappa\) may be required. Moreover, note that the complexity statements in [5] go further and analyze the worst-case dependence of \(\xi\) on the overall error \(\epsilon\), and additionally make the worst-case replacement \(\zeta \leq \mathcal{O}\left( n \right)\). The numerical values of \(\kappa\), \(\zeta\), and \(\xi\) are generally difficult to determine in advance for a specific application. The block-encoding queries can be executed in circuit depth \(\mathrm{polylog}(n+m,1/\epsilon)\), which can also be absorbed into the \(\tilde{\mathcal{O}}\) notation (although it is important to note that the circuit size is generally \(\mathcal{O}\left( n^2 \right)\)). If the input matrices are sparse or given in a form other than as a list of matrix entries, there may be other more efficient methods for block-encoding; in this case the parameter \(\zeta\) might be replaced with another parameter \(\alpha > 1\), whose value would depend on the block-encoding method. 2. Quantum algorithms based on the multiplicative weights update (MWU) method have been proposed for SDP [7, 8, 9, 1] and LP [9, 10]. The quantum algorithm closely follows the classical algorithm based on MWU to iteratively update a candidate solution to the program. Each iteration is carried out using quantum subroutines, including Gibbs sampling, as well as Grover search and quantum minimum finding [11, 9] (a direct application of Grover search). Let \(s\) denote the sparsity, that is, the maximum number of nonzero entries in any row or column of the matrices composing the problem input (thus, \(s \leq \max(m,n)\)). Then, the runtime has been upper bounded by

    \[\begin{align} \text{LP \cite{bouland2023zerosum}:} & \qquad \tilde{\mathcal{O}}\left(\sqrt{s}\left(\frac{rR}{\epsilon}\right)^{3.5}\right) \\ \text{SDP \cite{apeldoorn2018ImprovedQSDPSolving}:} & \qquad \tilde{\mathcal{O}}\left(s\sqrt{m}\left(\frac{rR}{\epsilon}\right)^{4}+s\sqrt{n}\left(\frac{rR}{\epsilon}\right)^5\right) \end{align}\]

    assuming sparse access to the input matrices, where \(r,R\) are the parameters related to the size of the primal and dual solutions, defined above. In [1], the input model was generalized to a "quantum operator input model," based on block-encodings where \(s\) is replaced by the block-encoding normalization factor \(\alpha\) in the runtime expressions. Note that it is possible the runtime for LP could be improved by applying the dynamic Gibbs sampling method of [12] together with the reduction from LP to zero-sum games in [10].

The runtime expressions for the QIPM approach and the MWU approach are not directly comparable, as the former depends on instance-specific parameters \(\kappa\), \(\zeta\), and \(\xi\), while the latter depends on instance-specific parameters \(r\) and \(R\). However, note that the explicit \(n\)-dependence is better in the case of MWU than QIPM, while the \(\epsilon\)-dependence is worse.

Existing error corrected resource estimates

Neither of the approaches for conic programs have garnered study at the level of error-corrected resource estimates. Reference [13] performed a resource analysis for a QIPM at the logical level, but did not analyze additional overheads due to error correction. The goal of that analysis was to completely compile the QIPM for SOCP into Clifford gates and \(T\) gates, and then to numerically estimate the parameters \(\kappa\), \(\zeta\), and \(\xi\) for the particular use case of financial portfolio optimization, which can be reduced to SOCP. A salient feature of the QIPM is that \(\mathcal{O}\left( n+m \right) \times \mathcal{O}\left( n+m \right)\) matrices of classical data must be repeatedly accessed by the QLSS via block-encoding, necessitating a large-scale quantum random access memory (QRAM) with \(\mathcal{O}\left( n^2 \right)\) qubits. Accordingly, for SOCPs with \(n=500\) and \(m=400\) (which are still easily solved on classical computers) it was estimated that 8 million logical qubits would be needed. The total number of \(T\) gates needed for the same instance size was on the order of \(10^{29}\), which can be distributed over roughly \(10^{24}\) layers.

We are not aware of an analogous logical resource analysis for the MWU approach to conic programming. Such an analysis would be valuable and should ideally choose a specific use case to be able to evaluate the size of all parameters involved.

Caveats

  • The QIPM approach requires a large-scale QRAM of size \(\mathcal{O}\left( n^2 \right)\). This is a necessary ingredient to retain any hope of a speedup, and for relevant choice of \(n\) the associated hardware requirements could be prohibitively large.
  • The QIPM approach has a weak case for a large asymptotic speedup: even under optimal circumstances, the asymptotic speedup over classical interior point methods is less than quadratic. See the article on the QIPM approach for more information.
  • The QIPM approach has a large constant prefactor that is estimated to be on the order of \(10^3\) coming from state-of-the-art QLSS [14, 15]. (It is possible this could be improved using alternative approaches to QLSS such as variable-time amplitude amplification [16]. See also [15].)
  • The MWU approach requires a medium-scale QRAM of size \(\mathcal{O}\left( R^2r^2/\epsilon^2 \right)\). This requirement can be avoided at the cost of an additional multiplicative overhead of \(\mathcal{O}\left( R^2r^2/\epsilon^2 \right)\).
  • The MWU approach has poor dependence on error \(\epsilon\); for SDPs it is \(\epsilon^{-5}\). Even at modest choices of \(\epsilon\), this may lead the algorithm to be impractical pending significant improvements.
  • A general caveat that applies to both approaches is that the appearance of instance-specific parameters makes it difficult to predict the performance of these algorithms for more specific applications.

Comparable classical complexity and challenging instance sizes

As in the quantum case, there are multiple distinct approaches in the classical case.

  1. Classical interior point methods (CIPMs): There exist fast IPM-based software implementations for solving conic programs, such as ECOS [17] and MOSEK [18]. These solvers can solve instances with thousands of variables in a matter of seconds on a standard laptop [17]. However, the runtime scaling is poor and scaling too far beyond this regime leads the solvers to be far less practical. The runtime of the best provably correct classical IPMs for the regime where the number of constraints is roughly equal to the number of degrees of freedom is
    \[\begin{align} \text{LP \cite{cohen2021LPsinMMtime}:} & \qquad \tilde{\mathcal{O}}\left(n^{\omega}\log(1/\epsilon)\right) \\ \text{SOCP \cite{monteiro2000SOCP}:} & \qquad \tilde{\mathcal{O}}\left(n^{\omega+0.5}\log(1/\epsilon)\right) \\ \text{SDP \cite{huang2022fasterIPM}:} & \qquad \tilde{\mathcal{O}}\left( n^{2\omega}\log(1/\epsilon)\right) \end{align}\]

    where \(\omega<2.37\) is the matrix multiplication exponent. It is plausible that, with some attention, the extra \(n^{0.5}\) factor for SOCP could be eliminated with modern techniques. Additionally, the runtime can be somewhat reduced when the number of constraints is much less than the number of degrees of freedom; for example, the \(n\)-dependence of the complexity of the CIPM for SDP in [19] can be as low as \(\tilde{\mathcal{O}}(n^{2.5})\) when there are few constraints. On practical instances, employing techniques for fast matrix multiplication is often not beneficial, and Gaussian Elimination-like methods are used, where \(\omega=3\). Note that, alternatively, by using iterative classical linear systems solvers [20], each \(n^{\omega}\) factor could be replaced by a factor of \(n\) at the cost of a linear dependence on \((\kappa \zeta)^2\), which could be superior if the matrices are well conditioned. 2. Classical MWU methods: a classical complexity statement for LPs is inferred from the reduction in [10] from LPs to zero-sum games and the classical analysis that appears there. For the SDP case, references in the classical literature appear only to examine specific subclasses of SDPs (e.g. [21, 22]). A general statement of the classical complexity for SDPs appears alongside the quantum algorithm in [9, Section 2.4].

    \[\begin{align} \text{LP \cite{apeldoorn2019QAlgorithmsForZeroSumGames}:} & \qquad \tilde{\mathcal{O}}\left(s\left(\frac{rR}{\epsilon}\right)^{3.5}\right) \\ \text{SDP \cite{apeldoorn2017QSDPSolvers}:} & \qquad \tilde{\mathcal{O}}\left(s\sqrt{nm}\left(\frac{rR}{\epsilon}\right)^{4}+sn\left(\frac{rR}{\epsilon}\right)^7\right) \end{align}\]
  2. Cutting-plane methods: these classical methods are used for SDPs and can outperform IPMs when the number of constraints is small. The best algorithm, based on [23, 24], has runtime \(\mathcal{O}\left( m(mn^2+n^{\omega}+m^2)\log(1/\epsilon) \right)\), which can be as low as \(\mathcal{O}\left( n^\omega \right)\) when \(m\) is small.

It is important to note that the algorithms with the best provable complexities may not be the ones that are most useful in practice.

Speedup

For both the IPM and the MWU approach, there can be at most a polynomial quantum speedup: upper and lower bounds scaling polynomially with \(n\) are known in both the classical and quantum cases [1]. The speedup of the QIPM method depends on the scaling of \(\kappa\) with \(n\), but the speedup cannot be more than quadratic. The speedup of the MWU method with respect to the \(n\)-scaling could be as much as quadratic, assuming the sparsity is constant. There is a possibility that the speedup could be larger in practice if the Gibbs sampling routine is faster in practice than its worst-case upper bounds suggest, perhaps by utilizing Monte Carlo–style approaches to Gibbs sampling.

Outlook

It is very plausible that an asymptotic polynomial speedup can be obtained in problem size using the MWU method for solving LPs or SDPs, but the speedup appears only quadratic, and an assessment of practicality depends on the scaling of certain unspecified instance-specific parameters. Similarly, the QIPM method could bring a subquadratic speedup but only under certain assumptions about the condition number of certain matrices. These quadratic and subquadratic speedups alone might be regarded as unlikely to yield practical speedups after error correction overheads and slower quantum clock speeds are considered. Future work should aim to find additional asymptotic speedups while focusing on specific practically relevant use cases that allow the unspecified parameters to be evaluated.

Bibliography

  1. Joran van Apeldoorn and András Gilyén. Improvements in quantum sdp-solving with applications. In Proceedings of the 46th International Colloquium on Automata, Languages, and Programming (ICALP), 99:1–99:15. 2019. arXiv: https://arxiv.org/abs/1804.05058. doi:10.4230/LIPIcs.ICALP.2019.99.

  2. Iordanis Kerenidis and Anupam Prakash. A quantum interior point method for lps and sdps. ACM Transactions on Quantum Computing, 2020. arXiv: https://arxiv.org/abs/1808.09266. doi:10.1145/3406306.

  3. Iordanis Kerenidis, Anupam Prakash, and Dániel Szilágyi. Quantum algorithms for second-order cone programming and support vector machines. Quantum, 5:427, 2021. arXiv: https://arxiv.org/abs/1908.06720. doi:10.22331/q-2021-04-08-427.

  4. Brandon Augustino, Tamás Terlaky, and Luis F Zuluaga. An inexact-feasible quantum interior point method for second-order cone optimization. Technical Report 21T-009, Department of Industrial and Systems Engineering, Lehigh University, 2022.

  5. Brandon Augustino, Giacomo Nannicini, Tamás Terlaky, and Luis F. Zuluaga. Quantum interior point methods for semidefinite optimization. Quantum, 7:1110, 9 2023. arXiv: https://arxiv.org/abs/2112.06025. URL: https://doi.org/10.22331/q-2023-09-11-1110, doi:10.22331/q-2023-09-11-1110.

  6. Baihe Huang, Shunhua Jiang, Zhao Song, Runzhou Tao, and Ruizhe Zhang. A faster quantum algorithm for semidefinite programming via robust ipm framework. arXiv: https://arxiv.org/abs/2207.11154, 2022.

  7. Fernando G. S. L. Brandão and Krysta M. Svore. Quantum speed-ups for solving semidefinite programs. In Proceedings of the 58th IEEE Symposium on Foundations of Computer Science (FOCS), 415–426. 2017. arXiv: https://arxiv.org/abs/1609.05537. URL: http://ieee-focs.org/FOCS-2017-Papers/3464a415.pdf, doi:10.1109/FOCS.2017.45.

  8. Fernando G. S. L. Brandão, Amir Kalev, Tongyang Li, Cedric Yen-Yu Lin, Krysta M. Svore, and Xiaodi Wu. Quantum sdp solvers: large speed-ups, optimality, and applications to quantum learning. In Proceedings of the 46th International Colloquium on Automata, Languages, and Programming (ICALP), 27:1–27:14. 2019. arXiv: https://arxiv.org/abs/1710.02581. doi:10.4230/LIPIcs.ICALP.2019.27.

  9. Joran van Apeldoorn, András Gilyén, Sander Gribling, and Ronald de Wolf. Quantum sdp-solvers: better upper and lower bounds. Quantum, 4:230, 2020. Earlier version in FOCS'17. arXiv: https://arxiv.org/abs/1705.01843. doi:10.22331/q-2020-02-14-230.

  10. Joran van Apeldoorn and András Gilyén. Quantum algorithms for zero-sum games. arXiv: https://arxiv.org/abs/1904.03180, 2019.

  11. Christoph Dürr and Peter Høyer. A quantum algorithm for finding the minimum. arXiv: https://arxiv.org/abs/quant-ph/9607014, 1996.

  12. Adam Bouland, Yosheb Getachew, Yujia Jin, Aaron Sidford, and Kevin Tian. Quantum speedups for zero-sum games via improved dynamic gibbs sampling. arXiv: https://arxiv.org/abs/2301.03763, 2023.

  13. Alexander M Dalzell, B David Clader, Grant Salton, Mario Berta, Cedric Yen-Yu Lin, David A Bader, Nikitas Stamatopoulos, Martin J A Schuetz, Fernando G S L Brandão, Helmut G Katzgraber, and others. End-to-end resource analysis for quantum interior point methods and portfolio optimization. PRX Quantum, pages to appear, 2023. arXiv: https://arxiv.org/abs/2211.12489.

  14. Pedro C.S. Costa, Dong An, Yuval R. Sanders, Yuan Su, Ryan Babbush, and Dominic W. Berry. Optimal scaling quantum linear-systems solver via discrete adiabatic theorem. PRX Quantum, 3:040303, 10 2022. arXiv: https://arxiv.org/abs/2111.08152. URL: https://link.aps.org/doi/10.1103/PRXQuantum.3.040303, doi:10.1103/PRXQuantum.3.040303.

  15. David Jennings, Matteo Lostaglio, Sam Pallister, Andrew T. Sornborger, and Yigit Subasi. Efficient quantum linear solver algorithm with detailed running costs. arXiv: https://arxiv.org/abs/2305.11352, 2023.

  16. Andris Ambainis. Variable time amplitude amplification and quantum algorithms for linear algebra problems. In Proceedings of the 29th Symposium on Theoretical Aspects of Computer Science (STACS), 636–647. 2012. arXiv: https://arxiv.org/abs/1010.4458. doi:10.4230/LIPIcs.STACS.2012.636.

  17. Alexander Domahidi, Eric Chu, and Stephen Boyd. Ecos: an socp solver for embedded systems. In 2013 European Control Conference (ECC), volume, 3071–3076. 2013. doi:10.23919/ECC.2013.6669541.

  18. Erling D. Andersen and Knud D. Andersen. The Mosek Interior Point Optimizer for Linear Programming: An Implementation of the Homogeneous Algorithm, pages 197–232. Springer US, Boston, MA, 2000. URL: https://doi.org/10.1007/978-1-4757-3216-0\_8, doi:10.1007/978-1-4757-3216-0\_8.

  19. Haotian Jiang, Tarun Kathuria, Yin Tat Lee, Swati Padmanabhan, and Zhao Song. A faster interior point method for semidefinite programming. In Proceedings of the 61st IEEE Symposium on Foundations of Computer Science (FOCS), volume, 910–918. 2020. arXiv: https://arxiv.org/abs/2009.10217. doi:10.1109/FOCS46700.2020.00089.

  20. Thomas Strohmer and Roman Vershynin. A randomized kaczmarz algorithm with exponential convergence. Journal of Fourier Analysis and Applications, 15(2):262–278, 2009. arXiv: https://arxiv.org/abs/math/0702226. doi:10.1007/s00041-008-9030-4.

  21. S. Arora, E. Hazan, and S. Kale. Fast algorithms for approximate semidefinite programming using the multiplicative weights update method. In Proceedings of the 46th IEEE Symposium on Foundations of Computer Science (FOCS), volume, 339–348. 2005. doi:10.1109/SFCS.2005.35.

  22. Sanjeev Arora and Satyen Kale. A combinatorial, primal-dual approach to semidefinite programs. In Proceedings of the 39th ACM Symposium on the Theory of Computing (STOC), 227–236. New York, NY, USA, 2007. Association for Computing Machinery. URL: https://doi.org/10.1145/1250790.1250823, doi:10.1145/1250790.1250823.

  23. Yin Tat Lee, Aaron Sidford, and Sam Chiu-wai Wong. A faster cutting plane method and its implications for combinatorial and convex optimization. In Proceedings of the 56th IEEE Symposium on Foundations of Computer Science (FOCS), 1049–1065. 2015. arXiv: https://arxiv.org/abs/1508.04874. doi:10.1109/FOCS.2015.68.

  24. Haotian Jiang, Yin Tat Lee, Zhao Song, and Sam Chiu-Wai Wong. An improved cutting plane method for convex optimization, convex-concave games, and its applications. In Proceedings of the 52nd ACM Symposium on the Theory of Computing (STOC), 944–953. New York, NY, USA, 2020. Association for Computing Machinery. arXiv: https://arxiv.org/abs/2004.04250. URL: https://doi.org/10.1145/3357713.3384284, doi:10.1145/3357713.3384284.

  1. Inequality constraints of the form \(Ax \leq b\) can be converted to linear equality constraints and positivity constraints by introducing a vector of slack variables \(s\) and imposing \(Ax + s = b\) and \(s_i \geq 0\) for all \(i\). An analogous trick is possible for SOCP and SDP. ↩