Variational Quantum Algorithms for Semidefinite Programming

A semidefinite program (SDP) is a particular kind of convex optimization problem with applications in operations research, combinatorial optimization, quantum information science, and beyond. In this work, we propose variational quantum algorithms for approximately solving SDPs. For one class of SDPs, we provide a rigorous analysis of their convergence to approximate locally optimal solutions, under the assumption that they are weakly constrained (i.e., $N\gg M$, where $N$ is the dimension of the input matrices and $M$ is the number of constraints). We also provide algorithms for a more general class of SDPs that requires fewer assumptions. Finally, we numerically simulate our quantum algorithms for applications such as MaxCut, and the results of these simulations provide evidence that convergence still occurs in noisy settings.


Introduction
Semidefinite programming (SDP) is one of the most important tools in optimization developed over the past few decades.One can see SDPs as a natural extension of the better known linear programs (LPs), in which the vector inequalities of LPs are replaced by matrix inequalities.One of the reasons underlying the importance of SDPs is their applicability to a broad range of problems, including approximation algorithms for combinatorial optimization [1], control theory [2], and sum-of-squares [3].Additionally, a variety of quantum information problems can be formulated as SDPs, including state discrimination [4,5], upper bounds on quantum channel capacity [6,7], and self-testing [8].
The power of SDPs lies with the fact that they can be solved efficiently in polynomial time using classical algorithms such as the celebrated interior-point method [9].Although SDPs can be solved efficiently using classical techniques, as the size of the input matrices increases, many first-order and second-order algorithms incur significant computational overhead due to the expensive gradient computation at each iteration.For this reason, it is imperative to design more efficient algorithms for solving SDPs.
Given the speed-ups of quantum algorithms over classical algorithms for a variety of problems [10,11,12,13,14], it is natural to ask if there exists a quantum algorithm that can solve SDPs efficiently.This question was positively answered in [15], wherein a quantum algorithm was proposed and proven to have a quadratic speedup over the classical Arora-Kale algorithm [16].Following this initial result, more efficient quantum algorithms for solving SDPs were later developed [17,18].
Although quantum algorithms have been theoretically proven to outperform known classical algorithms for many applications, a fault-tolerant quantum computer is required to reap their benefits.Currently, fault-tolerant quantum computers are not available, and we are instead in the Noisy Intermediate-Scale Quantum (NISQ) era [19].Google constructed a noisy quantum computer with just 54 qubits [20], but it remains an open challenge to design a fault-tolerant quantum computer that requires millions of qubits for successful operation.Some of the quantum information science research community is focused on determining the power of noisy quantum computers and is designing quantum algorithms that acknowledge the limitations of such devices [21,22,23,24], such as a finite number of gates and qubits, noisy gate execution, and rapid decoherence of qubits.
Variational quantum algorithms (VQAs) constitute an important class of NISQ-friendly algorithms [21,22].VQAs can be seen as hybrid quantum-classical algorithms that have a classical computer available for optimization, only calling a quantum subroutine for tasks that are not efficiently solvable by it.Hitherto, VQAs have been proposed for numerous computational tasks, for which there are also known quantum and classical algorithms [21].

Main Idea and Setup
In this paper, we consider the following three different kinds of SDPs, and we present variational quantum algorithms for all three of them: • General Form (GF): The general form of an SDP can be concisely written as follows: where C P S N , B P S M , and the map Φ is Hermiticity-preserving (Definition 3).
Here, the notation S N denotes the set of N ˆN Hermitian operators.For this case, we do not assume that the SDPs are weakly constrained; i.e., we do not assume that the dimension M of the constraint variable B is much smaller than the dimension N of the objective variable C. Additionally, to generalize it further, we consider an inequality-constrained problem.For solving these general SDPs, we propose a variational quantum algorithm (Algorithm 1).
• Standard Form (SF): Here, we consider the Hermiticity-preserving map Φ and the Hermitian operator B to have a diagonal form, as given in (15).Specifically, in this case, we set ΦpXq " diag pTrrA 1 Xs, . . ., TrrA M Xsq , where A 1 , . . ., A M P S N and b 1 , . . ., b M P R.This form of an SDP is well known in the convex optimization literature [28], and most exact or approximation algorithms designed for combinatorial optimization problems are based on this form.For this case, we assume that the SDPs are weakly constrained, i.e., N " M .We further categorize them based on the nature of the constraints as follows: -Equality Constrained Standard Form (ECSF): Here, we consider equality constraints, i.e., ΦpXq " B. For solving such SDPs, we propose a variational quantum algorithm (Algorithm 2) and establish its convergence rate and total iteration complexity.
-Inequality Constrained Standard Form (ICSF): Here, we consider inequality constraints, i.e., ΦpXq ď B, and we propose a variational quantum algorithm for this case (Algorithm 3).
In our methods, we first reduce these constrained optimization problems to their unconstrained forms by employing a series of identities.Second, we express the final unconstrained form as a function of expectation values of the input Hermitian operators, i.e., C, A 1 , . . ., A M P S N .
For solving these final unconstrained formulations using a gradient-based method, we need access to the full gradient of the objective function at each iteration.However, when the dimension N of the input Hermitian operators is large, the evaluation of the gradient of the objective function using a classical computer becomes computationally expensive.Therefore, we delegate this gradient computation to parameterized quantum circuits, and we use a technique called the parameter-shift rule [29,30,31] to evaluate the partial derivatives of the objective function with respect to circuit parameters.Please refer to Section 2.2 for more details on the parameter-shift rule.
We design variational quantum algorithms for solving these unconstrained optimization problems.Our methods provide bounds on the optimal values, due to the reduction in the search space, as well as the non-convex nature of the objective function landscape in terms of quantum circuit parameters.We do not assume that the final objective function is convex with respect to these parameters, as generally it is non-convex [32].In general, finding a globally optimal point for a non-convex function is known to be NP-hard [33, Section 2.1], and so an important question regards the type of solutions that we can guarantee in such a scenario.In the classical optimization literature, the notion of approximate stationary points (Definition 10) is considered when the objective function is non-convex.Therefore, in this paper, we focus on proving the convergence of our algorithms to approximate stationary points, as proving the same for a global optimal point is quite difficult.We would also like to emphasize that VQAs that have been proposed prior to our paper have not considered this notion, which is quite natural when considering non-convex objective functions.
Intuitively, stationary points are those for which the gradient of the function under consideration is equal to zero.Therefore, a stationary point can be a local maximum (including the global maximum), a local minimum (including the global minimum), or a saddle point.Consequently, when using a first-order solver such as gradient descent, one may get stuck at one of these points.However, it is often desirable to escape unwanted stationary points and move to a more favorable one, depending on the type of optimization problem (maximization or minimization).This can be achieved either through the use of higher-order solvers or by employing a noisy first-order solver [34].We focus on the latter approach in this paper, as computing higher-order derivatives can be computationally expensive compared to calculating the first-order gradient.Moreover, due to the inherent stochastic nature of variational quantum algorithms (VQAs), a first-order solver can inherently be noisy.This noise can aid in escaping unwanted stationary points and converging to a better stationary point.
The primary reason behind proposing different variational quantum algorithms for three different kinds of SDPs is the nature of the final unconstrained forms of these SDPs.In the general form of an SDP, as we do not consider any assumptions on the input matrices, the final unconstrained optimization problem turns out to be a non-convex-non-concave optimization problem (see (3)).Proving the convergence of the proposed algorithm to an approximate stationary point of this problem requires sophisticated analysis, which we leave for future work.In contrast, although the standard form of SDPs is a special case of the general form, we consider cases where these SDPs are weakly constrained.Studying the standard form of SDPs with this assumption is important because many SDPs of interest are actually large and weakly constrained.Due to this assumption, we observe that the final unconstrained forms turn out to be non-convex-concave optimization problems (see (4) and ( 6)).Non-convex-concave optimization [35,36,37,38,39] has a rich literature when compared to that of non-convex-non-concave optimization, as the latter is a harder problem than the former.Therefore, the design and convergence analysis of more sophisticated variational quantum algorithms for solving non-convex-concave optimization problems is possible.
Reformulation of General Form (GF) of SDPs-For the SDPs written as (1), we reduce this form to the following final unconstrained form expressed in terms of quantum circuit parameters θ 1 P r0, 2πs r 1 and θ 2 P r0, 2πs r 2 : where Γ Φ is the Choi operator (Definition 4) of the linear map Φ. See Section 3.1 for more details.For solving the above optimization problem, we propose a VQA called inexact Variational Quantum Algorithm for General Form (iVQAGF), in which we have two parameterized quantum circuits competing against each other to maximize/minimize the objective function, and then there is a classical optimizer that updates the parameters of these quantum circuits.
Reformulation of Equality Constrained Standard Form (ECSF) of SDPs-As stated before, for this type of SDP, we make an assumption on it being weakly-constrained, i.e., N " M .By exploiting this assumption, we design more sophisticated variational quantum algorithms, in which we need just one parameterized quantum circuit.For such SDPs, we arrive at the following final unconstrained form expressed in terms of quantum circuit parameters θ P r0, 2πs r : See Section 3.2 for more details.Throughout this paper, we use the notation to represent the expectation value of a Hermitian operator H with respect to |ϕpθqy.For solving the above optimization problem, we propose a VQA called inexact Variational Quantum Algorithm for Equality Constrained standard form (iVQAEC).We run iVQAEC on a classical computer, and at any step of the algorithm, the expectation value of a Hermitian operator is evaluated using a parameterized quantum circuit.One of our main results establishes the convergence rate and total iteration complexity of iVQAEC under the assumption of the SDP being weakly-constrained.
Reformulation of Inequality Constrained Standard Form (ICSF) of SDPs-For this case also we make the weakly-constrained assumption on SDPs.Here, we solve the dual problem instead of the primal problem.For this type of SDP, we arrive at the following final unconstrained form expressed in terms of quantum circuit parameters θ P r0, 2πs r : See Section 3.3 for more details.For solving this problem, we propose a VQA called inexact Variational Quantum Algorithm for Inequality Constrained standard form (iVQAIC).For this algorithm, we do not establish its convergence rate, which we leave for future work; instead, we prove a property of the objective function that is necessary for providing such a convergence analysis, i.e., smoothness of the objective function for a fixed ȳ.

Related Work
Recently, an approach to semidefinite programming on NISQ devices was proposed in [40], which is non-variational and called NISQ SDP solver (NSS) therein.We should note that this approach does not provide an efficient solution for a general SDP problem.Like our approach, the NSS approach also optimizes over a subset of the positive-semidefinite operator space.The NSS approach assumes the ability to prepare pure states in a set t|ψ i yu i .Using these states, the NSS approach then constructs a hybrid density matrix ansatz: X β " ř i,j β ij |ψ i yxψ j |, where the entries tβ ij u ij are stored on a classical computer.The NSS approach then transforms the original SDP into a low-dimensional SDP, which is solved by optimizing over these classical entries.Note that the NSS approach does not change the set of pure states for each iteration.As one can see, the NSS approach does not encompass the entire space of positive semidefinite operators; therefore, it is heuristic in nature, similar to ours.This means that there is no guarantee of convergence to the global optimal point.Additionally, another technique was proposed in [41] where the authors "quantized" the classical randomized cutting plane method for solving semidefinite programs.They used a quantum eigensolver subroutine in order to speedup the classical method.Their results indicate that the robustness of their approach against noise may be useful in implementing their method on NISQ devices.
Our approach here is complementary to both of these approaches because our algorithm is a variational quantum algorithm.Another important point to note here is that, unlike prior works that focus solely on solving SDPs in the standard form, we propose algorithms for solving SDPs in both the standard form and a form considered in [42,43].The latter form is prevalent in quantum information theory, as it is used to compute various relevant quantities like fidelity and trace distance.While one can convert this form to the standard form, working in the original form is more convenient, avoiding unnecessary conversion overhead.Comparing the aforementioned approaches to ours is an interesting direction for future work.

Preliminaries
In this section, we introduce some notations and definitions.
Notations: We denote the set of real and complex numbers by R and C, respectively.For a positive integer m, the notation rms denotes the set t1, . . ., mu.We use upper-case letters to denote matrices and bold lower-case letters to denote vectors (e.g., A is a matrix, and x is a vector).Let H denote a finite-dimensional Hilbert space of dimension N .This Hilbert space represents a system of n qubits, where n " rlog 2 N s.The notation LpHq denotes the set of linear operators acting on H.The set of N ˆN Hermitian or self-adjoint matrices is denoted by S N Ă LpHq.The notation S N `Ă S N denotes the set of N ˆN positive semidefinite (PSD) matrices, and D N Ă S N `denotes the set of N ˆN density matrices.Additionally, we use }¨} to represent the ℓ 2 norm of a vector, as well as the spectral or operator norm of a matrix (its largest singular value), and it should be clear from the context which is being used.Let TrrXs denote the trace of a matrix X, i.e., the sum of its diagonal terms.Let X J and X : denote the transpose and Hermitian conjugate (or adjoint) of the matrix X, respectively.The notation A ě B or A ´B ě 0 indicates that A ´B P S N `.For a multivariate function f : R n Ñ R, we use ∇f and ∇ 2 f to denote its gradient and Hessian, respectively.Let Bf p¨q Bx i denote the partial derivative of f with respect to i th component of the vector x.For a multivariate vector-valued function f : R n Ñ R m , we denote its Jacobian by J f p¨q.
Let H and H 1 be Hilbert spaces of dimensions N and N 1 , respectively.
Definition 3 (Hermiticity-preserving linear map).A linear map Φ : LpHq Ñ LpH 1 q is a Hermiticity-preserving map if ΦpXq P S N 1 for all X P S N .Equivalently, Φ is Hermiticity preserving if and only if ΦpX : q " ΦpXq : for all X P LpHq.
Definition 4 (Choi representation of a linear map).For every linear map Φ : LpHq Ñ LpH 1 q, its Choi representation Γ Φ is defined as The operator Γ Φ is also known as the Choi operator.Additionally, for every linear map Φ, the following holds @X P LpHq, Y P LpH 1 q: which is a direct consequence of the well known fact that ΦpXq " Tr 1 rΓ Φ pX J b Iqs, (11) with the partial trace over the first factor in the tensor-product space (see Ref. [43,Proposition 4.2] for a proof of (11)).
Lemma 5. Given a linear map Φ : LpHq Ñ LpH 1 q, its Choi operator Γ Φ is a Hermitian operator if and only if Φ is a Hermiticity-preserving map.
Proof.First, suppose that Φ is Hermiticity preserving.From (9), we can write where equality (a) follows from Definition 3. Now suppose that Γ Φ is Hermitian.Then it follows from (11) that ΦpXq is Hermitian if X is Hermitian, so that Φ is Hermiticity preserving.
We now recall some definitions related to the Lipschitz continuity and smoothness of a function.

Definition 6 (Lipschitz continuity).
A function f : R n Ñ R m is L-Lipschitz continuous if there exists L ą 0, such that }f pxq ´f px 1 q} ď L }x ´x1 } for all x, x 1 P R n .We say that L is a Lipschitz constant of f .For a univariate function f , suppose that the absolute value of its first derivative on an interval I is bounded from above by a positive real L, i.e., @x P I : | df pxq {dx| ď L. Then L is a Lipschitz constant of f .Lemma 7 (Lipschitz constant for a multivariate function).For a function f : R n Ñ R with bounded partial derivatives, the value L " ?n max i tsup x |Bf pxq{Bx i |u is a Lipschitz constant of f .Proof.See Appendix A.1.
As most of the objective functions that we deal with in this paper are non-convex in some parameters, it is vital to focus on local optimality rather than global optimality because finding globally optimal points of a non-convex function is generally intractable.Therefore, the notion of ϵ-stationary points is important for us.Intuitively, a point is ϵ-stationary if the norm of the gradient at that point is very small.Formally, we define an ϵ-stationary point as follows: Definition 10 (ϵ-stationary point).Let f : R n Ñ R m be a differentiable function, and let ϵ ě 0. A point x P R n is an ϵ-stationary point of f if }∇f pxq} ď ϵ.
The above definition applies to first-order stationary points.This definition is important because we use inexact first-order solvers that converge to approximate first-order stationary points, and such a definition acts as a stopping criterion.

Definition 11 (Polyak-Łojasiewicz (PL) Inequality).
A function f : R n Ñ R satisfies the PL inequality if, for some µ ą 0, the following holds for all x P R n : where f ˚is the globally optimal value of f .
In other words, the above inequality implies that every stationary point is a global minimum.

Semidefinite Programming
In this section, we recall some basic aspects of semidefinite programming [28,42].A semidefinite program is an optimization problem for which the goal is to optimize a linear function over the intersection of the positive semidefinite cone with an affine space.SDPs extend linear programs (LPs), such that the vector inequalities of LPs are generalized to matrix inequalities.
To begin with, recall the standard or canonical form of an SDP [28]: where C, A 1 , . . ., A M P S N and b 1 , . . ., b M P R. The standard form of an SDP is widely known for designing approximation algorithms for combinatorial optimization problems.
A more general form of an SDP, as considered in [42], is as follows: where C P S N , the map Φ is Hermiticity preserving (see Definition 7), B P S M , and p ˚is the optimal value of the program (16).The aforementioned form is known as the primal form of an SDP, and the corresponding dual form is given as where Y P S M `, the map Φ : is the adjoint of Φ (see Definition 2), and d ˚is the optimal value of the program (17).
The duality theorem of SDPs states that if both the primal and dual programs have feasible solutions, then the optimum of the primal program is bounded from above by the optimum of the dual program.Under a very mild condition that the primal has a feasible solution and the dual has a strictly feasible solution (or vice versa), strong duality holds; i.e., the duality gap (difference between p ˚and d ˚) is closed.In the optimization literature, this condition is well known as Slater's condition (see Theorem 1.18 of [42]).Throughout this paper, we assume that strong duality holds.

Variational Quantum Algorithms and Parameter-Shift Rule
Variational quantum algorithms are hybrid quantum-classical algorithms, designed for solving optimization tasks with an objective function of the following form: where ρ P D N is a density operator, tH k u k is a set of problem-specific Hermitian operators, i.e., H k P S N for all k, and tg k u k is a problem-specific set of functions [21,22].Additionally, each g k is a function of the expectation value of H k with respect to ρ.The corresponding optimization problem is as follows: min When the dimension N is large, the evaluation of the expectation values of Hermitian operators with respect to ρ is computationally intractable using classical algorithms.VQAs provide a quantum advantage because these hybrid algorithms attempt to circumvent this dimensionality problem by evaluating the expectation values of Hermitian operators using a quantum computer.Specifically, these algorithms utilize a parameterized quantum circuit to explore a problem-specific subspace of density operators and evaluate the expectation values with respect to these density operators.We discuss more about the quantum advantage of VQAs at the end of this section.First, let us discuss what we mean by a parameterized quantum circuit and how this circuit prepares a parameterized quantum state.
Parameterization: Let |0y RS denote the all-zeros state of systems R and S, each of which consists of n qubits.Let ρ S P D 2 n , and let U ρ RS be a quantum circuit that prepares a purification |ψy RS of ρ S when U ρ RS is applied to the initial state |0y RS .Here, the subscripts S and R are used to denote the system of interest and a reference system, respectively.VQAs simulate the space of density operators by parameterizing this quantum circuit as U RS pθq, where θ " pθ 1 , . . ., θ r q J P r0, 2πs r .Specifically, a VQA applies a parameterized quantum circuit U RS pθq to the initial pure state |0y RS to generate a parameterized pure state Let ρpθq S :" Tr R r|ψpθqyxψpθq| RS s denote the reduced density operator of |ψpθqy RS .By using the fact that TrrH k ρpθqs " xψpθq| RS pI R b H k q|ψpθqy RS and under the assumption that the parameterized circuit U RS pθq is fully expressive, the objective function in (18) and its associated optimization problem (19) can be written in terms of the parameter θ as follows: where it is implicit that H k acts on system S.
Assumption 1.We assume that the objective function in (22) is 'faithful,' which means that the minimum of this objective function corresponds to the optimal value of the problem (19).
A parameterized quantum circuit U pθq is also known as a variational ansatz, and its choice plays an important role in obtaining an approximation of the optimal value.Throughout this paper, we use the terms "parameterized quantum circuit" and "variational ansatz" interchangeably.Furthermore, there are problem-specific ansatzes as well as problem-independent ansatzes.For our case, we use a problem-independent ansatz having the following form: U pθq " U r pθ r qU r´1 pθ r´1 q ¨¨¨U 1 pθ 1 q, where each unitary U j pθ j q is written as with W j as an unparameterized unitary.
Assumption 2. We assume that the number of parameters, r, of a parameterized circuit is Oppolypnqq.
The above assumption is natural in the context of variational quantum algorithms, as we only have access to quantum circuits with short depth.Each Hermitian operator H k is arbitrary.In general, we can express a Hermitian operator as a weighted sum of tensor products of Pauli operators, i.e., H " where w i P R, σ i,j P tI, σ x , σ y , σ z u, and σ x , σ y , and σ z are the Pauli operators.From (26) we see that the expectation value of an arbitrary Hermitian operator is equal to a linear combination of the expectation values of each tensor product of Pauli operators.However, in general, this linear combination may contain many terms.
Assumption 3. We assume that the number of terms in (26) is polynomial in n; i.e., p " Oppolypnqq.
The above three assumptions are standard in the literature on variational quantum algorithms.Overall, VQAs use a quantum circuit with parameter θ to estimate the expectation value of a given Hermitian matrix, and they utilize a classical optimizer to solve the optimization problem min θPr0,2πs r Fpθq.In each round, a VQA updates the parameters of the quantum circuit according to a classical optimization algorithm.We can update these parameters according to gradient-free approaches that include Nelder-Mead [44], Simultaneous Perturbation Stochastic Approximation (SPSA) [45], and Particle Swarm Optimization [46].The main drawbacks of such gradient-free methods include slower convergence and less robustness against noise.On the contrary, if the evaluation of the gradient of a given objective function is not computationally expensive, then first-order methods like gradient descent are more suitable.
VQAs have an advantage as they do not use a classical method to evaluate gradients.Instead, they evaluate the partial derivatives of Fpθq with respective to each parameter using the same quantum circuit but with shifted parameters.It is possible to evaluate the partial derivatives on a quantum computer using the parameter-shift rule [29,30,31].Formally, we state the parameter-shift rule for the evaluation of the partial derivatives of xHy θ with respect to its j th component, i.e., θ j , as where êj is a unit vector with 1 as its j th element and 0 otherwise.From (27), we note that for the exact evaluation of the partial derivatives of xHy θ at a parameter value θ, we need to compute the expectation values of this operator at the shifted parameters exactly.However, the exact evaluation of the expectation value of an operator requires an infinite number of measurements on a quantum circuit.In reality, we perform only a limited number of measurements and then take the average of those values.Therefore, it is important to have an unbiased estimator of the expectation value.We assume that this is the case for all of our algorithms.Such an assumption is standard in the literature, as it acts as a good starting point for understanding the nature of the algorithms under such settings.
Unbiased Estimator: First, we consider a simple objective function consisting of a single expectation value term, i.e., Fpθq " xHy θ . ( Now, we define a k-sample mean unbiased estimator of this expectation value as follows: Definition 12 (k-sample mean unbiased estimator of expectation value).Given a quantum circuit U pθq with parameter θ P r0, 2πs r , we define u H k pθq as an average of k measurements of the observable H with respect to the pure state U pθq|0y.It is a k-sample mean unbiased estimator of (28) if the following holds: Next, we define an unbiased estimator of the partial derivatives of xHy θ .According to the parameter-shift rule in (27), the partial derivative of the expectation values of a Hermitian operator is a linear combination of its expectation values with shifted parameters.Setting g H,k j pθq :" it follows that g H,k j pθq is an unbiased estimator for that partial derivative because `Eru H k pθ `pπ{2qê j qs ´Eru H k pθ ´pπ{2qê j qs ˘( 32) We denote an unbiased estimator of the full gradient, i.e., ∇ θ xHy θ , as g H,k pθq.Furthermore, we evaluate the mean square error of g H,k pθq as follows, " " Var ´gH,k j pθq ¯(37) Var `uH k pθ `pπ{4qê j q ´uH k pθ ´pπ{4qê j ˘( 38) Var `uH k pθ `pπ{4qê j q ˘`Var `uH k pθ ´pπ{4qê j ˘( 39) where we denote the variance of a random variable X as VarpXq.The fourth equality follows from the definition of g H,k j pθq, given by (30).The fifth equality uses the fact that VarpX ´Y q " VarpXq `VarpY q if X and Y are independent random variables.For our case, the sample means u H k pθ `pπ{4qê j q and u H k pθ ´pπ{4qê j q are independent because they are computed using two different quantum circuit evaluations.For computing gradients in our algorithms, we take a sufficient number of samples, k, such that the mean square error is very small throughout the paper.This implies that the stochastic gradient (g H,k pθq) is almost equal to the exact gradient (∇ θ xHy θ ).
Potential quantum advantage of a VQA when N is large: According to (26), we consider a Hermitian operator H consisting of a sum of p weighted Pauli strings where p " Oppolypnqq (Assumption 3).We can measure a given Pauli string P i of arbitrary size with respect to a given quantum state in constant time.If the desired precision of the expectation value's estimate is ϵ, then we need to make Op1{ϵ 2 q repetitions of the procedure.Here the procedure consists of preparing the quantum state and then measuring the expectation value.As we have p Pauli strings, overall we need no more than Opmax i |w i | 2 p{ϵ 2 q repetitions to estimate the expectation value of H to precision ϵ.In contrast, if the evaluation of the expectation value of a Pauli string P i with respect to a general quantum state is conducted using known classical algorithms, it appears that we need Op2 n q time, as well as space.This is because the size of P i is 2 n ˆ2n , and the same goes for the size of the matrix needed to represent the quantum state.Therefore, if all the above mentioned assumptions hold, then we have an exponential advantage over the best known classical algorithms for estimating the expectation value of a Hermitian operator.A similar argument can be provided for estimating gradients because the partial derivatives of the expectation value with respect to quantum circuit parameters depend on the expectation value evaluated on shifted parameters.
3 Variational Quantum Algorithms for SDPs

General Form (GF) of SDPs
In this section, we consider the general form of an SDP, as given in (16).We write it concisely as follows: Next, we modify the above formulation as follows: " sup Equality (a) follows due to the equivalence between both problems.If we pick a primal PSD variable X that does not satisfy the constraint, i.e., ΦpXq ę B, then the inner minimization results in the value ´8.This is because, in such a case, there exists at least one negative eigenvalue of B ´ΦpXq.If we set Y " s|e i yxe i |, where |e i y is a unit vector in the negative eigenspace corresponding to that eigenvalue, then, due to the inner minimization, we can take the limit s Ñ 8.This in turn implies that inf Y ě0 tTrrpB ´ΦpXqqY su " ´8.In other words, the inner minimization forces the outer maximization to pick a feasible X and imposes an infinite penalty if chosen otherwise.Equality (b) follows by taking the infimum outside.This formulation is the Lagrangian of the original primal problem (1), where LpX, Y q :" TrrCXs `TrrBY s ´TrrY ΦpXqs (49) is a Lagrangian and Y is a dual PSD variable.Equality (c) follows from the definition of the Choi representation of the linear map Φ (see Definition 4), where Γ Φ is the Choi operator of Φ.According to Lemma 5, the Choi operator Γ Φ is Hermitian.
Equality (d) follows from the substitution: X " λρ and Y " µσ, where ρ P D N , σ P D M , λ " TrrXs, and µ " TrrY s.Now, we are interested in solving the following unconstrained optimization problem, i.e., the equality (d): This optimization problem is now expressed in terms of the expectation values of the Hermitian operators C J , B, and Γ Φ with respect to the density operators ρ, σ, and ρ b σ, respectively.

Variational Quantum Algorithm for SDPs in GF
As stated earlier, when the dimension N of the Hermitian operators is large, solving the problem (50) is generally intractable using a gradient-based classical algorithm.Therefore, we propose a variational quantum algorithm for the optimization problem (50).
First, we introduce a parameterization of the density operators, i.e., ρ and σ, by using parameterized quantum circuits.Second, we optimize the modified objective function of the problem (50) over the parameters of those quantum circuits using our variational quantum algorithm.Parameterization: Let ρ S 1 pθ 1 q be the density operator prepared by first applying the quantum circuit U R 1 S 1 pθ 1 q to the all-zeros state of the quantum system R 1 S 1 and then tracing out the system R 1 .Similarly, let σ S 2 pθ 2 q be the density operator prepared by first applying the quantum circuit U R 2 S 2 pθ 2 q to the all-zeros state of the quantum system R 2 S 2 and then tracing out the system R 2 .Here, θ 1 P r0, 2πs r 1 and θ 2 P r0, 2πs r 2 .Also, we set r 1 , r 2 " Oppolypnqq.Defining we have that Furthermore, the parameterized quantum circuits, i.e., U R 1 S 1 pθ 1 q and U R 2 S 2 pθ 2 q, are of the form shown in (24).As these quantum circuits generate purifications of density operators, the following equalities hold: where Due to the parameterization of the density operators ρ S 1 pθ 1 q and σ S 2 pθ 2 q, the problem (50) transforms into the following optimization problem: where and we optimize over the space of quantum circuit parameters θ 1 P r0, 2πs r 1 and θ 2 P r0, 2πs r 2 .As mentioned before, we assume that the objective function Gpθ 1 , λ, θ 2 , µq is faithful, which means that the global optimal value of the optimization problem (61) is equal to p ˚. Additionally, the objective function of (61) is generally non-convex as a function of the quantum circuit parameters θ 1 and θ 2 .Hence, the max-min problem (61) is a non-convex-non-concave optimization problem.Due to the fact that obtaining a globally optimal point of a non-convex-non-concave function is generally NP-hard [33, Section 2.1], we focus on finding first order ϵ-stationary points of (61).As a global optimal point is also a stationary point, if we use techniques to initialize quantum circuit parameters according to the problem at hand such that these parameters lie in the vicinity of a global optimal point, then our algorithm will converge to that point.
We propose a variational quantum algorithm for obtaining a first order ϵ-stationary point of (61).We call this algorithm inexact Variational Quantum Algorithm for General Form (iVQAGF), and the pseudocode for this algorithm is provided in Algorithm 1.It is an inexact version because we solve the subproblem involving a maximization (Step 4) to approximate stationary points instead of solving it until global optimality is reached, due to the non-convex nature of the objective function in terms of the quantum circuit parameters.Furthermore, iVQAGF is a hybrid quantum-classical algorithm, as we have two parameterized quantum circuits for estimating expectation values of Hermitian operators at any given step and a classical optimizer that updates the parameters of these quantum circuits as per the algorithm (see Figure 1).
At any step of the algorithm, the expectation value xHy θ of a Hermitian operator H is evaluated using a quantum circuit with parameter θ.Moreover, the partial derivatives of Gpθ 1 , λ, θ 2 , µq with respect to the parameters θ 1 and θ 2 depend on the partial derivatives of the expectation values of the Hermitian operators I b C J , I b B, and I b I b Γ Φ with respect to those parameters.We evaluate the partial derivative of the expectation value of a Hermitian operator with respect to the quantum circuit parameters using the parametershift rule (see (27)).We do not explicitly mention these quantum-circuit calls in the main algorithm.
# For any step, expectation values of observables and their gradients are evaluated using parameterized quantum circuits.

4:
Maximize Gp¨, ¨, θ k 2 , µ k q using a first order method such as gradient descent, where θ k 2 , µ k are fixed, to obtain pθ k`1 1 , λ k`1 q such that the following holds:

STOP and return Gpθ
estimators of expectation values of Hermitian operators, i.e., xI b C J y θ 1 , xI b By θ 2 , and Finally, as we have the unbiased estimators of all these partial derivatives, we have an unbiased estimator of the full gradient ∇Gpθ 1 , λ, θ 2 , µq of the function Gpθ 1 , λ, θ 2 , µq.

Equality Constrained Standard Form (ECSF) of SDPs
In this section, we shift our focus to the standard form of SDPs.Due to their specific problem structure and the aforementioned assumption of it being weakly constrained (i.e, N " M ), the design of more sophisticated variational quantum algorithms for solving them is possible.Moreover, we establish a convergence rate for one of the algorithms.We consider the standard form of SDPs as given in (15).Here, the Hermiticitypreserving map Φ and the Hermitian operator B have a diagonal form: where A 1 , . . ., A M P S N and b 1 , . . ., b M P R.
The equality constrained primal SDP is given as follows: where b " pb 1 , . . ., b M q J and Φ is the vector form of the original linear map, i.e., ΦpXq " pTrrA 1 Xs, . . ., TrrA M Xsq J .Taking this formulation into account, we write the optimization over the Lagrangian LpX, Y q in (49) as where LpX, yq :" TrrCXs `yJ pb ´ΦpXqq (67) and y P R M is the dual vector of the Lagrangian LpX, yq.As discussed earlier, the inner minimization with respect to y results in ´8 if X violates the constraint in (65).On the contrary, if there exists a PSD operator X that satisfies the constraints, then (66) reduces to sup Xě0 TrrCXs.This is because inf yPR M y J pb ´ΦpXqq " 0 in this case.Hence, there is an equivalence between (65) and (66), as indicated previously.
For the equality constrained problem (65), we consider the Augmented Lagrangian L c pX, yq as the objective function, which consists of a quadratic penalty term c 2 }b ´ΦpXq} 2 in addition to the terms of the original Lagrangian LpX, yq.Therefore, the optimization problem with the modified objective function can be written as follows: where c ą 0 is the penalty parameter and The method for optimization associated with the Augmented Lagrangian formulation is known as the Augmented Lagrangian Method (ALM), independently introduced in [47] and [48].This classical method involves the following three steps that iterate until convergence: 1. X k`1 :" arg max Xě0 L c pX, y k q, 2. Update the dual variable y according to y k`1 :" y k ´c pb ´ΦpXqq, 3. Update the penalty parameter c according to c k`1 :" c 0 µ k`1 , where µ ą 1.
Due to the equivalence between (65) and (68), iterating through the above steps until convergence leads to the primal optimal value.Hence, we can write ALM is a well-studied method for solving unconstrained optimization problems with convex objective functions.It is considered superior to methods that exclusively involve the original Lagrangian or exclusively involve the penalty term [49].ALM converges faster than the original Lagrangian method, as it involves a quadratic penalty term in addition to a linear penalty term of the original Lagrangian.Furthermore, it solves many issues associated with both these methods, such as ill-conditioning arising due to large values of the penalty parameter.Now, we make a mild assumption that one of the constraints is a trace constraint on the primal PSD operator X.Specifically, let A M " I and b M " λ.Hence, the constraint is TrrXs " λ.Subsequently, substituting X " λρ, where ρ P D N in (70), we obtain the following: " sup This optimization problem is now expressed in terms of the expectation values of the Hermitian operators C, A 1 , . . ., A M with respect to the density operator ρ.According to Assumption 3, the matrices C, A 1 , . . ., A M consist of polypnq terms.

Variational Quantum Algorithm for SDPs in ECSF
Solving the unconstrained optimization problem (72) using ALM is computationally intractable if the dimension N of the operators is large.Therefore, we propose a variational quantum algorithm to solve this problem.
As before, we first introduce a parameterization of the density matrix ρ and then optimize the modified objective function over the subspace of those parameters using our variational quantum algorithm.
Parameterization: Let U RS pθq be a parameterized quantum circuit with parameter θ P r0, 2πs r , and suppose that it acts on the all-zeros state of the quantum system RS and prepares a purification of the density operator ρ S pθq.The form of U RS pθq is given by (24).Moreover, the following equality holds because U RS pθq generates a purification of ρ S pθq: where H P tC, A 1 , . . ., A M u.
Due to the parameterization of the density operator ρ S pθq, the problem (72) transforms into the following optimization problem: where As mentioned before, we assume that the objective function L c pθ, yq is faithful, which means that the global optimal value of the optimization problem (4) is equal to p ˚.
We now focus on solving the optimization problem (4) using the modified Lagrangian L c pθ, yq.Note that, in general, the objective function provided above has a non-convex landscape with respect to the quantum circuit parameter θ.Also, it is linear in terms of the dual variable y.Therefore, the optimization problem (4) is a non-convex-concave optimization problem.For a convex-concave setting, the natural algorithm to consider is to first perform the maximization until a globally optimal point is reached and then update the dual variable, and repeat these steps until convergence (similar to ALM).On the contrary, for a non-convex-concave setting, solving the maximization until global optimality is reached is generally NP-hard [33,Section 2.1].Therefore, we focus on obtaining approximate stationary points of (4).As a global optimal point is also a stationary point, if we use techniques to initialize quantum circuit parameters according to the problem at hand, such that these parameters lie in the vicinity of a global optimal point, then our algorithm will converge to that point.Definition 14 (First order pϵ, cq-stationary point of ( 4)).For ϵ, c ą 0, a point θ P r0, 2πs r is a first order pϵ, cq-stationary point of (4) if there exists y P R M such that the following hold: The above condition acts as a measure of closeness to the first order pϵ, cq-stationary points of (4).Such a measure is useful to give a stopping criterion for the algorithm.
We now propose a variational quantum algorithm for obtaining a first order pϵ, cqstationary point of (4).We call this algorithm inexact Variational Quantum Algorithm for Equality Constrained standard form (iVQAEC).It is an inexact version because we solve the subproblem involving the maximization to approximate stationary points instead of solving it until global optimality is reached, due to the nonconcave nature of the objective function in terms of quantum circuit parameters.The pseudocode of iVQAEC is given by Algorithm 2.
We run iVQAEC on a classical computer, and at any step of the algorithm, the expectation value of a Hermitian operator H P tC, A 1 , . . ., A M u (i.e., xHy θ ) is evaluated using a quantum circuit with parameter θ (see Figure 2).Moreover, the partial derivative of L c pθ, yq with respect to the parameter θ depends on the partial derivatives of the expectation values of the Hermitian operators I b C, I b A 1 , . . ., I b A M with respect to θ.For evaluating the partial derivatives of the latter with respect to θ, we use the parameter-shift rule (see (27)).Moreover, we do not explicitly mention these calls to quantum circuits in the pseudocode of the algorithm.
Unbiased Estimators: As mentioned earlier in the previous section, due to the fact that evaluation of the expectation value of a Hermitian operator using a quantum circuit is stochastic in nature, we use unbiased estimators of expectation values and their corresponding partial derivatives (see (29)

Convergence Rate
In this section, we provide the convergence rate and total iteration complexity of iVQAEC in terms of the number of iterations of the for-loop at Step 3 of Algorithm 2 needed to arrive at an approximate stationary point of (75).
2: Initialization: # For any step, expectation values of observables and their gradients are evaluated using a parameterized quantum circuit.
3: for k " 1, 2, . . ., do 4: Maximize L c k p¨, y k q, where y k is kept fixed, to obtain θ k`1 such that the following holds: STOP and return where f : R d Ñ R is a continuously-differentiable non-convex function, A : R d Ñ R m is a nonlinear operator, and g : R d Ñ R is a convex function.The constrained optimization problem (79) is equivalent to the following unconstrained minimax formulation: where β ą 1, and L β px, yq is the augmented Lagrangian.As the classic ALM algorithm is suited for solving optimization problems involving convex objective functions, the authors proposed an inexact Augmented Lagrangian method, known as the iALM algorithm for obtaining an approximate first order stationary point of (80).It is an inexact version because the subproblem involving a maximization is over a non-convex function.They claim that if we use an inexact solver for that subproblem at each iteration, then their algorithm converges to an approximate first order stationary point (see Theorem 4.1 of [50]).They also provide the total iteration complexity of iALM in terms of the number of iterations of the for-loop (see Corollary 4.2 of [50]).
Our problem (75) also consists of an augmented Lagrangian term, which is non-convex in terms of the quantum circuit parameters θ, and gpθq " 0. Additionally, for solving our problem, we use an inexact solver for the subproblem (Step 5) of iVQAEC, and the update rule for modifying the learning rate η at each iteration (Step 6) is similar to their iALM algorithm.Now, the following theorem characterizes iVQAEC's convergence rate and total iteration complexity for finding an approximate stationary point of the problem (75) in terms of the number of iterations of the for-loop (Step 3).The proof of this theorem is provided by the proofs of Theorem 4.1 and Corollary 4.2 of [50].
Theorem 15 (Convergence rate and total iteration complexity of iVQAEC).For integers k 1 ě k 0 ě 2, consider the interval K :" tk 0 , . . ., k 1 u, and let tθ k u kPK be the output sequence of iVQAEC on the interval K, where θ k P r0, 2πs r for all k P K. Suppose that f pθq -xI bCy θ and Apθq -λΦpθq ´b have Lipschitz constants L f and L A .Also suppose that the function L c p¨, yq is smooth for every y P R M (i.e., its gradient is L c,y -Lipschitz continuous).Additionally, suppose that there exists ν ą 0 such that for every k P K.If an inexact solver is used in Step 5 of the algorithm, then θ k is a first order pϵ k , c k q-stationary point of (75) with for every k P K, where y max pθ 1 , y 0 , η 1 q is given as Additionally, if we use the Accelerated Proximal Gradient Method (APGM) from [51] as an inexact solver for the subproblem (Step 5), then the algorithm finds a first-order pϵ T , c T q-stationary point, after T calls to the solver, where with Q " Qpf, A, η 1 q, and µ is a constant as mentioned in Algorithm 2.
The inequality given by (82) is known as the PL inequality for minimizing }Apθq} 2 [52].This regularity condition is used in the proof of Theorem 4.1 of [50] for obtaining a bound on › › Apθ k q › › for all k P K.For our purposes, we assume that this condition holds for all θ P r0, 2πs r .Then, Theorem 15 is a direct consequence of the following statements: (A) Smoothness of L c p¨, yq: The function L c p¨, yq is smooth for every y P R M (i.e., its gradient is L c,y -Lipschitz continuous).
(B) Lipschitz continuity of f pθq and Apθq: The functions f pθq and Apθq are Lipschitz continuous with constants L f and L A , respectively.
We prove the statements (A) and (B) in Lemmas 17 and 18, respectively.Before we prove these statements, we state the following lemma that concerns Lipschitz continuity of the function h : r0, 2πs r Ñ R defined as and its gradient ∇hpθq, where O P S N .
Proof.The proof is given in Appendix A.3.
Next, we formally write the statements (A) and (B) as lemmas.
Lemma 17 (Smoothness of L c p¨, yq).For all y P R M , there exists L c,y ą 0 such that the gradient of the function L c p¨, yq : r0, 2πs r Ñ R is L c,y -Lipschitz continuous.
Proof.The proof is given in Appendix A.4.
Lemma 18 (Lipschitz continuity of f pθq and Apθq).The function f : r0, 2πs r Ñ R, where f pθq " xI b Cy θ , and the linear map Apθq : r0, 2πs r Ñ R M are L f -Lipschitz and L A -Lipschitz continuous, respectively, for some L f , L A ą 0.
Proof.The proof directly follows from the Lipschitz continuity of hpθq (Lemma 16).

Inequality Constrained Standard Form (ICSF) of SDPs
In this section, we consider the following inequality constrained primal form of SDPs: Here we take vector forms of the Hermiticity-preserving linear map Φ and the Hermitian operator B because both have a diagonal form.Furthermore, we assume that the last constraint is the trace constraint on the primal PSD variable X.As such, we set A M " I and TrrXs ď b M .Additionally, for this case also, we make an assumption that the SDPs are weakly constrained, i.e., N " M .The corresponding dual form is given as follows: where y " py 1 , . . ., y M q J is a dual variable.We can determine the dual map Φ : as follows, using the definition of the adjoint of a linear map (see Definition (8)): which implies that We derive unconstrained formulations for the primal and dual forms of inequality constrained SDPs.First, by taking the primal form of SDPs into account, we introduce slack variables and convert the inequality constrained problem (87) to the following equality constrained problem: TrrCXs where z " pz 1 , . . ., z M q J is a vector of slack variables and z 1 , . . ., z M ě 0. Now, with the problem (92) being an equality constrained problem, we can use iVQAEC (see Algorithm 2) to solve it.Next, we focus on the dual form of inequality constrained SDPs as follows: " inf Now, consider that inf tě0 tt : tI ě Hu " maxt0, λ max pHqu " max # 0, sup where H is a Hermitian operator, λ max pHq is the maximum eigenvalue of H, and the supremum is over the set D N of density operators.Using this fact, we convert the constrained problem (95) into the following unconstrained problem: where we set ȳ :" py 1 , . . ., y M ´1q J , (99) Furthermore, as the first term is independent of ρ, we have that We can use Sion's minimax theorem [53] and interchange the supremum and infimum because the set D N of density operators is compact and convex, the objective function is convex with respect to ȳ for all ρ P D N , and it is quasi-concave with respect to ρ P D N for all ȳ ě 0. Hence, The above problem is not differentiable at some parameter values due to the max operation in the objective function.In order to smoothen the sharp corners that result from this operation, we modify the above problem in the following way by introducing γ ą 0: In the above, we used the scaled version of the log-sum-exp function [28, Section 3.1.5]as an approximation to the original max function in (102).This approximation can be controlled by varying the parameter γ because the following holds for all x, y P R: By using this approximation, the objective function in (103) does not have any sudden corners where the partial derivatives can change drastically.In fact, the objective function of (103) is infinitely differentiable.It is also equivalent to the previous objective function in (102) in the limit γ Ñ 8.For all finite values of γ, the value d ˚is bounded from above by d 1 .

Variational Quantum Algorithm for SDPs in ICSF
Solving the unconstrained optimization problem (103) using a gradient based classical technique is computationally expensive if the dimension N of the operators is large.Hence, in this section, we propose a variational quantum algorithm to solve this problem.First, we introduce a parameterization of density operators using parameterized quantum circuits, in order to efficiently estimate the expectation values of the Hermitian operators C, A 1 , . . ., A M ´1.Second, we optimize the resulting objective function using our variational quantum algorithm.
Parameterization: Let U RS pθq be a parameterized quantum circuit with a parameter θ P r0, 2πs r , and let it act on the all-zeros state of the quantum system RS and prepare a purification of the density operator ρ S pθq.The structure of U RS pθq is given by (24).Furthermore, the following equality holds because U RS pθq generates a purification of ρ S pθq: Algorithm 3 iVQAIC(C, tA i u M i"1 , tb i u M i"1 , η, ϵ, γ) 1: Input: Hermitian operators pC, tA i u M i"1 q, scalars tb i u M i"1 , learning rate η ą 0, precision ϵ ą 0, constant γ ą 0.  Maximize F γ pθ, ȳk q, where ȳk is fixed, to obtain θ k`1 such that the following holds: STOP and return F γ pθ k`1 , ȳk`1 q where H P tC, A 1 , . . ., A M ´1u.
Due to the parameterization of the density operator ρ S pθq, the problem (103) transforms into the following optimization problem: where and we optimize over the space of quantum circuit parameters θ P r0, 2πs r .As mentioned before, we assume that the objective function F γ pθ, ȳq is faithful, which means that the global optimal value of the optimization problem (106) is equal to d 1 .
The objective function of (106) is in general non-convex with respect to the quantum circuit parameter θ.On the contrary, it is concave in ȳ.Hence, the optimization problem (106) is a non-convex-concave optimization problem.For such a setting, finding a globally optimal solution is generally NP-hard [33,Section 2.1].Therefore, we focus on obtaining approximate stationary points of (106).As a global optimal point is also a stationary point, if we use techniques to initialize quantum circuit parameters according to the problem at hand such that these parameters lie in the vicinity of a global optimal point, then our algorithm will converge to that point.Definition 19 (First order ϵ-stationary points of (106)).For ϵ ą 0, a point pθ, ȳq P r0, 2πs r ˆRM´1 is a first order ϵ-stationary point of (106) if the following holds: We use the above condition as a stopping criterion for our variational quantum algorithm.In order to obtain a first order ϵ-stationary point of (106), we propose a variational quantum algorithm (see Algorithm 3).We call this algorithm inexact Variational Quantum Algorithm for Inequality Constrained standard form (iVQAIC).It is an inexact version because we solve the subproblem involving the maximization to an approximate stationary point instead of solving it until global optimality is reached, due to the non-convex nature of the objective function in terms of quantum circuit parameters.Its pseudocode is provided in Algorithm 3, and it is depicted in Figure 3.
We run iVQAIC on a classical computer and utilize a parameterized quantum circuit U RS pθq with a parameter θ P r0, 2πs r for estimating the expectation values of the Hermitian operators.Moreover, the partial derivative of F γ pθ, ȳq with respect to the quantum circuit parameters θ depends on the partial derivatives of the expectation values of the Hermitian operators I b C, I b A 1 , . . ., I b A M ´1 with respect to θ.In order to compute the partial derivatives of the latter with respect to θ, we use the parameter-shift rule (see (27)).Note that we do not explicitly mention these quantum-circuit calls in the pseudocode of the algorithm.
For this case, we do not state a theorem that indicates the convergence rate for finding approximate stationary points of (106) in terms of the number of iterations in the for-loop of the algorithm.Rather, we prove a property of the objective function that is necessary for providing such a convergence analysis, i.e., smoothness of Fp¨, ȳq for a fixed ȳ ě 0. We state it formally in the following lemma: Lemma 20 (Smoothness of F γ p¨, ȳq).For all ȳ ě 0, there exists L γ, ȳ ą 0 such that the gradient of the function F γ p¨, ȳq : r0, 2πs r Ñ R is L γ, ȳ-Lipschitz continuous.
Proof.The proof is given in Appendix A.5.

Numerical Simulations
For validating our reformulations of SDPs and verifying the convergence of their respective algorithms to approximate stationary points, we randomly generated SDPs such that they contain a valid feasible region.To verify the convergence of the algorithms proposed in this paper, we focus on three cases based on the number M of constraints and the dimension N of the input Hermitian operators, and whether the problem is an equality or inequality constrained problem: 1. N « M : Here, we consider a well-known and extensively studied MaxCut problem.
First, we briefly recall what the problem is and its SDP relaxation in the next subsection.This problem is taken into account because the number of constraints and the dimension of the matrices are equal (N " M ).For this case, we evaluated the performance of iVQAGF (see Algorithm 1), as it does not make the weakly-constrained assumption on SDPs. Figure 4 shows the convergence of iVQAGF for solving randomly generated SDP instances of the MaxCut problem.Furthermore, we performed these simulations for different dimensions of the input matrices.Specifically, we considered N P t8, 16, 32u.
2. N " M : We divide this case further according to the type of constraints: (2a) equality constraints and (2b) inequality constraints.For randomly generated equalityconstrained problems, we report the performance of iVQAEC (see Algorithm 2).Similarly, we analyse the performace of iVQAIC (see Algorithm 3) for randomly generated instances of an inequality constrained problem.Here also we consider N P t8, 16, 32u. Figure 5 and Figure 6 showcase the convergence of iVQAEC and iVQAIC, respectively, for randomly generated instances of their respective problems.
We take into account Assumption 3 while creating the input Hermitian matrices.We assume that the Pauli string decomposition of these input matrices are provided beforehand.Additionally, for the analysis of iVQAGF for solving MaxCut, we assume that we are provided with the Pauli string decomposition of C J and the Choi operator of the linear map Φ, i.e., Γ Φ .
We executed our algorithms using Pennylane's Python libraries where we set QASM simulator of the Qiskit Python package as a backend.* Pennylane is an open-source Python library developed by Xanadu for differential programming of quantum computers [54].Similarly, Qiskit is an open-source package/interface developed by IBM to interact with the underlying quantum computer [27].Additionally, the QASM simulator simulates a real IBM Quantum Backend, which is actually noisy in nature due to gate errors and decoherence.We integrated Pennylane and Qiskit's QASM simulator using the Pennylane-Qiskit plugin.
We use the Strong Entangling Layers template of Pennylane as our variational ansatz, where each layer consists of Oppolypnqq single-qubit rotations and entangling gates.We then repeat this layer Oppolypnqq number of times, where n is the number of qubits.Therefore, the overall gate complexity is Oppolypnqq.
In order to assess the convergence of our methods to a globally optimal point, we initialize the quantum circuit parameter such that it lies in the convex region of that globally optimal point.We then report the results and their associated analyses on how well our algorithms perform on the QASM noisy simulator and compare the results with a noiseless simulator.We also report the time complexity (number of iterations of the for-loops) of our algorithms, i.e., how fast our algorithms converge to an actual globally optimal point.
First, for the sake of completeness, we recall the definition of a cut and a MaxCut of a given graph.

Definition 21 and MaxCut).
A cut is a bi-partition W of the vertex set V of a graph G " pV, Eq, where |V | " N .An edge pi, jq P E is part of the cut set if its vertices i and j lie in separate partitions.A MaxCut is the largest cut possible of a graph G.
The problem of finding a MaxCut of a graph can be formulated as a Quadratic Integer Program (QIP) [55]: subject to x i P t´1, 1u; @i P V. (109) Here, the optimal value of the above optimization is the optimal cut size (i.e., the number of the edges in the MaxCut).This quadratic integer program is in general computationally intractable [56].However, there exist LP and SDP relaxations for the above program [57].
The above mentioned LP formulation is not good enough as it is a 1/2-approximation to the original problem.It is well known that this can be extended to an SDP formulation in which an algorithm proposed in [57] gives a 0.879-approximation of the original QIP in (109): subject to X ě 0, X ii " 1, @i P V. (111) For our case, we numerically simulate iVQAGF for solving the aforementioned SDP relaxation of a MaxCut problem.Recasting the constraints of the above MaxCut SDP as constraints of the general form of SDP, we can write them as ΦpXq " B, where and A i " |iyxi|.Therefore, the Choi operator for this linear map is given as Additionally, the objective function can be written as TrrCXs, where C " L{4 and L is the Laplacian matrix for the graph G.
The Pauli decomposition of the above Choi operator consists of N " 2 n Pauli strings.This gives the impression that we need to compute 2 n expectation values in order to compute the expectation value of Γ Φ .However, all these Pauli strings can be constructed from just I and σ z operators, because it is clearly a diagonal matrix.This implies that all the Pauli strings are mutually commuting, and their expectation values can be estimated simultaneously [58].Therefore, we just need to compute a single expectation value to evaluate the expectation value of Γ Φ .
The shaded regions in Figure 4, 5, and 6 signify the variance in the values for different runs (specifically 20) of the algorithm.The y-axis measures the difference between the cost function value at each iteration and the actual optimal value evaluated with the CVXPY package [59].The x-axis shows the number of outer iterations of algorithms, i.e., the time needed to converge to the actual optimal value.
The numerical simulations demonstrate that all three algorithms indeed converge to their respective optimal values approximately.This numerical evidence suggests that all three proposed algorithms work well in practice.In other words, convergence to the optimal parameters in the presence of noise showcases noise resilience of our variational quantum algorithms.

Conclusion
In this paper, we proposed variational quantum algorithms for solving semidefinite programs.We considered three constrained formulations of SDPs, which were first converted to unconstrained forms by employing a series of reductions.When the dimension N of the input Hermitian operators of SDPs is large and for these unconstrained forms, the computation of the objective function's gradient is difficult when using known classical techniques.To address this problem, we utilized parameterized quantum circuits to estimate these gradients.We also established the convergence rate and total iteration complexity of one of our proposed VQAs.Finally, we numerically simulated our variational quantum algorithms for different instances of SDPs, and the results of these simulations provide evidence that convergence still occurs in noisy settings.
The estimation of the gradients using parameterized quantum circuits is stochastic in nature.In this paper, we assumed that we have unbiased estimators of these gradients, and the variance of these estimators is also small.Therefore, it remains open to study the effect of the variance of these estimators on the convergence rate of our algorithms.where the inequality paq follows from the triangle inequality, inequality pbq follows from (114)-(115), and the last inequality pcq follows from the fact that }¨} 1 ď ? 2 }¨} for the two-variable case.Therefore, L " ? 2 maxtL X , L Y u is a Lipschitz constant for f .The proof given above for two variables can be easily extended to a function f with an n-variable input, using the fact that }¨} 1 ď ?n }¨}, where L " ?n max i tL i u i is a Lipschitz constant given in terms of (121) Here, x " px 1 , . . ., x n q J .

A.2 Proof of Lemma 8
By hypothesis, each component f i of the vector-valued function f : R n Ñ R m is L i -Lipschitz continuous.Hence, the following holds for all x, x 1 P R n and i P t1, . . ., mu: Now consider the following: " Hence, L " a ř m i"1 L i is a Lipschitz constant for the vector-valued function f .

A.3 Proof of Lemma 16
According to Lemma Bθ i in terms of the linear combination of hpθ `pπ{4qê i q and hpθ ´pπ{4qê i q.
A.4 Proof of Lemma 17 The last inequality follows from (142) and (137).We see that L i γ, ȳ is bounded from above by a positive number for a fixed ȳ.Therefore, the Lipschitz constant L γ, ȳ is also bounded from above by a positve number because according to Lemma 8, we have L γ, ȳ " b ř r i"1 L i γ, ȳ.

Figure 2 :
Figure 2: This figure depicts the iVQAEC algorithm in which we utilize one parameterized quantum circuit, i.e., U RS pθq.

2 :
Initialization: ȳ1 " 0.# For any step, expectation values of observables and their gradients are evaluated using a parameterized quantum circuit.

Figure 3 :
Figure 3: This figure depicts the iVQAIC algorithm in which we utilize one parameterized quantum circuit, i.e., U RS pθq.

Figure 4 :Figure 5 :
Figure 4: Convergence of iVQAGF for three randomly generated MaxCut-SDP instances with different numbers of vertices in the graph: N P t8, 16, 32u.

Figure 6 :
Figure 6: Convergence of iVQAIC for three separate cases of randomly generated inequality constrained semidefinite programs with nonempty feasbile regions: N P t8, 16, 32u.