Solving Conic Programs over Sparse Graphs using a Variational Quantum Approach: The Case of the Optimal Power Flow

TL;DR

This paper introduces a variational quantum approach to solve large-scale conic programs over sparse graphs, demonstrated on the optimal power flow problem, offering a potentially more scalable alternative to classical methods.

Contribution

It develops a hybrid quantum-classical primal-dual framework using parameterized quantum circuits for solving conic programs, including OPF, with improved efficiency for large sparse graphs.

Findings

Successfully applied to IEEE 57-node power system

High-quality solutions obtained with quantum approach

Framework applicable to various conic problems and machine learning constraints

Abstract

Conic programs arise broadly in physics, quantum information, machine learning, and engineering, many of which are defined over sparse graphs. Although such problems can be solved in polynomial time using classical interior-point solvers, the computational complexity scales unfavorably with graph size. In this context, this work proposes a variational quantum paradigm for solving conic programs, including quadratically constrained quadratic programs (QCQPs) and semidefinite programs (SDPs). We encode primal variables via the state of a parameterized quantum circuit (PQC), and dual variables via the probability mass function of a second PQC. The Lagrangian function can thus be expressed as scaled expectations of quantum observables. A primal-dual solution can be found by minimizing/maximizing the Lagrangian over the parameters of the first/second PQC. We pursue saddle points of the Lagrangian in a hybrid fashion. Gradients of the Lagrangian are estimated using the two PQCs, while PQC parameters are updated classically using a primal-dual method. We propose permuting the primal variables so that related observables are expressed in a banded form, enabling efficient measurement. The proposed framework is applied to the OPF problem, a large-scale optimization problem central to the operation of electric power systems. Numerical tests on the IEEE 57-node power system using Pennylane's simulator corroborate that the proposed doubly variational quantum framework can find high-quality OPF solutions. Although showcased for the OPF, this framework features a broader scope, including conic programs with numerous variables and constraints, problems defined over sparse graphs, and training quantum machine learning models to satisfy constraints.

Figures

Tables

Table 1. Table 1: Average relative errors and constraint violation statistics across problem instances: (a) Number of constraint violations per instance; (b) maximum constraint violation [%]; (c) average constraint violation [%].

Model	${𝐱}_g$ [%]	$𝝀$ [%]	(a)	(b)	(c)
QCQP-PD	14.95	19.30	25.87	13.73	0.49
QCQP-EG	13.67	19.20	24.53	12.29	0.39
${\text{QCQP}}_{\theta }\text{-PD}$	11.82	12.87	11.73	13.03	0.23
${\text{QCQP}}_{\theta }\text{-EG}$	7.62	12.17	11.53	11.86	0.21

Table 2. Table 2: Tested PQC architectures.

No.	architecture of 1 layer	$L$ (primal)	$P$	error (primal)	$L$ (dual)	$Q$	error (dual)
1	${𝐑}_x-$CX	20	120	$9\times {10}^{-1}$	35	315	$8\times {10}^{-5}$
2	${𝐑}_y-$CX	20	120	$5\times {10}^{-3}$	35	315	$2\times {10}^{-5}$
3	${𝐑}_z-$CX	20	120	1	35	315	$2\times {10}^{-3}$
4	${𝐑}_x-$CX$-{𝐑}_y-$CX	10	120	$3\times {10}^{-4}$	18	324	$7\times {10}^{-5}$
5	${𝐑}_x-$CX$-{𝐑}_z-$CX	10	120	$3\times {10}^{-4}$	18	324	$7\times {10}^{-5}$
6	${𝐑}_y-$CX$-{𝐑}_z-$CX	10	120	$2\times {10}^{-4}$	18	324	$7\times {10}^{-5}$
7	${𝐑}_x-$CX$-{𝐑}_y-$CX$-{𝐑}_z-$CX	6	108	$3\times {10}^{-3}$	12	324	$7\times {10}^{-5}$
8	${𝐑}_x-{𝐑}_y-{𝐑}_z-$CX	7	126	${10}^{-3}$	12	324	$4\times {10}^{-6}$

Equations

$$\displaystyle\vartheta(\mathcal{G}):=\max_{\mathbf{X}\succeq 0}\penalty 10000\$$

$$\text{penalty}10000\mathrm{Tr}(\mathbf{X})=1,$$

$$\text{penalty}10000{\mathbf{X}}_{i,j}=0\forall (i,j)\in \mathcal{E},$$

$$\displaystyle P^{*}\coloneqq\min_{\mathbf{x}\in\mathbb{C}^{N}}\penalty 10000\$$

$$\text{penalty}10000{\mathbf{x}}^{†}{\mathbf{M}}_m\mathbf{x}\le b_m,\text{for}\text{penalty}10000\text{penalty}10000m=1:M$$

$$\mathcal{L}(\mathbf{x};\mathbf{\lambda })={\mathbf{x}}^{†}{\mathbf{M}}_0\mathbf{x}+\sum _{m=1}^M{\lambda }_m({\mathbf{x}}^{†}{\mathbf{M}}_m\mathbf{x}-b_m).$$

$$D^{\ast }:=\underset{\mathbf{\lambda }⪰0}{\max }\underset{\mathbf{x}\in {\mathbb{C}}^N}{\min }\text{penalty}10000\mathcal{L}(\mathbf{x};\mathbf{\lambda }).$$

$$\mathcal{L}(\mathbf{x}^{*};\boldsymbol{\lambda})\leq\mathcal{L}(\mathbf{x}^{*};\boldsymbol{\lambda}^{*})\leq\mathcal{L}(\mathbf{x};\boldsymbol{\lambda}^{*})\quad\penalty 10000\$$

$${\mathbf{x}}^{t+1}$$

$${\mathbf{\lambda }}^{t+1}$$

$$\frac{\partial \mathcal{L}}{\partial {\lambda }_m}={\mathbf{x}}^{†}{\mathbf{M}}_m\mathbf{x}-b_m,$$

$${\nabla }_{\mathbf{x}}\mathcal{L}=2\left({\mathbf{M}}_0+\sum _{m=1}^M{\lambda }_m{\mathbf{M}}_m\right)\mathbf{x}$$

$$\displaystyle\min_{\mathbf{x}\in\mathbb{C}^{N}}\penalty 10000\$$

$$\text{penalty}10000{\mathbf{x}}^{†}\mathbf{x}=1.$$

$$\mathbf{x}(\mathbf{\theta })=\alpha |\mathbf{\psi }(\mathbf{\theta })⟩$$

$$\displaystyle P_{\theta}^{*}\coloneqq\min_{\boldsymbol{\theta},\alpha>0}\penalty 10000\$$

$$\displaystyle\mathrm{s.to}\penalty 10000\$$

$$F_m(\mathbf{\theta }):=⟨\mathbf{\psi }(\mathbf{\theta })|{\mathbf{M}}_m|\mathbf{\psi }(\mathbf{\theta })⟩,m=0:M.$$

$$F_m(\mathbf{\theta })=\mathrm{Tr}({\mathbf{M}}_m\mathbf{\rho }(\mathbf{\theta })),m=0:M.$$

$$D_{\theta }^{\ast }:=\underset{\mathbf{\lambda }\ge 0}{\max }\text{penalty}10000\underset{\mathbf{\theta },\alpha >0}{\min }\mathcal{L}(\mathbf{\theta },\alpha ;\mathbf{\lambda })$$

$$\mathcal{L}(\mathbf{\theta },\alpha ;\mathbf{\lambda })$$

$${\lambda }_m(\mathbf{\varphi })={\beta }^2\cdot |{\xi }_m(\mathbf{\varphi }){|}^2,\text{for}\text{penalty}10000m=1:M$$

$$\mathcal{L}(\mathbf{\theta },\alpha ;\mathbf{\varphi },\beta )={\alpha }^2{F}_0(\mathbf{\theta })+{\alpha }^2{\beta }^{2}F(\mathbf{\theta },\mathbf{\varphi })-{\beta }^{2}G(\mathbf{\varphi }),$$

$$F(\mathbf{\theta },\mathbf{\varphi }):=\sum _{m=1}^M|{\xi }_m(\mathbf{\varphi }){|}^2{F}_m(\mathbf{\theta })$$

$$G(\mathbf{\varphi }):=\sum _{m=1}^M{b}_m|{\xi }_m(\mathbf{\varphi }){|}^2.$$

$$D_{\theta ,\varphi }^{\ast }:=\underset{\mathbf{\varphi },\beta >0}{\max }\underset{\mathbf{\theta },\alpha >0}{\min }\mathcal{L}(\mathbf{\theta },\alpha ;\mathbf{\varphi },\beta ).$$

$${\mathbf{\theta }}^{t+1}$$

$${\alpha }^{t+1}$$

$${\mathbf{\varphi }}^{t+1}$$

$${\beta }^{t+1}$$

$$F_0(\mathbf{\theta })=⟨\mathbf{\psi }(\mathbf{\theta })|{\mathbf{M}}_0|\mathbf{\psi }(\mathbf{\theta })⟩.$$

$$G(\mathbf{\varphi })=⟨\mathbf{\xi }(\mathbf{\varphi })|\mathbf{S}|\mathbf{\xi }(\mathbf{\varphi })⟩,$$

$$F(\mathbf{\theta },\mathbf{\varphi })=⟨\mathbf{\xi }(\mathbf{\varphi }),\mathbf{\psi }(\mathbf{\theta })|\mathbf{M}|\mathbf{\xi }(\mathbf{\varphi }),\mathbf{\psi }(\mathbf{\theta })⟩$$

$$\mathbf{M}=\sum _{m=1}^M{\mathbf{e}}_m{\mathbf{e}}_m^{\top }\otimes {\mathbf{M}}_m$$

$$\mathbf{U}(\mathbf{\theta })={\mathbf{W}}_{P+1}\prod _{p=1}^P\exp (-j{\theta }_p{\mathbf{G}}_p){\mathbf{W}}_p,$$

Full text

††preprint: APS/123-QED

I Introduction

Conic programs defined over sparse graphs are ubiquitous in physics, quantum information, machine learning, and engineering. Many of these problems naturally take the form of semidefinite programs (SDPs) or nonconvex quadratically-constrained quadratic programs (QCQPs) that can be lifted to SDPs. In many of these formulations, the number of constraints scales with the graph size, and thus, although each constraint involves a small subset of variables, the Lagrangian function contains numerous terms. Representative examples include finding the ground-state energy of an Ising spin system [4], computing graph-theoretic quantities such as the Lovász theta number [42], verifying deep neural network properties [53], and networked optimization in engineering, such as the optimal power flow problem [3, 38]. Although classical algorithms for solving SDPs, such as interior-point methods [8], are guaranteed to find a globally optimal solution, their computational cost may scale unfavorably with the problem dimension.

Variational quantum algorithms (VQAs) have emerged as a versatile and promising paradigm for leveraging noisy intermediate-scale quantum (NISQ) processors to solve complex optimization problems. By employing parameterized quantum circuits (PQCs) optimized via classical routines, these hybrid quantum-classical algorithms are designed to comply with the constraints of current quantum hardware. Nonetheless, their application to constrained optimization remains underdeveloped. Existing methodologies are mostly limited to problems with binary variables and/or linear constraints, leaving a gap for the more general yet powerful class of conic programs, such as QCQPs and SDPs. Recently, Reference [51] used two PQCs to encode both primal and dual variables (Lagrange multipliers) for solving SDPs with many constraints. Nevertheless, it assumes efficient measurement of observables and relies on a dual-decomposition method that entails a nested variational optimization.

I.1 Contributions

This work proposes a novel VQA framework for solving large-scale conic programs with numerous constraints defined over sparsely connected graphs. The contributions are on three fronts:

c1) Optimizing a doubly parameterized Lagrangian function. Building on the modeling proposed in [51], where primal variables are captured by the state of a primal PQC and the dual variables are captured through the probability mass function (PMF) associated with a dual PQC, we develop a VQA seeking a saddle point of the associated Lagrangian function. We aim to minimize the Lagrangian function with respect to the parameters of the primal PQC and maximize it with respect to the parameters of the dual PQC. This is accomplished in a hybrid fashion. First, the two PQCs are used to measure gradients of the Lagrangian with respect to the PQC parameters. A classical computer subsequently uses gradient information to iteratively update PQC parameters using a primal-dual method. In contrast to the dual decomposition [51, 24], which requires solving an inner variational optimization problem to optimality at each iteration, our method performs a single gradient descent step on the primal variables and a single gradient ascent step on the dual variables per iteration. The proposed primal-dual updates apply directly to a broader class of optimization problems that can be solved via VQAs. Moreover, the proposed method enables the training of quantum machine learning (QML) models that satisfy constraints. Such functionality is essential for endowing QML models with safety and stability features.

c2) Efficient measurement of graph-induced observables. Although the concerned observables are defined over sparse Hermitian matrices, standard LCU (linear combination of unitaries)-based decomposition methods do not yield qubit-efficient measurement protocols. To address this challenge, we utilize the graph structure of the problem instance and the extended Bell measurements (XBM) [33], to design a qubit-efficient protocol for measuring gradients of the Lagrangian function associated with the QCQP/SDP. Specifically, the XBM method partitions the entries of an $N\times N$ Hermitian measurement matrix into groups, allowing entries within each group to be measured simultaneously using at most $\mathcal{O}(\log N)$ additional quantum gates. If a matrix can be decomposed into $\tilde{\mathcal{O}}(\log N)$ such groups, the XBM method is qubit-efficient. Hence, we propose to permute the graph’s nodes so that the related observables are banded. Finding the node permutation that yields the optimal number of groups is an NP-hard problem [18]. Instead, we adopt the reverse Cuthill-McKee (RCM) algorithm [14], a greedy heuristic algorithm of linear complexity $\mathcal{O}(N)$. This is effective because banded matrices are known to feature a small number of XBM groups [33]. The devised workflow can be instrumental in effectively measuring quantum observables associated with sparse connectivity graphs arising in applications over natural or engineered networked systems, including the representative examples discussed earlier.

c3) Application to the Optimal Power Flow (OPF) problem. We illustrate the efficacy of the proposed framework on the OPF problem. OPF can be posed as an optimization problem over thousands of variables and constraints. If power flows need to be modeled precisely through the alternating-current (AC) power flow equations, the OPF can be formulated as a QCQP. Although this QCQP is not convex, the OPF can be relaxed to a convex SDP and solved to global optimality under certain conditions [3, 38]. Nonetheless, the complexity of modern-day electric power grids raises scalability and optimality challenges. In this context, VQA offers a promising alternative, provided that the required quantum measurements can be performed efficiently. Two key points for successfully adopting the XBM method to the OPF are: i) The sparsity pattern (the positions of nonzero entries) of OPF observables is determined by the graph defined by the nodes and transmission lines of the power system, which is sparse; and ii) According to numerical tests, nodes in a power system can be permuted via the RCM algorithm so that OPF observables can be partitioned in $\mathcal{O}((\log N)^{3})$ groups.

I.2 Literature review

Handling constrained optimization problems using VQAs. A well-studied exemplar of VQAs is the variational quantum eigensolver (VQE), which aims to find the smallest eigenvalue of a quantum observable by employing a PQC [52]. The quantum approximate optimization algorithm (QAOA) is a variant of VQE with a problem-dependent PQC and a diagonal observable targeting combinatorial problems [17]. While VQE and QAOA were originally devised for unconstrained problems, extensions to optimization problems with constraints have been proposed in [43, 24, 39, 51, 11, 25, 6, 7]. References [43, 11] convert constrained problems to unconstrained ones by penalizing constraint violations on the objective function. However, selecting a suitable penalty parameter is a nontrivial task. Unreasonably large values can cause ill-conditioning, while unreasonably small values can yield infeasible solutions. In [25], the QAOA’s ansatz is adapted to ensure that the output state stays within the feasible subspace; yet, this strategy applies only to binary problems with linear constraints. References [6] propose a hybrid method for solving SDPs. The idea is to approximate the matrix variable of the original SDP by a semidefinite matrix of smaller dimension. The resultant SDP of reduced dimension is solved classically using standard interior point-based methods. Quantum measurements are used to translate the cost and constraint functions of the original to the compressed SDP formulation. Nonetheless, the solution of the compressed SDP may be infeasible or suboptimal for the original SDP.

An alternative approach to cope with constrained problems is via Lagrangian duality. Reference [24] seeks a saddle point of the Lagrangian function with respect to the PQC parameters and the dual variables. To tackle problems with numerous constraints, reference [51] encodes dual variables through a parameterized quantum state, in addition to the primal variables. Both [24] and [51] solve the associated dual problem using the iterative method of dual decomposition. Nevertheless, each iteration of dual decomposition involves solving a VQE-type problem to optimality, which is computationally overwhelming. In [39], this drawback is circumvented by using the primal-dual update method, which only requires estimating gradients of the Lagrangian function. Nonetheless, reference [39] applies only to binary and linear programs. In this work, we extend [39] and [51] to deal with a doubly variational VQA formulation of QCQP/SDPs, which we showcase via the OPF problem.

Measuring quantum observables. Measuring quantum observables is key to the effective implementation of VQAs. A prominent measurement method decomposes a Hermitian matrix into a linear combination of unitaries (LCU), and measures the resultant unitary observables via the destructive swap test [59, 21]. This method is practically relevant only if the resultant unitaries can be implemented efficiently using quantum gates. This is the case, for example, if the unitaries can be expressed as Pauli strings (tensor products of Pauli matrices), and the number of Pauli strings scales polynomially with the number of qubits, $n$; see [46, 12, 49]. Unfortunately, a general matrix may involve $\mathcal{O}(4^{n})$ Pauli strings [36]. The number of Pauli strings grows exponentially with $n$, even for some structured sparse matrices, such as $k$-banded matrices that involve $\mathcal{O}(kn2^{n})$ Pauli strings [36]. One possible strategy to reduce the number of measurements is to group Pauli strings into sets of terms that pairwise commute. If two Hermitian matrices commute, they can be jointly diagonalized, and thus, measured simultaneously [60, pg. 43]. To this end, various grouping methods have been proposed, including qubit-wise commutativity grouping [62, 63], general commutativity grouping [23], and unitary partitioning [29, 64]. Nevertheless, partitioning Pauli strings to achieve the minimum number of groups is an NP-hard problem [23].

Classical shadowing is an alternative technique for measuring observables [28]. It is particularly effective when the measurement matrix can be expressed as local Pauli strings, i.e., Pauli strings that act non-trivially on a small number of qubits. For observables that do not admit this form, however, the number of measurements scales linearly with the squared Frobenius norm of the matrix. Unfortunately, for connectivity matrices of sparse graphs, this quantity generally grows linearly with the network size.

Another recent method for measuring observables involves extended Bell measurements [33]. This method decomposes the measurement matrix into matrices whose associated observables are simultaneously measurable using qubit-efficient rotation unitaries. Nevertheless, the number of induced observables may reach $2^{n}$ for a general sparse measurement matrix [33]. Our key idea is to permute nodes of the underlying graph once and beforehand, so that the associated quantum observables can be measured efficiently.

Quantum computing for the OPF. A previous quantum computing effort to the OPF problem has focused on the classical primal-dual interior point method, wherein the Harrow–Hassidim–Lloyd (HHL) algorithm serves as a linear system solver in each Newton step [1]. The idea of utilizing the HHL algorithm as a linear system solver has also been explored for solving the nonlinear AC power flow equations, which constitute a subset of constraints of the AC OPF problem [54, 19, 41]. Nonetheless, the computational complexity of the HHL algorithm does not scale well for linear systems involving the Jacobian matrix of the AC power flow equations, as established in [50]. Reference [27] integrates a PQC as an intermediate layer of a classical neural network, trained in a supervised fashion to predict OPF solutions. However, it is not clear which measurement outputs are encoded into the subsequent classical hidden layer, nor how the parameter-shift rule is applied under such a design.

I.3 Notation

Column vectors are denoted by boldface lower-case letters. Matrices are denoted by boldface upper-case letters. The symbol ⊤ stands for transposition, and † for conjugate transposition. The operator $\operatorname{dg}(\mathbf{x})$ returns a diagonal matrix with vector $\mathbf{x}$ on its main diagonal, while $\operatorname{dg}(\mathbf{X})$ returns a vector with the diagonal entries of matrix $\mathbf{X}$. The set of real-valued symmetric matrices of dimension $M$ is denoted by $\mathbb{S}^{M}$. Logarithmic functions are of the binary base. The imaginary unit is denoted by $\iota\coloneqq\sqrt{-1}$. The notation $m=1:M$ means that index $m$ takes values from $1$ to $M$. Symbol $\left\|\mathbf{H}\right\|$ denotes the spectral norm of a matrix $\mathbf{H}$.

II Primal-dual iterations for QCQP

We focus on optimization problems defined over sparse graphs, often arising from physical models as nonconvex QCQPs or SDPs. The latter may appear either as native formulations or as convex relaxations of the QCQPs. Prototypical examples of such problems include:

• •

Ising ground-state energy: Minimizing a quadratic Hamiltonian $H(\mathbf{s})=\sum_{i<j}J_{ij}s_{i}s_{j}+\sum_{i}h_{i}s_{i}$ over spin configurations $\mathbf{s}\in\{\pm 1\}^{N}$ can be expressed as a nonconvex QCQP. The coupling matrix $\mathbf{J}$ is typically sparse, reflecting local interactions in physical lattices or graphical models. The resulting QCQP admits an SDP relaxation obtained by lifting the quadratic $\pm 1$ constraints $s_{i}^{2}=1$; see [4].

• •

Lovász theta number: For a graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with $N=|\mathcal{V}|$ vertices, the Lovász theta number can be computed as the optimal value of an SDP over an $N\times N$ matrix variable:

[TABLE]

where $\mathbf{1}$ is the all-ones vector. Each edge $(i,j)\in\mathcal{E}$ yields one constraint, so the constraint matrices inherit the sparsity pattern of $\mathcal{G}$. While commonly presented as an SDP, the Lovász theta number admits an equivalent QCQP formulation [42].

• •

Neural network verification: Certifying the robustness of neural networks (NN) against adversarial perturbations entails bounding the worst‑case margin between classes over a permissible input set. In networks using the rectified linear unit (ReLU) as activation function, the input-output relationship of the ReLU $z=\max(0,x)$ admits an exact reformulation including a quadratic constraint, namely $z\geq 0$, $z\geq x$, and $z(z-x)=0$ [53]. Based on this reformulation, the NN verification task can be cast as a nonconvex QCQP over the high-dimensional vector of neuron activations. The constraint matrices are highly sparse and structured, reflecting the layered connectivity of the underlying computational graph [53].

• •

Optimal power flow: The AC OPF problem minimizes the electricity generation cost to meet load demand while satisfying grid engineering constraints. The problem can be posed as a nonconvex QCQP over nodal voltages, where the quadratic constraints represent Kirchhoff’s laws and operational limits. The sparsity pattern is governed by the power grid topology through transmission lines, as each grid node is directly connected to only a few other nodes [3].

To streamline the presentation, we focus on the QCQP formulation and discuss its extension to the SDP setting in Section III. We consider the QCQP:

[TABLE]

where $\{\mathbf{M}_{m}\in\mathbb{C}^{N\times N}\}_{m=0}^{M}$ are sparse Hermitian matrices determined by the underlying problem graph, and the number of constraints $M$ typically scales with the network size $N$.

Although (QCQP) is generally nonconvex, it can be relaxed to an SDP, which can be solved by standard interior-point solvers. In some application domains, this convex relaxation is successful in the sense that the optimal SDP variable is rank-one and a globally optimal $\mathbf{x}$ can be readily recovered from the SDP solution. Nonetheless, the memory and time requirements imposed by SDP solvers may be challenging to meet as $N$ increases. This motivates us to explore a VQA for QCQP/SDPs.

Toward developing a PQC-based solver for large-scale QCQP/SDPs, let us discuss how (QCQP) could be solved on a classical computer. Among the different approaches, we present the primal-dual (PD) decomposition as a representative example of a first-order optimization algorithm. The PD method relies on Lagrangian duality [8]. Each constraint indexed by $m$ in (QCQP) is associated with a non-negative dual variable $\lambda_{m}$. Let vector $\boldsymbol{\lambda}\in\mathbb{R}^{M}_{+}$ collect all dual variables. The related Lagrangian function is a function of both $\mathbf{x}$ and $\boldsymbol{\lambda}$:

[TABLE]

The corresponding dual problem is posed as

[TABLE]

A primal/dual solution $(\mathbf{x}^{*},\boldsymbol{\lambda}^{*})\in\mathbb{C}^{N}\times\mathbb{R}^{M}_{+}$ of (2) satisfies

[TABLE]

for all $(\mathbf{x},\boldsymbol{\lambda})\in\mathbb{C}^{N}\times\mathbb{R}^{M}_{+}$. Such a pair is called a saddle point of the Lagrangian function. In essence, $\mathbf{x}^{*}$ is a minimizer of $\mathcal{L}(\mathbf{x};\boldsymbol{\lambda}^{*})$ and $\boldsymbol{\lambda}^{*}$ is a maximizer of $\mathcal{L}(\mathbf{x}^{*};\boldsymbol{\lambda})$. For problems featuring strong duality (in which case $P^{*}=D^{*}$), a point is an optimal primal/dual point if and only if it is a saddle point of the Lagrangian function [8, p. 239]. Furthermore, at a primal/dual solution $(\mathbf{x}^{*};\boldsymbol{\lambda}^{*})$, the value of the Lagrangian satisfies $\mathcal{L}(\mathbf{x}^{*};\boldsymbol{\lambda}^{*})=P^{*}=D^{*}$. Finding an exact saddle point of the Lagrangian imposes stringent conditions on $\mathcal{L}(\mathbf{x};\boldsymbol{\lambda})$. However, approximate saddle points can be reached under specific settings using gradient-based methods, such as the PD method [48].

According to the PD method, during iteration $t$, the primal and dual variables are updated as:

[TABLE]

where $\mu_{x}^{t}$ and $\mu_{\lambda}^{t}$ are positive step sizes, and $[\lambda]_{+}\coloneqq\max\{\lambda,0\}$. The update in (3a) is a gradient descent step on $\mathcal{L}$ in terms of the primal variable. The update in (3b) is a projected gradient ascent step on $\mathcal{L}$ in terms of the dual variable.

One may attempt to solve (QCQP) using the PD iterations in (3). In this case, computing the Lagrangian gradients incurs complexity $\mathcal{O}(sN^{2})$, assuming each row of $\mathbf{M}_{m}$ has at most $s$ nonzero entries. To see this, note that the partial derivative of $\mathcal{L}$ with respect to $\lambda_{m}$ is

[TABLE]

which can be computed in $\mathcal{O}(sN)$ operations. Therefore, computing the entire gradient $\nabla_{\boldsymbol{\lambda}}\mathcal{L}$ entails complexity $\mathcal{O}(sNM)=\mathcal{O}(sN^{2})$ given that $M$ is proportional to $N$. Similarly, computing the gradient

[TABLE]

incurs $\mathcal{O}(sN^{2})$ operations. For large-scale problems, such as solving AC OPF for contemporary power systems with more than $N=50,000$ nodes, this computational cost quickly becomes prohibitive. In this context, can a VQA attain reduced computational complexity? The next section proposes a variational model for QCQP/SDPs.

III A doubly variational QCQP/SDP

Variational quantum computing is an algorithmic tool for solving high-dimensional problems in a parameterized form using a PQC. The well-known VQE solves a specific form of (QCQP) [52]:

[TABLE]

The optimal cost in (4) is the smallest eigenvalue of $\mathbf{M}_{0}$. Moreover, the eigenvector corresponding to the smallest eigenvalue of $\mathbf{M}_{0}$ is the minimizer. The key idea of VQE is to model the original high-dimensional vector $\mathbf{x}$ as a quantum state $\ket{\boldsymbol{\psi}}$ using $\log N$ qubits. Suppose for now that $N$ is a power of 2. The constraint in (4) can be omitted as the quantum state has unit Euclidean norm anyway. This quantum state is generated upon applying a PQC on an initial state $\ket{\mathbf{0}}_{\log N}$. The PQC is modeled by a unitary matrix $\mathbf{U}(\boldsymbol{\theta})$ parameterized by vector $\boldsymbol{\theta}\in\mathbb{R}^{P}$ with $P\ll N$. The PQC generates the state $\ket{\boldsymbol{\psi}(\boldsymbol{\theta})}=\mathbf{U}(\boldsymbol{\theta})\ket{\mathbf{0}}$, and then measures the expectation $F_{0}(\boldsymbol{\theta})=\braket{\boldsymbol{\psi}(\boldsymbol{\theta})|\mathbf{M}_{0}|\boldsymbol{\psi}(\boldsymbol{\theta})}$ and its gradient $\nabla_{\boldsymbol{\theta}}F_{0}(\boldsymbol{\theta})$. Subsequently, a classical computer receives the gradient information and updates $\boldsymbol{\theta}$ to minimize $F_{0}(\boldsymbol{\theta})$ using gradient descent or other classical optimization methods.

We propose solving the (QCQP) in a similar manner. Because the primal variable of the (QCQP) does not have a unit Euclidean norm, we parameterize it as [57]:

[TABLE]

where $\alpha>0$ is an auxiliary optimization variable. The PQC used to generate $\ket{\boldsymbol{\psi}(\boldsymbol{\theta})}$ will be henceforth termed the primal PQC or PQCp.

Given the parameterization in (5), the problem (QCQP) can be expressed in the variational form:

[TABLE]

where expectations are defined as

[TABLE]

These expectations can be alternatively expressed via the density operator $\boldsymbol{\rho}(\boldsymbol{\theta})$ associated with the PQCp state as

[TABLE]

If PQCp is characterized by a pure state $\ket{\boldsymbol{\psi}(\boldsymbol{\theta})}$, its density operator is $\boldsymbol{\rho}(\boldsymbol{\theta})=\ket{\boldsymbol{\psi}(\boldsymbol{\theta})}\!\bra{\boldsymbol{\psi}(\boldsymbol{\theta})}$. Generally, a density operator is a positive semidefinite (PSD) and unit-trace matrix [49]. Under this interpretation, (QCQPθ) is a variational SDP over the density operator $\boldsymbol{\rho}(\boldsymbol{\theta})$.

As in VQE, the original (QCQP) over $\mathbf{x}$ has been replaced by the variational problem (QCQPθ) over the parameter vector $\boldsymbol{\theta}$. Attempting to solve (QCQPθ) using Lagrangian duality entails solving the dual problem:

[TABLE]

over $\boldsymbol{\lambda}\in\mathbb{R}^{M}_{+}$. The related Lagrangian function is

[TABLE]

Unfortunately, because the dual variable $\boldsymbol{\lambda}$ has length $\mathcal{O}(N)$, the curse of dimensionality remains.

To bypass this difficulty, dual variables can also be represented using a variational model. A similar idea was proposed in [51] for capturing the primal and dual variables of an SDP. To model $\boldsymbol{\lambda}$ variationally, we introduce a second PQC, henceforth termed the dual PQC or PQCd. The dual PQC operates on $\log M$ qubits, supposing again that the number of constraints, $M$, is a power of two. The dual PQC is modeled by a unitary matrix $\mathbf{V}(\boldsymbol{\phi})$ parameterized by vector $\boldsymbol{\phi}\in\mathbb{R}^{Q}$ such that $Q\ll M$. The state of PQCd is $\ket{\boldsymbol{\xi}(\boldsymbol{\phi})}=\mathbf{V}(\boldsymbol{\phi})\ket{\mathbf{0}}$; see Fig. 1.

With the aid of PQCd, the $m$-th entry of $\boldsymbol{\lambda}$ is parameterized as

[TABLE]

where $\beta>0$ is an auxiliary variable and $\xi_{m}(\boldsymbol{\phi})$ is the $m$-th entry of the state $\ket{\boldsymbol{\xi}(\boldsymbol{\phi})}$. Like $\alpha$ in (5), variable $\beta$ is introduced for scaling purposes. Unlike (5), however, dual variables are related to the squared magnitudes of quantum state entries. This ensures that $\lambda_{m}(\boldsymbol{\phi})$ takes real nonnegative values.

The Lagrangian function related to this doubly parameterized QCQP can be expanded into three terms:

[TABLE]

where

[TABLE]

and

[TABLE]

The variational dual problem is defined as

[TABLE]

One may wonder why we do not parameterize the primal and dual variables jointly using a single PQC operating on $\log N+\log M$ qubits. Unfortunately, under that design, we would not be able to update the PQC parameters in a meaningful way. Having two PQCs parameterized separately, as depicted in Figure 1, allows us to minimize the Lagrangian function in (10) over $(\boldsymbol{\theta},\alpha)$ and maximize it over $(\boldsymbol{\phi},\beta)$.

The standard Lagrangian dual function is known to be concave with respect to the dual variables, even if the primal problem is nonconvex. That is not the case for the dual function $\min_{\boldsymbol{\theta},\alpha>0}\mathcal{L}(\boldsymbol{\theta},\alpha;\boldsymbol{\phi},\beta)$ in (13). This is because the dual variables are now parameterized, and $\mathcal{L}(\boldsymbol{\theta},\alpha;\boldsymbol{\phi},\beta)$ is not concave in $(\boldsymbol{\phi},\beta)$.

IV Solving the doubly variational QCQP/SDP

We propose seeking a saddle point of the Lagrangian function in (10) using a PD algorithm. The $t$-th iteration of the algorithm comprises four steps:

[TABLE]

where $\mu^{t}_{\theta}$, $\mu^{t}_{\alpha}$, $\mu^{t}_{\phi}$, and $\mu^{t}_{\beta}$ are positive step sizes. Steps (14a)–(14b) constitute gradient descent steps to update the primal variables, while (14c)–(14d) are gradient ascent steps to update the dual variables. The updates of (14) run on a classical computer. The suggested workflow is shown in Fig. 1.

The gradients in (14) can be estimated with the help of the two PQCs. The key observation is that each term of the Lagrangian function can be expressed as an expectation over an observable. More specifically, the first term can be computed using an observable operating on PQCp through the expectation:

[TABLE]

The third term can be computed using an observable operating on PQCd through the expectation:

[TABLE]

where $\mathbf{S}\coloneqq\operatorname{dg}(\{b_{m}\}_{m=1}^{M})$ is a diagonal matrix with the OPF parameters $b_{m}$ on its main diagonal. Due to the parameterization in (9), the Lagrangian of the doubly parameterized problem is quadratic in $\ket{\boldsymbol{\xi}(\boldsymbol{\phi})}$, although the Lagrangian in (8) was linear in $\boldsymbol{\lambda}$.

The second term of the Lagrangian can be expressed as the expectation over an observable operating jointly on PQCp and PQCd, as shown in the ensuing lemma proved in Appendix A.2.

Lemma 1**.**

The second summand in (11) can be computed as the expectation

[TABLE]

defined by the $MN\times MN$ Hermitian matrix

[TABLE]

and $\mathbf{e}_{m}$ is the $m$-th column of $\mathbf{I}_{M}$.

Lemma 1 asserts that $F(\boldsymbol{\theta},\boldsymbol{\phi})$ is an expectation applied to the composite state of the two PQCs. This is important as it bypasses the need to measure each one of the numerous constraint expectations $\{F_{m}(\boldsymbol{\theta})\}_{m=1}^{M}$ of the OPF separately. This powerful feature was first identified in [51] for solving doubly variational SDPs on a quantum computer. In this work, we adopt this idea to the setting of solving QCQPs and SDPs over sparse graphs. We also propose a method for efficiently measuring the observables associated with sparse graphs. We uniquely leverage the underlying graph structure to efficiently measure the three observables using the two PQCs. We defer this discussion to Section V. If the observables can be measured efficiently, their gradients can be measured efficiently as well via the so-called parameter-shift rule (PSR). We review this rule and adapt it to the PD iterations in (14).

Let us focus on the gradient $\nabla_{\boldsymbol{\theta}}F_{0}(\boldsymbol{\theta}^{t})$. Suppose that the unitary matrix modeling PQCp takes the form

[TABLE]

where each $\mathbf{G}_{p}$ is a single-qubit Hermitian generator with two distinct eigenvalues $\pm r$, and $\{\mathbf{W}_{p}\}_{p=1}^{P+1}$ is a set of fixed gates. Then, the partial derivative of $F_{0}(\boldsymbol{\theta}^{t})$ with respect to the $p$-th entry of $\boldsymbol{\theta}$ can be computed by evaluating $F_{0}(\boldsymbol{\theta})$ at two values of $\boldsymbol{\theta}$ [47, 58]:

[TABLE]

where $\mathbf{e}_{p}$ is the $p$-th column of the identity matrix $\mathbf{I}_{P}$. If $\mathbf{G}_{p}$ is a Pauli rotation from the set $\frac{1}{2}\{\sigma_{x},\sigma_{y},\sigma_{z}\}$, we get $r=\frac{1}{2}$ and each $\theta^{t}_{p}$ is shifted by $\pm\frac{\pi}{2}$.

The expectation value $F_{0}(\boldsymbol{\theta})$ cannot be evaluated exactly; it can only be estimated via quantum measurements: To estimate $F_{0}(\boldsymbol{\theta})$ for a particular $\boldsymbol{\theta}$ with precision $\epsilon$, PQCp has to be sampled $S=\mathcal{O}(\epsilon^{-2})$ times. To obtain the complete gradient $\nabla_{\boldsymbol{\theta}}F_{0}(\boldsymbol{\theta}^{t})$, the previous process is repeated $2P$ times and involves a total of $2PS$ measurement samples.

The gradient of the expectation value in (16) with respect to $\boldsymbol{\phi}$ can be computed similarly by running PQCd by shifting the values of $\boldsymbol{\phi}^{t}$ for $2Q$ times or $2QS$ measurement samples. The gradient $\nabla_{\boldsymbol{\theta}}F(\boldsymbol{\theta}^{t},\boldsymbol{\phi}^{t})$ can be computed by running both PQCs while only shifting the entries of $\boldsymbol{\theta}^{t}$; and vice versa for computing the gradient $\nabla_{\boldsymbol{\phi}}F(\boldsymbol{\theta}^{t},\boldsymbol{\phi}^{t})$. It should be clear by now that, due to the PSR, evaluating gradients of expectations is equivalent to measuring expectations. Section V discusses how to measure the observables induced from a sparse graph.

Although the plain PD method has been extensively studied for convex/concave and nonconvex/concave saddle-point problems [31, 61, 35, 40], its convergence is not fully understood in the nonconvex/nonconcave setting, like the one in (13). Recently, a variation of the PD method, termed the extragradient (EG) method [34], has been shown to converge under certain conditions for nonconvex/nonconcave saddle-point problems [16]. Compared to the plain PD method, each EG iteration doubles the number of gradient evaluations of the Lagrangian function. More specifically, let vector $\mathbf{z}^{t}\coloneqq[\boldsymbol{\theta}^{t};\alpha^{t};\boldsymbol{\phi}^{t};\beta^{t}]$ collect the primal and dual variables at iterate $t$, and define the vector:

[TABLE]

Then, the EG iterations can be expressed as

[TABLE]

for a step size $\mu_{z}^{t}>0$. In essence, the EG method first computes an intermediate point $\bar{\mathbf{z}}^{t}$ based on the plain PD method, and then updates $\mathbf{z}^{t}$ to $\mathbf{z}^{t+1}$ upon evaluating the gradient operator at the intermediate point. Despite doubling the number of gradient evaluations, EG iterates feature favorable convergence properties as established analytically in Section VI and corroborated numerically in Section VIII.

Before concluding this section, it is worth noting that, as with other VQAs, training the primal-dual PQC pair could be impeded by the so-called barren plateaus issue, according to which gradients vanish exponentially fast in terms of the number of qubits and depth of PQCs [45, 37]. This phenomenon can occur when designing VQAs for large-scale problems and/or using overly expressive PQCs. A related concern is that the PQC architectures that are both sufficiently expressive and free from barren plateaus are often efficiently simulated by classical algorithms [10]. Interestingly, recent results show that the so-called dynamic PQCs [15], which incorporate early measurements followed by measurement-conditional unitaries, and quantum recurrent embedding neural networks [30] can provably avoid barren plateaus while retaining high expressivity. Moreover, dynamic PQCs inherit the worst-case classical hardness of the universal unitary circuits, precluding the possibility of classical simulatability. Investigating the performance of the doubly variational QCQP/SDP on dynamic PQCs or quantum recurrent embedding neural networks for large-scale problems is an interesting future direction.

V Measuring Observables Induced by Sparse Graphs

This section reviews the extended Bell measurement (XBM) method proposed in [33] and explains how it can be adapted to efficiently measure observables induced from sparse graphs.

V.1 Extended Bell measurement (XBM) method

Consider measuring the expectation value $\bra{\boldsymbol{\psi}(\boldsymbol{\theta})}\mathbf{M}_{m}\ket{\boldsymbol{\psi}(\boldsymbol{\theta})}$ associated with the $N\times N$ Hermitian matrix $\mathbf{M}_{m}$. The matrix $\mathbf{M}_{m}$ is conventionally expressed as

[TABLE]

The XBM method proposes the alternative expansion:

[TABLE]

where $\oplus$ denotes the bitwise exclusive OR (XOR) operator if indices $i$ and $c$ are represented in binary form. To establish the equivalence between (23) and (24), we need to show that every pair $(i,j)$ corresponds to a unique pair $(i,i\oplus c)$, which indeed holds for $c=i\oplus j$. Due to this alternative indexing, the matrix $\mathbf{M}_{m}$ can be decomposed into $N$ matrices $\{\mathbf{M}_{m}^{c}\}_{c=0}^{N-1}$ defined as

[TABLE]

For example, the matrix $\mathbf{M}_{m}^{0}$ carries the diagonal entries of $\mathbf{M}_{m}$ because $i\oplus i=0$. Figure 2 uses color coding to identify the index pairs $(i,j)$ corresponding to the same $c$ for an observable over three qubits [33]. We will henceforth say that all index pairs $(i,j)$ satisfying $i\oplus j=c$ for a specific $c$ belong to the same color $c$.

Because $i\oplus j=j\oplus i$, every $\mathbf{M}_{m}^{c}$ is Hermitian. Therefore, the matrices

[TABLE]

are also Hermitian for all $c$, and $\check{\mathbf{M}}_{m}^{0}=\mathbf{0}$.

Based on (24), the expectation value associated with $\mathbf{M}_{m}$ can be expressed as a sum of $2N-1$ observables:

[TABLE]

Because $\hat{\mathbf{M}}_{m}^{0}$ is diagonal, the expectation value $\braket{\boldsymbol{\psi}|\hat{\mathbf{M}}_{m}^{0}|\boldsymbol{\psi}}$ can be measured in the computational basis. Upon measuring $\ket{\boldsymbol{\psi}}$, the binary vector outcome $\ket{i}$ is observed with probability $|\psi_{i}|^{2}$. Define a random variable taking the value $\braket{i|\hat{\mathbf{M}}_{m}^{0}|i}$ when outcome $\ket{i}$ is observed. Then, the expectation $\braket{\boldsymbol{\psi}|\hat{\mathbf{M}}_{m}^{0}|\boldsymbol{\psi}}$ in (26) is the expected value of this random variable.

The remaining expectation values in (26) are not as straightforward to measure because the matrices $(\hat{\mathbf{M}}_{m}^{c},\check{\mathbf{M}}_{m}^{c})$ are non-diagonal for $c=1,\ldots,N-1$. The advantage of XBM’s color decomposition is that these matrices can be diagonalized by unitaries corresponding to qubit-efficient circuits. This key result from [33] is summarized in the following lemma.

Lemma 2** ([33]).**

Consider the color index $c\geq 1$. Let index $k_{c}$ denote the most significant bit taking the value of one in the binary representation of $\ket{c}=\ket{c_{\log N-1}\cdots c_{1}c_{0}}$. Then, the matrices $\hat{\mathbf{M}}_{m}^{c}$ and $\check{\mathbf{M}}_{m}^{c}$ admit the eigenvalue decompositions:

[TABLE]

where

• i)

unitary matrix $\mathbf{U}_{c}$ can be implemented using one Hadamard gate and at most $(\log N-1)$ CNOT gates;

• ii)

unitary matrix $\mathbf{S}_{c}$ applies a phase gate $\mathbf{S}$ on qubit $k_{c}$ and acts trivially on all other qubits; and

• iii)

the diagonal eigenvalue matrices are

[TABLE]

with $\mathbf{X}_{k_{c}}$ applying a NOT gate on qubit $k_{c}$.

The practical value of Lemma 2 is that the expectation value

[TABLE]

can be measured by transforming $\ket{\boldsymbol{\psi}}$ to $\mathbf{U}_{c}^{\top}\ket{\boldsymbol{\psi}}$ and then measuring the latter in the computational basis using the diagonal matrix $\hat{\mathbf{\Lambda}}_{m}^{c}$. The transformation is qubit-efficient because the unitary $\mathbf{U}_{c}$ involves at most $\log N$ elementary gates. The expectation value $\braket{\boldsymbol{\psi}|\check{\mathbf{M}}_{m}^{c}|\boldsymbol{\psi}}$ can be measured similarly, with at most $\log N+1$ elementary gates; $\log N$ for $\mathbf{U}_{c}$ plus one for $\mathbf{S}_{c}$. A critical feature is that the unitary $\mathbf{U}_{c}$ diagonalizing $(\hat{\mathbf{M}}_{m}^{c},\check{\mathbf{M}}_{m}^{c})$ depends solely on the color (locations of nonzero entries) of $(\hat{\mathbf{M}}_{m}^{c},\check{\mathbf{M}}_{m}^{c})$, and not on the values of those entries. Only the eigenvalue matrices $(\hat{\mathbf{\Lambda}}_{m}^{c},\check{\mathbf{\Lambda}}_{m}^{c})$ depend on those values, and can be computed without the need for an eigenvalue decomposition.

Although each expectation value in (26) can be estimated efficiently, there are $(2N-1)$ of them. Unfortunately, each one of them requires a different quantum circuit. Therefore, the expectation value $\braket{\boldsymbol{\psi}|\mathbf{M}_{m}|\boldsymbol{\psi}}$ can be measured efficiently only if all nonzero entries of $\mathbf{M}_{m}$ are covered by a few colors. Let $\mathcal{C}$ denote the set of colors for which either $\hat{\mathbf{M}}_{m}^{c}$ or $\check{\mathbf{M}}_{m}^{c}$ is nonzero, so that $\mathbf{M}_{m}$ is decomposed as $\mathbf{M}_{m}=\sum_{c\in\mathcal{C}}\mathbf{M}_{m}^{c}$. For the XBM protocol to be qubit-efficient, the number of colors $C=|\mathcal{C}|$ should scale polynomially in $\log N$. Reference [33] establishes that banded matrices occupy relatively few colors. Specifically, it is shown that a $k$-banded matrix can be decomposed using $C=\mathcal{O}(k\log N)$ rather than $N$ colors. For the example of Figure 2, a matrix with bandwidth $k=1$ requires 4 colors.

Based on Lemma 2, if the expectation value $F_{m}(\boldsymbol{\theta})$ occupies $C$ colors, it can be decomposed into $2C-1$ expectations:

[TABLE]

where $\ket{\boldsymbol{\psi}_{c}(\boldsymbol{\theta})}=\mathbf{U}_{c}^{\top}\ket{\boldsymbol{\psi}(\boldsymbol{\theta})}$ and $\mathbf{\Lambda}_{m}^{c}=\hat{\mathbf{\Lambda}}_{m}^{c}$ for $c=1,\ldots,C$; and $\ket{\boldsymbol{\psi}_{c}(\boldsymbol{\theta})}=\mathbf{U}_{c}^{\top}\mathbf{S}_{c}\ket{\boldsymbol{\psi}(\boldsymbol{\theta})}$ and $\mathbf{\Lambda}_{m}^{c}=\check{\mathbf{\Lambda}}_{m}^{c}$ for $c=C+1,\ldots,2C-1$. Each one of the observables involved can be diagonalized by a different unitary matrix.

We next describe how to effect a measurement matrix $\mathbf{M}$ that uses a few colors.

V.2 Node permutation for qubit-efficient measurements

To solve the doubly parameterized QCQP/SDP, we need to evaluate the Lagrangian function in (10) by measuring the observable in Lemma 1 so that we can compute its gradients over $(\boldsymbol{\theta},\alpha;\boldsymbol{\phi},\beta)$. We next adopt the XBM protocol to efficiently measure the related observables in (10).

According to (18), observable matrix $\mathbf{M}$ has the $N\times N$ matrices $\{\mathbf{M}_{m}\}_{m=1}^{M}$ as its diagonal blocks. Despite being block-diagonal, matrix $\mathbf{M}$ may still occupy up to $N$ colors. Critically, we do not necessarily have to work with $\mathbf{M}$. Instead, we can permute variables to effect a new measurement matrix that is amenable to more efficient measurements. Suppose we permute vector $\mathbf{x}$ by a permutation matrix $\mathbf{P}$ to get $\mathbf{P}\mathbf{x}$. This entails reordering the nodes of the underlying graph, which is a trivial task that can be performed before solving (QCQP). The observable in Lemma 1 can be expressed as

[TABLE]

Ideally, we would like to design a single $\mathbf{P}$ to minimize the maximum bandwidth across all $\mathbf{P}\mathbf{M}_{m}\mathbf{P}^{\top}$. For many conic programs defined over sparse graphs, the constraint matrices $\mathbf{M}_{m}$ are dictated by a common underlying graph. In such problems, there often exists a representative matrix $\mathbf{Y}$, for instance, the graph Laplacian matrix, whose sparsity pattern captures the union of the sparsity patterns of all $\mathbf{M}_{m}$. Therefore, we propose designing $\mathbf{P}$ to minimize the bandwidth of $\mathbf{P}\mathbf{Y}\mathbf{P}^{\top}$. This is because minimizing the bandwidth of $\mathbf{P}\mathbf{Y}\mathbf{P}^{\top}$ is equivalent to minimizing the maximum bandwidth across all $\mathbf{P}\mathbf{M}_{m}\mathbf{P}^{\top}$. In the OPF problem (detailed in Section VII), the representative matrix is the node admittance matrix $\mathbf{Y}$ of the power network. We henceforth assume such a representative matrix $\mathbf{Y}$ is available or can be constructed for the problem at hand.

Unfortunately, minimizing the bandwidth of a matrix is an NP-hard problem even if the matrix is sparse [18]. Nonetheless, heuristic approaches, such as the reverse Cuthill-McKee (RCM) algorithm [14], provide reasonably good solutions in linear complexity $\mathcal{O}(L_{e})=\mathcal{O}(N)$, where $L_{e}$ is the number of edges in the underlying graph. The RCM algorithm has been used in classical computing to reduce the bandwidth of sparse matrices arising in various domains, including power systems, where it is applied to the node admittance matrix $\mathbf{Y}$ [55, 20, 13]. Albeit permuting $\mathbf{M}_{m}$’s using the RCM algorithm takes $\mathcal{O}(N)$ operations, the crucial point is that this permutation can be performed once for a given graph topology. In applications where the graph structure remains fixed while other problem parameters vary, the same permutation can be reused across multiple problem instances. For example, in the OPF problem, the power system topology changes infrequently and affects only a few edges of the graph. As long as the network topology remains unchanged, an instance of the OPF is characterized by parameters $b_{m}$, which vary frequently. In contrast to the permutation step, OPF instances over the same topology must be solved repeatedly, every few minutes during real-time operation and potentially millions of times in long-term planning.

We numerically validated the effect of node permutation on the bandwidth of $\mathbf{M}$. We ran the RCM algorithm on the 66 benchmark power system models of the pglib dataset [2]. Figure 3 shows the bandwidth of $\mathbf{M}$ before and after node permutation. Because the RCM solution depends on the initial node, we conducted 200 Monte Carlo runs per power system, with the initial node drawn uniformly at random per run. Out of the 200 runs, we retained the permutation that yielded the smallest bandwidth. Evidently, permuting nodes reduces the bandwidth of $\mathbf{M}$ dramatically. Upon data fitting and cross-validation, the reduced bandwidth was numerically found to scale as $(\log N)^{3}$. We also empirically tested the effect of node permutation on the number of colors in $\mathbf{M}$ for these power systems. Figure 4 shows that node permutation significantly reduces the number of colors, too. If the bandwidth of $\mathbf{P}\mathbf{M}\mathbf{P}^{\top}$ is $k=\mathcal{O}((\log N)^{3})$ per data fitting in Figure 3, the analysis of [33] upper bounds the number of its colors as $C=\mathcal{O}(k\log N)=\mathcal{O}((\log N)^{4})$. However, the number of colors in $\mathbf{P}\mathbf{M}\mathbf{P}^{\top}$ was found to scale as $(\log N)^{3}$.

We have so far focused on measuring the second observable of the Lagrangian function in (10). The third observable is diagonal and is thus straightforward to measure. The first observable $F_{0}(\boldsymbol{\theta})$, however, is defined over the non-diagonal matrix $\mathbf{M}_{0}$. Fortunately, the sparsity pattern of $\mathbf{M}_{0}$ is a subset of the sparsity pattern of $\mathbf{Y}$. Therefore, if nodes are permuted so that $\mathbf{M}$ occupies a few colors $C$, the same permutation works well also for $\mathbf{M}_{0}$. In other words, matrix $\mathbf{P}\mathbf{M}_{0}\mathbf{P}^{\top}$ would occupy $C$ colors or fewer.

Since permuting the underlying graph can yield measurement matrices with much fewer colors, we suggest the following workflow for solving doubly parameterized QCQPs over sparse graphs: i) Identify the underlying graph and construct a representative matrix $\mathbf{Y}$ whose sparsity pattern is the union of the sparsity patterns of all $\mathbf{M}_{m}$. For the OPF, $\mathbf{Y}$ is the node admittance matrix; ii) Feed the sparsity pattern of $\mathbf{Y}$ into the RCM algorithm to find a near-optimal node permutation matrix $\mathbf{P}$; iii) Formulate (QCQP) based on the permuted nodes; iv) Solve (QCQP) using the VQA of Section III; and v) If needed, apply the reverse permutation to the primal solution of (QCQP) to recover the original node ordering. It is henceforth assumed that primal variables and measurement matrices have been permuted already. Therefore, we will henceforth use the lighter notation $\mathbf{M}$ and $\mathbf{M}_{m}$ instead of $\mathbf{P}\mathbf{M}\mathbf{P}^{\top}$ and $\mathbf{P}\mathbf{M}_{m}\mathbf{P}^{\top}$.

V.3 Implementation details

We next explain how to estimate the Lagrangian gradients needed in (14) using the primal/dual PQC pair. We start with $\nabla_{\alpha}\mathcal{L}$ in (14b), which entails estimating $F_{0}(\boldsymbol{\theta})$ and $F(\boldsymbol{\theta},\boldsymbol{\phi})$. To simplify the presentation of measuring $F_{0}(\boldsymbol{\theta})$, suppose $\mathbf{M}_{0}$ occupies all $C$ colors used in $\mathbf{Y}$. According to (27), measuring $F_{0}(\boldsymbol{\theta})$ requires $2C-1$ circuits. Each rotated primal circuit $c$ involves the same primal PQC $\mathbf{U}(\boldsymbol{\theta})$ followed by a different XBM unitary $\mathbf{U}_{c}$. These $2C-1$ circuits can be executed sequentially or in parallel, as shown in Figure 5. The sequential implementation operates on $\log N$ qubits and requires compiling $2C-1$ circuits in a sequential fashion. The parallel implementation is $2C-1$ times faster but operates on $\log N(2C-1)$ qubits and requires replicating the primal PQC $2C-1$ times.

Measuring observable $F(\boldsymbol{\theta},\boldsymbol{\phi})$ is more complicated as it involves estimating the double summation

[TABLE]

Interestingly, for a specific index $c$, the corresponding sum $\sum_{m=1}^{M}|\xi_{m}(\boldsymbol{\phi})|^{2}F^{c}_{m}(\boldsymbol{\theta})$ can be evaluated in two steps. We first measure the dual PQC and observe $\ket{m}$ with probability $|\xi_{m}(\boldsymbol{\phi})|^{2}$. Given the observed $\ket{m}$, we can now easily estimate the expectation $F_{m}^{c}(\boldsymbol{\theta})$ corresponding to a diagonal observable.

Measuring $\nabla_{\beta}\mathcal{L}$ in (14d) is trivial as it depends on $F(\boldsymbol{\theta},\boldsymbol{\phi})$ and $G(\boldsymbol{\phi})$.

The gradient $\nabla_{\boldsymbol{\theta}}\mathcal{L}(\boldsymbol{\theta})$ can be estimated based on the PSR. For example, the partial derivative of $\mathcal{L}$ with respect to $\theta_{p}$ can be found as

[TABLE]

where $\boldsymbol{\theta}^{\pm}_{p}\coloneqq\boldsymbol{\theta}\pm\frac{\pi}{2}\mathbf{e}_{p}$. According to (30), estimating $\nabla_{\boldsymbol{\theta}}\mathcal{L}$ requires running $2C-1$ rotated primal PQCs, $2P$ times each for $\{\boldsymbol{\theta}^{+}_{p},\boldsymbol{\theta}^{-}_{p}\}_{p=1}^{P}$, while the dual PQC is set at $\boldsymbol{\phi}$.

Similarly, evaluating $\nabla_{\boldsymbol{\phi}}\mathcal{L}$ in (14c) requires running $2C-1$ rotated primal circuits set at $\boldsymbol{\theta}$, and the dual PQC for $2Q$ values for $\{\boldsymbol{\phi}^{+}_{q},\boldsymbol{\phi}^{-}_{q}\}_{q=1}^{Q}$ as

[TABLE]

where $\boldsymbol{\phi}^{\pm}_{q}\coloneqq\boldsymbol{\phi}\pm\frac{\pi}{2}\mathbf{e}_{q}$.

Altogether, each PD iteration requires running $\left(2C-1\right)\left(2P+1\right)$ rotated primal circuits and $2Q+1$ dual circuits. Each circuit is run for $S$ samples. According to [52], the number of PQC parameters $(P,Q)$ is expected to be $\tilde{\mathcal{O}}(\log N)$. If the number of colors grows approximately as $C=\mathcal{O}((\log N)^{3})$, each PD iteration in (14) incurs a computational complexity of $\tilde{\mathcal{O}}(\log N)$, which is more efficient compared to $\mathcal{O}(sN^{2})$ of the PD iteration in (3) running on a classical computer.

VI Convergence analysis

This section studies the convergence of the EG iterates in (22). For a general optimization problem, if the Lagrangian function is convex in the primal variables, concave in the dual variables, and has a Lipschitz gradient in primal and dual variables, the EG iterates converge to a saddle point of the Lagrangian function; see [34]. For nonconvex/nonconcave saddle-point problems, however, establishing convergence of the EG iterates is challenging. Reference [16] has recently shown that, under certain conditions, the EG iterates converge to an approximately stationary point of the Lagrangian function. Of course, a stationary point may not necessarily be a saddle point of the Lagrangian. Similar convergence claims have also been established in the stochastic setting, where EG updates rely on noisy yet unbiased gradient estimates in lieu of the actual Lagrangian gradients. Given that the Lagrangian function in (13) is nonconvex in $(\boldsymbol{\theta},\alpha)$ and nonconcave in $(\boldsymbol{\phi},\beta)$, we adopt the convergence proof from [16] for the EG iterates in (22) and provide an upper bound on the number of measurement samples across all EG iterations. To this end, the following lemma shows that the operator $\mathbf{g}(\mathbf{z})$ is $L$-Lipschitz continuous; see Appendix A.3 for a proof.

Lemma 3**.**

The vector-valued function $\mathbf{g}(\mathbf{z})$ defined in (21) is Lipschitz continuous with the Lipschitz constant

[TABLE]

where $\bar{\alpha}$ and $\bar{\beta}$ are upper bounds on the optimal values of $\alpha$ and $\beta$, respectively.

The convergence analysis relies on the ensuing assumption [16].

Assumption 1**.**

There exists a stationary point $\mathbf{z}^{*}$ of $\mathcal{L}(\mathbf{z})$ such that:

[TABLE]

for all $\mathbf{z}\in\mathbb{R}^{P+Q+2}$ and some parameter $\rho\in\left[0,\frac{1}{4\sqrt{2}L}\right)$, where $L$ is defined in (3).

Another challenge for our problem is that the EG updates in (22) rely on noisy estimates of $\mathbf{g}(\mathbf{z})$ because observables are measured only based on samples. Specifically, as discussed earlier, each PD iteration in (14) utilizes $S$ measurement samples per rotated primal/dual circuit to evaluate the Lagrangian gradients. Similarly, each step (22a) and (22b) of the EG updates requires $S$ measurement samples per circuit to estimate $\mathbf{g}(\mathbf{z})$ as:

[TABLE]

where $\hat{\mathbf{g}}_{s}(\mathbf{z})$ denotes the estimate of $\mathbf{g}(\mathbf{z})$ based on the single sample $s$. Under this stochastic setting, we adapt the result from [16] to bound the number of iterations $T$ and measurement samples $S$ such that

[TABLE]

The expectation in (35) is taken over the randomness of $\hat{\mathbf{g}}(\mathbf{z})$. Suppose now that the final output $\hat{\mathbf{z}}$ of the EG method is drawn by uniformly sampling at random from the EG iterates $\{\bar{\mathbf{z}}^{t}\}_{t=0}^{T-1}$. Then, it is easy to see that $\hat{\mathbf{z}}$ satisfies

[TABLE]

where the expectation is now taken over the noisy gradients as well as the sampling process to draw the final output $\hat{\mathbf{z}}$. The convergence result is formalized in the following theorem from [16].

Theorem 1** ([16]).**

Let $\mathbf{g}(\mathbf{z})$ be an $L$-Lipschitz operator that satisfies Assumption 1. Let $\hat{\mathbf{g}}(\mathbf{z})$ be an unbiased estimator of $\mathbf{g}(\mathbf{z})$ such that the variance of $\hat{\mathbf{g}}(\mathbf{z})$ is bounded as

[TABLE]

where $S$ is the number of measurement samples per rotated primal/dual circuit. Given an arbitrary initial point $\mathbf{z}^{0}$, run the EG iterations in (22) for $T$ times with the step size $\mu^{t}_{z}=\frac{1}{2\sqrt{2}L}$. Accordingly, let $\{\mathbf{z}^{t}\}_{t=1}^{T}$ and $\{\bar{\mathbf{z}}^{t}\}_{t=0}^{T-1}$ denote the sequences of points generated by the EG iteration. Select a point $\hat{\mathbf{z}}$ uniformly at random from $\{\bar{\mathbf{z}}^{t}\}_{t=0}^{T-1}$. If the number of iterations is selected as

[TABLE]

and each step (22a) and (22b) uses

[TABLE]

samples per iteration, then $\mathbb{E}\!\left[\left\|\hat{\mathbf{g}}(\hat{\mathbf{z}})\right\|_{2}\right]\leq\epsilon$.

We next expound on whether the Lagrangian function $\mathcal{L}(\mathbf{z})$ satisfies the requirements of Theorem 1. First, we were unable to verify Assumption 1, so it is taken as given. Second, the $L$-Lipschitz continuity of the operator $\mathbf{g}(\mathbf{z})$ has been shown in Lemma 3. Third, the estimator $\hat{\mathbf{g}}(\mathbf{z})$ is unbiased because Lagrangian gradients $\mathcal{L}(\mathbf{z})$ are computed via the PSR, and estimates of quantum observables based on sample averages are known to be unbiased estimators [26]. Lastly, the bound on the variance in (37) is established in the ensuing lemma shown in Appendix A.4.

Lemma 4**.**

The parameter $\sigma^{2}$ upper bounding the variance of the gradient estimator in (37) can be found as:

[TABLE]

Having determined the Lipschitz constant $L$ for $\mathbf{g}(\mathbf{z})$ and the variance bound $\sigma^{2}$ for $\hat{\mathbf{g}}(\mathbf{z})$, the following proposition provides an upper bound on the total number of measurement samples running on all circuits.

Proposition 1**.**

Given an arbitrary point $\mathbf{z}^{0}$, run the EG iteration in (22) for $T$ times with the step size $\mu^{t}_{z}=\frac{1}{2\sqrt{2}L}$. Accordingly, let $\{\mathbf{z}^{t}\}_{t=1}^{T}$ and $\{\bar{\mathbf{z}}^{t}\}_{t=0}^{T-1}$ denote the sequences of points generated by the EG iteration in (22). Choose a point $\hat{\mathbf{z}}$ uniformly from $\{\bar{\mathbf{z}}^{t}\}_{t=0}^{T-1}$. If the operator $\mathbf{g}(\mathbf{z})$ satisfies Assumption 1, then after

[TABLE]

measurement samples across all iterations and for all circuits, it holds that

[TABLE]

The proposition can be established by multiplying the number of primal/dual circuits $((2P+1)(2C-1)+2Q+1)$, the number of iterations in (38), and twice the number of samples in (39). Here, we upper bound the factor $(8+\sqrt{2}L\rho)$ appearing in the numerator of (39) by $8.25$ for $\rho\in\left[0,\frac{1}{4\sqrt{2}L}\right)$.

Proposition 1 establishes that the total number of measurement samples scales polynomially with the number of colors $C$, the number of PQC parameters $(P,Q)$, the Lipschitz constant $L$, and the variance bound $\sigma^{2}$. Numerical tests with the OPF in Section V.2 suggest that the number of colors scales as $C=\mathcal{O}((\log N)^{3})$, a scaling we expect to extend to other problems with similar graph sparsity. Additionally, the parameter counts $P$ and $Q$ typically scale as $\tilde{\mathcal{O}}(\log N)$ [52]. The constants $L$ and $\sigma^{2}$ depend on problem-specific bounds $\bar{\alpha}$ and $\bar{\beta}$. If these bounds scale as $\tilde{\mathcal{O}}(\log N)$, then the total number of measurement samples, and hence the overall complexity, becomes $\tilde{\mathcal{O}}(\log N)$, offering a potential advantage over the PD method running on a classical computer. The latter requires $\mathcal{O}(sN^{2})$ operations multiplying with the number of PD iterations in (3). The following section models the OPF as a QCQP and expounds how the sparsity structure of the power network aligns with the measurement protocol presented in Section V.

VII Optimal power flow as a QCQP

The OPF is a large-scale, nonconvex optimization problem. Given forecasts of electric load demands for the next control period, the goal of the OPF is to schedule generators and flexible demand most economically while complying with engineering constraints imposed by the power transmission system and the generation units. This section states the OPF and poses it as a QCQP.

Let us first review a power system model; see [32] and references therein. An electric power system can be represented by an undirected connected graph $\mathcal{G}\coloneqq\{\mathcal{N},\mathcal{E}\}$. The vertex set $\mathcal{N}\coloneqq\{1,\ldots,N\}$ comprises $N$ nodes. Nodes correspond to points where electric power is produced by generators or consumed by loads. The edge set $\mathcal{E}\coloneqq\{\ell=(n,m):n,m\in\mathcal{N}\}$ consists of $L_{e}\coloneqq|\mathcal{E}|$ transmission lines, each connecting two nodes. Because each node is connected only to a few other nodes, the number of lines, $L_{e}$, is typically a small multiple of $N$.

The power system can be described by a vector $\mathbf{v}\in\mathbb{C}^{N}$ of the complex AC voltages (voltage phasors) experienced at all nodes. Because many other quantities of interest can be expressed as functions of nodal voltages, the vector $\mathbf{v}$ is typically selected as the state of the power system. For example, the vector of AC currents, $\mathbf{i}$, injected into all nodes can be expressed as $\mathbf{i}=\mathbf{Y}\mathbf{v}$, where $\mathbf{Y}\in\mathbb{C}^{N\times N}$ is the node admittance matrix, which can be thought of as a complex-valued Laplacian matrix of the power system graph $\mathcal{G}$. Other grid quantities, such as the nodal power injections, line power flows, squared voltage magnitudes, and squared line magnitudes, can all be expressed as quadratic functions of $\mathbf{v}$, as $\mathbf{v}^{\dagger}\mathbf{M}_{m}\mathbf{v}$ for some Hermitian matrix $\mathbf{M}_{m}$. If $v_{n}$ and $i_{n}$ are the AC voltage and current at node $n$, respectively, the complex power injected into the power network through node $n$ is defined as $s_{n}=v_{n}i_{n}^{*}=p_{n}+\iota q_{n}$, whose real part $p_{n}$ is termed active power, and imaginary part $q_{n}$ is reactive power.

The cost function of the OPF is usually a linear function of some of the active power injections. To comply with generation capacities and load demands, the OPF imposes upper and lower limits on active and reactive power injections in the form of

[TABLE]

Nodal voltage magnitudes are constrained within upper and lower limits as

[TABLE]

Line current magnitudes are upper-limited by line capacities as

[TABLE]

Since all quantities mentioned above are quadratic in $\mathbf{v}$, the OPF can be posed as a quadratically constrained quadratic program (QCQP) over $\mathbf{v}$ as [38, 32]: where $b_{m}$, for all $m$, is a given real-valued parameter determined by lower/upper limits and each $\mathbf{M}_{m}$ is a known Hermitian matrix. For example, if constraint \underaccent{\bar}{p}_{n}\leq p_{n} is indexed as the $m$-th constraint, then b_{m}=-\underaccent{\bar}{p}_{n}. The OPF involves quadratic equality constraints, e.g., when \underaccent{\bar}{p}_{n}=\bar{p}_{n} for a subset of nodes $n\in\mathcal{N}$. Each of these equality constraints corresponds to two inequality constraints in (QCQP). Overall, the number of OPF constraints can be several times the network size, such as $M\simeq 8N$.

The precise form of $\mathbf{M}_{0}$ and $\{(\mathbf{M}_{m},b_{m})\}_{m=1}^{M}$ is delineated in Appendix A.1. It is worth stressing that every matrix $\mathbf{M}_{m}$ is highly sparse. This is because every node of a power system is directly connected to only a few other nodes. The number of nonzero entries per matrix $\mathbf{M}_{m}$ is upper bounded by a small constant much smaller than $N$. The $(i,j)$-th off-diagonal entry across every matrix $\mathbf{M}_{m}$ is nonzero only if nodes $i$ and $j$ are connected through a transmission line, that is, if $(i,j)\in\mathcal{E}$. Consequently, the sparsity pattern (location of nonzero entries) for each $\mathbf{M}_{m}$ is determined by a small subset of $\mathcal{E}$. Therefore, the union of the sparsity patterns across all matrices $\mathbf{M}_{m}$ coincides with the sparsity pattern of the nodal admittance matrix $\mathbf{Y}$. The sparsity of the matrices ${\mathbf{M}_{m}}$ combined with node permutation allows the number of colors $C$ to scale favorably as $\tilde{\mathcal{O}}((\log N)^{3})$, as demonstrated in Section V.3. This efficient coloring enables the measurement protocol of Section V.

Nonetheless, the overall measurement complexity derived in Proposition 1 also depends on the Lipschitz constant $L$ and the variance bound $\sigma^{2}$ whose scaling is unfavorable for OPF. This is because the magnitudes of the entries in $\mathbf{v}$ in (QCQP) are constrained to lie in the range of $[0.9,1.1]$, so that $\bar{\alpha}=1.1\sqrt{N}$. Furthermore, as OPF entails equality constraints for nodes connected to loads, each of these equality constraints is converted to two inequality constraints upon (QCQP). Since these inequalities should be active at optimality, the complementary slackness of the Karush–Kuhn–Tucker conditions written for (QCQP) implies that the associated dual variables are generally nonzero [8, Ch. 5]. Because the number of load nodes in real-world power systems is typically proportional to $N$, it follows that $\bar{\beta}$ may scale with $\sqrt{N}$ too. Overall, in practice, parameters $L$ in (3) and $\sigma^{2}$ in (4) could grow as $N$ and $N^{2}$, respectively, which indicates the total number of measurement samples in (41) could scale with $N^{4}$. This large number of measurement samples ultimately presents a challenge to attaining quantum advantage. A similar undesirable dependence of the measurement sample complexity on the scaling prefactors has also been encountered in [56], where the objective function entails expectations multiplied by prefactors exponentially in $\log N$, like the one in (10).

VIII Numerical Tests

The performance of the proposed doubly variational quantum QCQP was validated for the OPF problem (hereafter referred to as Q-OPF) using the IEEE 57-node power system benchmark [2]. For this benchmark, we have $N=57$ primal variables and $M=422$ dual variables. Therefore, the number of qubits needed to capture $\mathbf{v}(\boldsymbol{\theta})$ and $\boldsymbol{\lambda}(\boldsymbol{\phi})$ is 6 and 9, respectively. To simplify the implementation, nodes hosting generators were assumed to have no load, i.e., $p^{d}_{n}=q^{d}_{n}=0$ for all $n\in\mathcal{N}_{g}$. Reactive power demands were set so that load nodes have $q^{d}_{n}=0.33\times p^{d}_{n}$ for all $n\in\mathcal{N}_{l}$. We generated 15 different instances of the OPF by scaling original load demands $\{(p^{d}_{n},q^{d}_{n})\}_{n\in\mathcal{N}_{\ell}}$ by factors drawn uniformly distributed within $[0.90,1.05]$, independently per node and problem instance. Each of the 15 OPF instances was solved to global optimality using MATPOWER, a power system toolbox in MATLAB [65].

Selecting PQC architecture. To select the architecture of PQCp, we numerically tested different architectures. Each candidate architecture was validated over all 15 problem instances. For each problem instance $k$, we found the OPF solution $\mathbf{v}_{k}^{*}$ using MATPOWER, and then solved the variational problem

[TABLE]

with $\ket{\boldsymbol{\psi}(\boldsymbol{\theta}_{k})}$ generated by the specified PQCp architecture. Assuming the previous cost is optimized to global optimality, it measures the alignment between the generated quantum state and the actual OPF solution. The cost is inspired by [9] and is lower-bounded by zero. We tested architectures with different gate types and number of layers $L$, as detailed in Appendix A.5. We selected the PQC that attained the smallest average cost across all 15 OPF problem instances. A similar procedure was followed to determine the PQCd architecture. The optimal PQCp architecture consists of 10 layers, with each layer comprising the sequence $\mathbf{R}_{y}(\theta)-\text{CX}-\mathbf{R}_{z}(\theta^{\prime})-\text{CX}$ gates, resulting in $P=2\times 10\times 6=120$ parameters. The optimal PQCd architecture is formed by 35 layers of $\mathbf{R}_{y}(\phi)-\text{CX}$ gates, yielding $Q=1\times 35\times 9=315$ parameters. The average fitting errors obtained were $2\times 10^{-4}$ and $2\times 10^{-5}$ for PQCp and PQCd, respectively. This demonstrates that the selected PQC pairs are sufficiently expressive.

Algorithm details. PQC parameters $\boldsymbol{\theta}$ and $\boldsymbol{\phi}$ were initialized uniformly at random within the range $[0,2\pi]$. Because voltages are normally around 1 per unit, the scaling variable $\alpha$ was initialized to $\sqrt{N}$; the variable $\beta$ was initialized to $2\left|\mathcal{N}_{\ell}\right|$. Step sizes were set according to the exponentially decaying rule as $\mu^{t}_{\theta}=0.015\times 0.99985^{t}$, $\mu^{t}_{\phi}=0.01\times 0.99985^{t}$, and $\mu^{t}_{\alpha}=\mu^{t}_{\beta}=10^{-5}\times 0.999^{t}$. Both PD and EG iterations were terminated when both conditions $\left\|\boldsymbol{\theta}^{t}-\boldsymbol{\theta}^{t-1}\right\|_{2}\leq 10^{-6}$ and $\left\|\boldsymbol{\phi}^{t}-\boldsymbol{\phi}^{t-1}\right\|_{2}\leq 10^{-6}$ were satisfied. Simulation scripts were written in Python and run on Pennylane’s exact simulator [5].

We first explored the performance of Q-OPF in finding AC OPF solutions. Rather than evaluating errors on OPF solutions in terms of voltage phasors $\mathbf{v}$, we focused on finding the active power injection and voltage magnitude at all generator nodes. This is because power system operators are primarily interested in generator setpoints. Define the vector of generator setpoints $\mathbf{x}_{g}=[\mathbf{p}_{g};\mathbf{v}_{g}]^{\top}$ with $\mathbf{p}_{g}$ and $\mathbf{v}_{g}$ collecting generator active power injections and voltage magnitudes, respectively. Generator active power injections can be computed from (45b) given $\mathbf{v}$ and the corresponding load demand $p_{n}^{d}$. In terms of dual variables, we focus on Lagrange multipliers $\lambda_{m}$ associated with power balance and line limit constraints because these are used in electricity markets to compute the so-called locational marginal prices.

We compared solving (QCQP) and (QCQPθ) in terms of the relative errors $\left\|\mathbf{x}_{g}-\mathbf{x}_{g}^{*}\right\|_{2}/\left\|\mathbf{x}_{g}^{*}\right\|_{2}$ and $\left\|\boldsymbol{\lambda}-\boldsymbol{\lambda}^{*}\right\|_{2}/\left\|\boldsymbol{\lambda}^{*}\right\|_{2}$ of the found generator setpoints $\mathbf{x}$ and dual variables $\boldsymbol{\lambda}$ from their global optimal solutions found using MATPOWER. The entries of the dual vector $\boldsymbol{\lambda}$ correspond to the Lagrangian multipliers associated with 280 power balance and line limit constraints. The non-convex (QCQP) was solved classically using the PD iteration in (3) and the EG iteration in (22) over voltage phasors $\mathbf{v}$ and dual variables $\boldsymbol{\lambda}$. Voltages were initialized at the flat voltage profile as $\mathbf{v}^{0}=\mathbf{1}$, while the entries of $\boldsymbol{\lambda}$ were initialized independently at random based on the standard normal distribution and scaled by $2\left|\mathcal{N}_{\ell}\right|$. The related step sizes were set as $\mu_{v}=\mu_{\lambda}=10^{-3}\times 0.9999^{t}$, and PD iterates were terminated when $\left\|\mathbf{v}^{t}-\mathbf{v}^{t-1}\right\|_{2}\leq 10^{-6}$ and $\left\|\boldsymbol{\lambda}^{t}-\boldsymbol{\lambda}^{t-1}\right\|_{2}\leq 10^{-6}$ were met.

Table 1 shows the relative errors for generator setpoints and dual variables attained by solving (QCQPθ) and (QCQP) via the PD and EG iterations. Among the four models, solving (QCQPθ) using the EG iteration achieved the smallest relative errors on $\mathbf{x}_{g}$ and $\boldsymbol{\lambda}$. Comparing algorithms when solving either (QCQPθ) or (QCQP), the EG method consistently outperforms the PD method. Similarly, fixing the algorithm to be either EG or PD, solving (QCQPθ) using a hybrid quantum-classical computer offers smaller errors compared to solving (QCQP) using the PD method.

Given that the EG iteration outperforms the PD iteration, Fig. 6 demonstrates the dual variables obtained by solving (QCQPθ) and (QCQP) by the EG method. The Q-OPF approach can approximate dual variables reasonably well. This test corroborates that Q-OPF is capable of finding near-optimal AC OPF solutions.

Feasibility. Does the power system safely operate under the generator setpoints computed by Q-OPFs? To answer this question, for each instance, given the input load demands and the setpoints $\mathbf{x}_{g}$ obtained by Q-OPF, we found an AC power flow solution via MATPOWER. The obtained voltages were used to evaluate all 222 inequality constraints in (45). Violations in line current and generator constraints were normalized by their maximum limits. We employed three constraint violation metrics: a) the average number of constraint violations surpassing a normalized magnitude of $10^{-6}$ across OPF instances; b) the maximum constraint violation incurred across instances; c) the violations averaged over all constraints and instances. As shown in Table 1, solving (QCQPθ) using the EG iteration attains the best constraint violation metrics among the tested models.

Figure 7 shows the relative errors in terms of the Lagrangian function $\left|\mathcal{L}-P^{*}_{\theta}\right|/P^{*}_{\theta}$ for the last EG iteration. The relative errors of the Q-OPF are consistently below 1.5% and smaller than those of QCQP in 14 out of 15 instances. This test indicates that under the Q-OPF’s generator setpoints, all physical constraints of the network are nearly satisfied.

IX Conclusions

This work proposes a doubly variational quantum approach to solve QCQP/SDPs, where the observables are induced by underlying sparse graphs. Primal/dual variables are modeled by scaling the state and the PMF corresponding to two PQCs. These two PQCs are trained using exclusively the gradient estimates of the Lagrangian function. The training procedure aims to seek an approximate stationary point of the Lagrangian function. By adopting the XBM method and leveraging graph sparsity, the related observables are judiciously permuted into a banded form, enabling efficient evaluation of the Lagrangian function. The proposed algorithmic framework is applied to the OPF problem. Numerical tests on the IEEE 57-node benchmark power system using PennyLane’s quantum simulator have demonstrated that: i) the EG iteration performs markedly better than the PD iteration; ii) solving (QCQPθ) on PQCs offers smaller errors in terms of the generator setpoint and dual variables, as well as smaller constraint violations compared to solving (QCQP) on a classical computer using PD; and iii) the Q-OPF found generator setpoints and dual variables with small errors, and under the generator setpoints obtained by the Q-OPF, the network constraints are nearly satisfied.

The proposed doubly variational quantum paradigm sets a solid foundation for exploring exciting directions: d1) Thus far, the PQC parameters $(\boldsymbol{\theta},\boldsymbol{\phi})$ have been optimized to solve one problem instance, and parameters $b_{m}$ have solely appeared on the diagonal of the observable $\mathbf{S}$. Nevertheless, vector $\mathbf{b}$ can be encoded as parameters of another unitary $\mathbf{T}(\mathbf{b})$, under which the corresponding primal and dual states are $\ket{\boldsymbol{\psi}(\boldsymbol{\theta},\mathbf{b})}=\mathbf{T}(\mathbf{b})\mathbf{U}(\boldsymbol{\theta})\ket{\mathbf{0}}$ and $\ket{\boldsymbol{\xi}(\boldsymbol{\phi},\mathbf{b})}=\mathbf{T}(\mathbf{b})\mathbf{V}(\boldsymbol{\phi})\ket{\mathbf{0}}$. Measuring the states $\ket{\boldsymbol{\psi}(\boldsymbol{\theta},\mathbf{b})}$ and $\ket{\boldsymbol{\xi}(\boldsymbol{\phi},\mathbf{b})}$ and utilizing the proposed doubly variational algorithm for the Lagrangian in (10) evaluated across instances can solve multiple problem instances; d2) Prudently construct an alternative objective function for (QCQP) to escape the problem of measurement samples scaling with $\bar{\alpha}$ and $\bar{\beta}$; d3) Designing network-informed PQCs and/or testing the Q-OPF using dynamic PQCs and quantum recurrent embedding neural networks will be imperative for dealing with large-scale networks; and d4) The Q-OPF can be a building block to investigate other computational tasks, such as multi-period OPF, which entails the time-coupled ramp constraints on generators.

Acknowledgements.

VK and TVL acknowledge support from the US National Science Foundation under grant no. 2412947, and the US Office of Naval Research under grant N000142412614. MMW acknowledges support from the National Science Foundation under grant no. 2329662 and from the Cornell School of Electrical and Computer Engineering.

Appendix A Appendix

A.1 Power system modeling

A transmission line connecting node $n$ to node $m$ is represented by its series conductance $G_{nm}$ and series susceptance $B_{nm}$. Given $G_{nm}$ and $B_{nm}$, power injections at node $n$ are computed through nodal voltages based upon the following quadratic power flow equations:

[TABLE]

[TABLE]

where $(v_{n}^{r},v_{n}^{i})$ and $(v_{m}^{r},v_{m}^{i})$ are the real and imaginary components of the voltages at nodes $n$ and $m$, respectively. Define the $N\times N$ node admittance matrix as $\mathbf{Y}\coloneqq\mathbf{G}+\iota\mathbf{B}$, where $\mathbf{G}$ collects $G_{nm}$ and $\mathbf{B}$ collects $B_{nm}$ of the network lines. Then the power flow equations can be compactly represented as

[TABLE]

where $\mathbf{M}_{p_{n}}$ and $\mathbf{M}_{q_{n}}$ are $N\times N$ Hermitian matrices taken the forms as

[TABLE]

[TABLE]

Similarly, the squared voltage magnitude at node $n$ is represented as

[TABLE]

where $\mathbf{M}_{v_{n}}\coloneqq\mathbf{e}_{n}\mathbf{e}_{n}^{\top}$. Given a line $(n,m)\in\mathcal{E}$ with the series admittance $Y_{nm}$, the complex current flowing on from node $n$ to node $m$ is provided as $i_{nm}=Y_{nm}(v_{n}-v_{m})$. This poses $|i_{nm}|=\mathbf{v}^{\dagger}\mathbf{M}_{i_{nm}}\mathbf{v}$ where $\mathbf{M}_{i_{nm}}\coloneqq\left|Y_{nm}\right|(\mathbf{e}_{n}-\mathbf{e}_{m})(\mathbf{e}_{n}-\mathbf{e}_{m})^{\top}$.

The power injection $(p_{n},q_{n})$ consists of a dispatchable constituent $(p^{g}_{n},q^{g}_{n})$ and an inflexible constituent $(p^{d}_{n},q^{d}_{n})$ such that $p_{n}=p^{g}_{n}-p^{d}_{n}$ and $q_{n}=q^{g}_{n}-q^{d}_{n}$. The former represents the power dispatch of a generator or a flexible load placed at node $n$, and the latter expresses the inelastic load supplied at node $n$. Let nodes with dispatchable power injections belong to a subset $\mathcal{N}_{g}\subseteq\mathcal{N}$. The remaining nodes hosted inelastic loads form the subset $\mathcal{N}_{\ell}=\mathcal{N}\backslash\mathcal{N}_{g}$. Node $n=1$ is referred to as the angle reference node, under which it follows $\mathbf{v}^{\dagger}\mathbf{e}_{1}\mathbf{e}_{1}^{\top}\mathbf{v}=1$, where $\mathbf{e}_{1}$ is the first column of the identity. To simplify the exposition, each node in $\mathcal{N}_{g}$ is assumed to host only one dispatchable unit. Given the inflexible loads $\{p^{d}_{n},q^{d}_{n}\}_{n\in\mathcal{N}}$, the OPF’s objective is to minimize the dispatch cost of generators and flexible load while respecting the resource and network limits. The OPF is formulated as the ensuing QCQP over nodal voltages and power of generators:

[TABLE]

Constraints (45b)–(45e) capture the power flow equations at generator and load nodes. Constraints (45f)–(45g) confine limits of generators and flexible loads. Constraint (45h) limits the squared voltage magnitudes within the specified range, and constraint (45i) represents the squared magnitude of the angle reference node. Constraint (45j) enforces the line current magnitude within the line ratings. By substituting $\{p^{g}_{n},q^{g}_{n}\}$ in (45b)–(45c) into the objective of (45), the OPF model is represented as a QCQP in terms of the nodal voltage vector $\mathbf{v}$.

A.2 Proof of Lemma 1

The claim follows readily by substituting $\mathbf{M}$ into the observable as

[TABLE]

A.3 Proof of Lemma 3

Quantum observables are trigonometric functions of $(\boldsymbol{\theta},\boldsymbol{\phi})$. Moreover, the Lagrangian function $\mathcal{L}(\mathbf{z})$ is quadratic in $(\alpha,\beta)$. Therefore, the operator $\mathbf{g}(\mathbf{z})$ is continuously differentiable with respect to $\mathbf{z}$. Let $\mathbf{J}(\mathbf{z})\coloneqq\nabla_{\mathbf{z}}\mathbf{g}(\mathbf{z})$ denote the Jacobian matrix of $\mathbf{g}(\mathbf{z})$. If the spectral norm of $\mathbf{J}(\mathbf{z})$ is upper bounded by a constant $L$, then $\mathbf{g}(\mathbf{z})$ is $L$-Lipschitz continuous [22, Lemma 2.6]. Based on matrix norm inequalities, the spectral norm of $\mathbf{J}(\mathbf{z})$ can be bounded as:

[TABLE]

where $\left\|\mathbf{J}(\mathbf{z})\right\|_{1}$ is the maximum absolute column sum of $\mathbf{J}(\mathbf{z})$, and $\left\|\mathbf{J}(\mathbf{z})\right\|_{\infty}$ is the maximum absolute row sum. Because $\mathcal{L}(\mathbf{z})$ has continuous second-order partial derivatives, Schwarz’s theorem predicates that $\frac{\partial^{2}\mathcal{L}}{{\partial z_{i}\partial z_{j}}}=\frac{\partial^{2}\mathcal{L}}{{\partial z_{j}\partial z_{i}}}$, so that $\left\|\mathbf{J}(\mathbf{z})\right\|_{1}=\left\|\mathbf{J}(\mathbf{z})\right\|_{\infty}$ and (46) simplifies as

[TABLE]

where $\mathbf{e}_{j}$ is the $j$-th column of the identity.

Consider first the first $P$ rows of $\mathbf{J}(\mathbf{z})$ corresponding to the partial Hessian matrix $\nabla_{\boldsymbol{\theta},\mathbf{z}}^{2}\mathcal{L}(\mathbf{z})$. More specifically, the absolute row sum for row $j\in\{1,\ldots,P\}$ is

[TABLE]

Applying the triangle inequality on each summand of the first term in (48) yields

[TABLE]

The second-order partial derivatives of quantum expectations can be computed according to the PSR as [44]:

[TABLE]

for all $m$, $p$, and $j$.

Based on the latter bounds, we get that:

[TABLE]

The second term of (48) can be upper bounded similarly using the PSR for first-order partial derivatives:

[TABLE]

The summands in the third term of (48) can be upper bounded as:

[TABLE]

for all $q$. Then, the third term of (48) is bounded as

[TABLE]

The fourth term in (48) can be bounded as

[TABLE]

Combining (A.3)–(A.3) yields

[TABLE]

For the remaining rows of $\mathbf{J}(\mathbf{z})$ corresponding to the partial Hessian matrices $\nabla_{\alpha,\mathbf{z}}^{2}\mathcal{L}(\mathbf{z})$, $-\nabla_{\boldsymbol{\phi},\mathbf{z}}^{2}\mathcal{L}(\mathbf{z})$, and $-\nabla_{\beta,\mathbf{z}}^{2}\mathcal{L}(\mathbf{z})$, the norms $\left\|\mathbf{e}_{j}^{\top}\mathbf{J}(\mathbf{z})\right\|_{1}$ can be upper bounded using similar arguments. A key point here is that $\left\|\mathbf{M}\right\|=\max_{m}\left\|\mathbf{M}_{m}\right\|$ because $\mathbf{M}$ in (18) is block diagonal. To avoid repetition and tedious algebra, we omit the detailed derivations and present the final bounds.

[TABLE]

[TABLE]

[TABLE]

According to (47), the Lipschitz constant $L$ of $\mathbf{g}(\mathbf{z})$ is the maximum of the right-hand sides of (A.3)–(A.3). If $\alpha\geq 4$, the right-hand side (RHS) of (A.3) is larger than the RHS of (A.3). The condition $\alpha\geq 4$ holds trivially if the problem size is $N>32$. Moreover, if $\beta\geq 2$, the RHS of (A.3) is larger than that of (A.3). The condition $\beta\geq 2$ is expected to occur for problems with many constraints. Therefore, we are left with the RHS of (A.3) and (A.3), whose maximum can be compactly expressed as

[TABLE]

Then, the claim of Lemma 3 follows readily.

A.4 Proof of Lemma 4

We first upper bound the variance of the sample-average estimators of the three terms in (10). Note that the three terms in (10) can be decomposed into quantum expectations across $c$. This is obvious for $F_{0}(\boldsymbol{\theta})$ and $G(\boldsymbol{\phi})$. Regarding $F(\boldsymbol{\theta},\boldsymbol{\phi})$, we obtain from (V.3):

[TABLE]

Because the three expectations in (10) are defined over different states and observables, we consider a general expectation $E\coloneqq\braket{\boldsymbol{\zeta}|\mathbf{N}|\boldsymbol{\zeta}}$, where $\ket{\boldsymbol{\zeta}}$ and $\mathbf{N}$ are dimensionally compatible and $\mathbf{N}$ is a Hermitian matrix with the same color decomposition as matrices $\mathbf{M}_{m}$’s. The expectation $E$ can be color-decomposed as

[TABLE]

where $\mathbf{N}^{c}=\mathbf{U}_{c}\mathbf{\Lambda}^{c}\mathbf{U}_{c}^{{\dagger}}$ is the eigendecomposition of $\mathbf{N}^{c}$ and $\ket{\boldsymbol{\zeta}_{c}}=\mathbf{U}_{c}^{{\dagger}}\ket{\boldsymbol{\zeta}}$.

Upon measuring $\ket{\boldsymbol{\zeta}_{c}}$ in the computational basis using $S$ measurement samples, the binary outcome $\ket{i}$ is observed with probability $|\zeta^{i}_{c}|^{2}\coloneqq|\braket{i|\boldsymbol{\zeta}_{c}}|^{2}$. For each sample $s$, define a random variable $\hat{E}_{s}^{c}$, taking the value $\braket{i|\mathbf{\Lambda}^{c}|i}$ when outcome $\ket{i}$ is observed while sampling $\ket{\boldsymbol{\zeta}_{c}}$. Subsequently, the sample-average estimator of $E$ can be expressed as

[TABLE]

The variance of $\hat{E}$ is bounded in the next lemma.

Lemma 5**.**

The variance of the estimator in (58) can be computed as:

[TABLE]

and can be upper bounded as

[TABLE]

Proof.

The variance of each $\hat{E}_{s}^{c}$ can be computed and upper bounded as follows:

[TABLE]

Because the random variables $\hat{E}_{s}^{c}$ are independent across $c$ and $s$, we get that

[TABLE]

for any $s$. Substituting $\operatorname{Var}(\hat{E}_{s}^{c})$ above yields the expression in (59). The final bound in (60) follows readily. ∎

We are now ready to upper bound the variance of the estimate $\hat{\mathbf{g}}(\mathbf{z})$ per Theorem 1.

Proof of Lemma 4.

The RHS of (37) can be expanded as

[TABLE]

Among four terms in (61), the terms $\operatorname{Var}(\hat{\nabla}_{\alpha}\mathcal{L}(\mathbf{z}))$ and $\operatorname{Var}(\hat{\nabla}_{\beta}\mathcal{L}(\mathbf{z}))$ are simpler to bound. Let $\hat{F}_{0}(\boldsymbol{\theta})$ and $\hat{F}(\boldsymbol{\theta},\boldsymbol{\phi})$ denote the sample-based estimates of $F_{0}(\boldsymbol{\theta})$ and $F(\boldsymbol{\theta},\boldsymbol{\phi})$, respectively. Then, we can express the first term in (61) as:

[TABLE]

It follows from (60) that

[TABLE]

Regarding $\operatorname{Var}(\hat{F}(\boldsymbol{\theta},\boldsymbol{\phi}))$, because $\mathbf{M}$ is a block-diagonal matrix with blocks $\mathbf{M}_{m}$, it follows that

[TABLE]

Substituting (63)–(64) into (A.4) provides

[TABLE]

Let $\hat{\partial}_{\theta_{p}}F_{0}(\boldsymbol{\theta})$ and $\hat{\partial}_{\theta_{p}}F(\boldsymbol{\theta},\boldsymbol{\phi})$ denote the sample-based estimates of $\partial_{\theta_{p}}F_{0}(\boldsymbol{\theta})$ and ${\partial}_{\theta_{p}}F(\boldsymbol{\theta},\boldsymbol{\phi})$, respectively. The second expectation in (61) entails computing

[TABLE]

Analogously, the third and fourth expectations in (61) can be upper bounded as

[TABLE]

[TABLE]

Summing (65)–(68), the RHS of (37) is bounded as:

[TABLE]

∎

A.5 PQC architectures

The tested PQC architectures are reported in Table 2.