Quantum Algorithm For Testing Convexity of Function
Nhat A. Nghiem Tzu-Chieh Wei



nhatanh.nghiemvu@stonybrook.edu
Department of Physics and Astronomy, State University of New York at Stony Brook, Stony Brook, NY 11794-3800, USA
C. N. Yang Institute for Theoretical Physics, State University of New York at Stony Brook, Stony Brook, NY 11794-3840, USA


Department of Physics and Astronomy, State University of New York at Stony Brook, Stony Brook, NY 11794-3800, USA
C. N. Yang Institute for Theoretical Physics, State University of New York at Stony Brook, Stony Brook, NY 11794-3840, USA
§ ABSTRACT 
Functions are a fundamental object in mathematics, with countless applications to different fields, and are usually classified based on certain properties, given their domains and images. An important property of a real-valued function is its convexity, which plays a very crucial role in many areas, such as thermodynamics and geometry. Motivated by recent advances in quantum computation as well as the quest for quantum advantage, we give a quantum algorithm for testing convexity of polynomial functions, which appears frequently in multiple contexts, such as optimization, machine learning, physics, etc. We show that quantum computers can reveal the convexity property superpolynomially faster than classical computers with respect to number of variables. As a corollary, we provide a significant improvement and extension on quantum Newton's method constructed in earlier work of Rebentrost et al [New J. Phys. 21 073023 (2019)]. We further discuss our algorithm in a broader context, such as potential application in the study of geometric structure of manifold, testing training landscape of variational quantum algorithm and also gradient descent/Newton's method for optimization. 

Quantum Algorithm For Testing Convexity of Function
    Tzu-Chieh Wei
    Received 16 July 2024; accepted 04 September 2024
=====================================================





§ INTRODUCTION

Quantum computers hold great promise to solve difficult computational problems that lie beyond the reach of classical computers. The underlying power of quantum computers is due to two intrinsic properties of quantum mechanics: superposition and entanglement. Tremendous efforts have been made to exploit the potential of quantum computers in various contexts. Some early pioneering works <cit.> showed that quantum computers could probe properties of a blackbox function with a single query usage. The breakthrough work of Shor <cit.> showed that quantum computers can factorize a given integer number superpolynomially faster than their classical counterpart, which has been recently improved in <cit.>. Grover <cit.> later showed that a quadratic speedup is achieved for unstructured database search. Further quantum speedup has been showcased in a wide array of problems, such as simulating quantum systems <cit.>, solving linear systems <cit.>, supervised and unsupervised learning  <cit.>, principle component analysis <cit.>, topological data analysis <cit.>, learning from experiments <cit.>, etc. As a whole, these developments have ignited an exciting view towards the application of quantum computers, as well as triggering efforts in experimental realization of fault-tolerance devices as to bring quantum computation steps closer to reality. 
  

Given these successes, the question of the repertoire of tasks quantum computers can excel in is still worthy of pursuing. An important model that has been central to investigating quantum advantage is the so-called blackbox model. In such a model, we have access to the blackbox with an unknown structure, that computes some Boolean functions, e.g., accepting some input strings and then outputting a Boolean variable (0 or 1). As the structure is unknown to us, the goal is to extract properties of such functions with minimum resource, e.g., the number of access to the corresponding blackbox. In fact, some early works <cit.> showed that quantum computers can reveal hidden properties of given functions using a minimum number of queries much smaller compared to classical counterparts. More relevant to our work, Jordan <cit.> considered the numerical gradient estimation problem given the blackbox that computed some function explicitly and showed that a single query is sufficient to reveal the gradient at a given point, up to some desired accuracy. Inspired by such a line of pursuit, we consider the potential advantage of quantum computers in the topic of functional analysis, where our object of interest is a (multivariate) function, which is a very basic object in mathematics. The intriguing point is that a function can possess very rich analytical properties, and thus the problem is very appealing to explore from a computational perspective. A particular property that we focus on in our work is the convexity, which captures the shape of the function in some domain (see Fig. <ref> and <ref> below).

A convex function has a peculiar feature that a local minimum is also a global minimum (see Fig. <ref> for simple illustration). Additionally, if a function is convex in some domain, then a minimum is easily obtained, e.g., by the gradient descent method, which makes it very useful in optimization areas.  

In this work, we aim to tackle the challenge of testing the convexity of some polynomial function. We begin with a simple case, which is a homogeneous polynomial of even degree (to be defined later) and subsequently, building upon such homogeneous polynomial, we generalize the construction to arbitrary polynomial type. The kind of homogeneous polynomial of even degree was also considered in <cit.>, where the authors proposed the quantum gradient descent and quantum Newton's method for finding local minima. In fact, our work is quite inspired by them, as the structure of a given function makes it simpler to compute two important quantities: the gradient and Hessian. Building upon <cit.>, we show that the ability to obtain the analytical form of Hessian translates into the ability to test convexity by examining the sign of Hessian's spectrum, and that quantum computers can achieve that goal with superpolynomially less cost than their classical counterpart in relative to the number of variables included in the given polynomial function. Along the way, as a corollary, we show how to construct the Hessian more efficiently than the original method in <cit.>, thus providing a significant improvement on quantum Newton's method that also appeared in <cit.>. First, our improved quantum Newton's method work on arbitrary polynomial, instead of homogeneous one of even degree. Second, our method reduces the complexity dependence on (inverse of) error tolerance from polynomial to polylogarithmic. The complexity dependence on the degree of given function is also reduced by a power of 4. Additionally, at each step of the Newton's method, the number of copies of (block encoding of) temporal solution required in our work is polynomially less than (by a power of 5) that of <cit.>. We further mention three potential subjects that might be useful for our framework: differential geometry, variational quantum algorithm and gradient descent as well as Newton's method for finding minima of objective function.  


The structure of the paper is as follows. First, in Section <ref>, we provide an overview of our objectives, with the corresponding assumptions and criteria for convexity in Section <ref>. Our main algorithm is outlined in details in Section <ref>. Remarks and further discussions are given in Section <ref>, where we show that our method can be generalized to an arbitrary polynomial and showcase some potential applications of our method in the context of differential geometry and variational quantum algorithm, as well as improving the quantum Newton's method originally introduced in <cit.>. Appendix  <ref> contains the necessary recipes that underlie our work. Appendix <ref> provides a proof of Lemma <ref>.



§ OVERVIEW



 §.§ Overview of the Problem and Assumptions


We consider a multivariate function f: ℝ^n ⟶ℝ which is a homogeneous polynomial of an even degree 2p (p ∈ℤ). Let x = (x_1, x_2, ..., x_n); as shown in <cit.>, such a polynomial admits a tensor algebraic decomposition:

    f(x) = 1/2x^T ⊗⋯⊗x^T A  x⊗⋯⊗x,

where A is an s-sparse matrix of dimension n^p × n^p. 

Throughout this work, we assume that the oracle's access to A is given in a similar fashion to that in <cit.>. The knowledge of the polynomial is obtained by the oracle's access to matrix A, which can be formally decomposed as:

    A = ∑_α = 1^K A_1^α⊗⋯⊗ A_p^α,

where each A_i^α is a matrix of size n× n (i= 1,2,...,p) and K is some natural number. While the above general decomposition is not required explicitly in our work (i.e., we do not assume the oracle has direct access to the submatrices A_i's), the expression is helpful as the important quantities, such as the gradient and Hessian, admit analytical forms. More specifically, the gradient can be written as:

    ∇ f(x) = D(x) ,

where D is specified as 

    D() =  ∑_α=1^K∑_j=1^p ( ∏_i=1, i≠ j^p  ^T A_i^α) A_j^α,

and the Hessian matrix of function f can be written as: 

    () = 2 ∑_α =1^K  ∑_j,k = 1, j≠ k^p   ∏_i=1, i≠ j,k^p (^T A_i^α)   A_k^α^T A^α_j + D().

 The special formulation above provides us an alternative way to compute the Hessian and gradient at a given point  as follows. Denote ^T as ρ_x, then we can write the gradient as

    D() = _1,2,...,p-1( M_D  ( (^T)^⊗ p-1⊗_n  )  ),

where 

    M_D = ∑^p_m=1 M_m,

and each M_m (m=1,..,p) can be obtained from A via permutation of entries as explained below. Recall that 

    A = ∑_α = 1^K A_1^α⊗⋯⊗ A_m^α⊗⋯⊗ A_p^α,

and M_m is defined similarly to A except having the p-th matrix swapped with the m-th matrix, i.e.,

    M_m = ∑_α = 1^K A_1^α⊗⋯⊗ A_p^α⊗⋯⊗  A_m^α.

In a similar fashion, the Hessian can be expressed alternatively as:

    () = _1,2,...,p-1( (M_H + M_D)   ( (^T) ^⊗ p-1⊗_n )   ),

where 

    M_H = 2 ∑_j ≠ k^p Θ_jk,

and each Θ_jk (for j,k = 1,...,p) can be obtained from A via permutation of matrices, e.g., the j,k-th matrices are swapped to p-1 and p-th ones, respectively, 

    Θ_jk = ∑_α = 1^K A_1^α⊗⋯⊗ A_p-1^α⊗⋯⊗ A_p^α⊗⋯⊗ A_j^α⊗ A_k^α.

Therefore, the oracle access to A allows us to obtain entries of all Θ_jk respectively. 



We remark that in the above, we discuss homogeneous polynomials of even degree as a primary target due to its simplification of gradient and Hessian formulation. As mentioned in <cit.>, inhomogeneity can be inserted. For example, given a homogeneous polynomial f_ homo, one can multiply it with c^T where c is some n-dimensional vector to obtain an inhomogeneous polynomial. In the following, we describe our framework with a homogeneous polynomial of even degree and then provide a generalization in Section <ref>. 

Let ||·|| denotes the operator norm. As specified in <cit.>, we have that the norm || D || ≤ p ||A|| and || || ≤ p^2 ||A||. Without loss of generalization, we can set the norm of A, ||A|| = 1, so that p^2 ||A|| ≤ p^2, and thus we can guarantee that ||D|| < p and |||| ≤ p^2. Throughout this work, we assume this condition holds for convenience, as rescaling by some factor does not change the nature of the problem, e.g., the convexity of function and particularly the (asymptotical) running time of our method. 



 §.§ Dissecting Convexity of Function




We are interested in the convexity of f in some domain 𝒟⊂ℝ^n. By trivially redefining the function (e.g., with a coordinate shift), we can choose 𝒟 to be a hypersphere with radius 1 for simplicity. The following criterion is a well-known result in mathematics and can be found in many standard literature. 


Criterion: Over some domain 𝒟, if the Hessian matrix of a given function f is positive-semidefinite at every point, then the function f is convex. 


The application of the above criterion is straightforward. We compute the Hessian matrix of f at chosen points in the region and check its spectrum. The positive-semidefiniteness of a matrix is equivalent to all eigenvalues being non-negative. Therefore, the sign of the spectrum is necessary and sufficient for revealing convexity. As we mentioned, the Hessian matrix  can be written down explicitly, which is convenient for spectral analysis. As the dimension of  grows with the number of variables (there are n variables), finding the full spectrum of  is computationally demanding. Furthermore, the only information that matters is the sign of the smallest eigenvalue, as we only need to see if it is positive or not. Therefore, we propose to find the minimum eigenvalue of . If such an eigenvalue is non-negative, then all other eigenvalues are non-negative, which implies that  is positive-semidefinite and hence, the corresponding f is convex within the domain 𝒟. In the following section, we construct a quantum algorithm to first find and then verify the sign of the smallest eigenvalue, thereby dissecting the convexity of the function f, given the details of f and related assumptions from the previous section Sec. <ref>. 



§ QUANTUM ALGORITHM 


We begin with a remark that all details regarding relevant definitions as well as useful tools are provided in Appendix <ref>. Here proceed to describe our main result.



 §.§ Constructing Hessian  for a Single Point


In Sec. <ref>, we first introduce the general idea behind our work: finding the minimum eigenvalue of . Our first task is to produce the block encoding of  that can be used for further analysis. In order to achieve this goal, we recall the following important formulation:

    () = _1,2,...,p-1( (M_H+M_D)   ( (^T)^⊗ p-1⊗_n )   ),

where M_H = 2 ∑_j ≠ k^p Θ_jk 
 
and M_D = ∑_m=1^p M_m. The oracle access to entries of each M_m and H_jk allows us to use Lemma <ref> to produce the ϵ-approximated block encoding of M_H/(sp^2) and M_D/(sp) (where s is the sparsity of A) in time: 
𝒪(  p(plog(n) + log^2.5(1/ϵ) )  ). 
Thus, it is quite simple to obtain the block encoding of M_D/(sp^2) from M_D/(sp)  by multiplying with a factor 1/p (see Lemma <ref>). Lemma <ref> then allows us to obtain the ϵ-approximated block encoding of (M_H+M_D)/(2sp^2). 

The next recipe that we need includes the following simple relations.

    _1 (A ⊗ B) (^T ⊗) = (^T A ) B, 
          (^T ⊗) (A ⊗ B) (^T ⊗) = ^T ⊗ (^T A ) B.

Using these properties plus the tensor structure of M_D and M_H, we can show that:

    (^T) ^⊗ p-1⊗_n  (M_H + M_D)  (^T) ^⊗ p-1⊗_n    = (^T) ^⊗ p-1⊗().

The reason comes from the fact that

    (^T)^⊗ p-1 = ^⊗ p-1 (^T)^⊗ p-1

and that

    (^T)^⊗ p-1⊗_n  (M_H+M_D)  ^⊗ p-1⊗_n = ().

For now, we assume that we have a block encoding of (^T), as subsequently, we will show how to produce such a block encoding and generalize it to deal with multiple points chosen from the domain 𝒟, as required by the convexity criterion. Lemma <ref> allows us to use p-1 block encodings of (^T) and a trivial block encoding of _n to obtain the block encoding of (^T)^⊗ p-1⊗_n. Then Lemma <ref> (combined with the simple relations derived above) yields the block encoding of (^T)^⊗ p-1⊗()/(2sp^2). 



 §.§ Constructing “multi-points” Hessian




  §.§.§ For homogeneous polynomial of even degree

We remark that the above construction yields the block encoding of the operator (^T)^⊗(p-1)⊗ ()/(2sp^2), which contains the tensor product of ^T and Hessian matrix  at a given point . As it will become clearer later, in order to take advantage of quantum parallelism for analyzing the Hessian spectrum at multiple points, we need to adjust the above procedure to construct what we call a “multi-point” Hessian. 


For further clarity and to avoid confusion with the above construction, we set the following notation for subsequent discussions. Let 𝒩 be the number of points of consideration; _i (∈ℝ^n) is the i-th point; the corresponding Hessian of f evaluated at _i is then (_i). The first goal is to produce the block encoding of the following operator ⊕_i^𝒩 (_i_i^T)^⊗p-1⊗(_i)/(2sp^2), which has the matrix representation as:

    1/2sp^2 [  (_1_1^T)^⊗p-1⊗(_1)                   ⋯                   ⋯                   ⋯;                   ⋯  (_2_2^T)^⊗p-1⊗(_2)                   ⋯                   ⋯;                   ⋯                   ⋯                   ⋯                   ⋯;                   ⋯                   ⋯                   ⋯ (_𝒩_𝒩^T)^⊗p-1⊗(_𝒩 );                     ].

Now, we outline the procedure that produces the desired block encoding. First, we have the following lemma:


   Given block encoding of (M_H+M_D)/(2sp^2) (as constructed in Sec. <ref>), the block encoding of operator ⊕_i^𝒩 (M_H + M_D)/(2sp^2) can be prepared with extra 𝒪(1) cost.

Proof: To show this, we need to add an extra register of dimension 𝒩. The resulting tensor product of _𝒩 and unitary block encoding of (M_H+M_D)/(2sp^2) produces the block encoding of _𝒩⊗ (M_H+M_D)/(2sp^2), which is exactly ⊕_i^𝒩 (M_H + M_D)/(2sp^2) by a simple algebraic property. 



Next, we introduce the following crucial recipe:


    The block encoding of the operator ⊕_i=1^𝒩 (_i_i^T)^⊗ p-1⊗_n can be prepared in time 𝒪(p(log(n)).

Proof: To prove the above lemma, we first consider a log(n𝒩) qubits system which resides a Hilbert space ℋ of dimension n𝒩. Let U be some unitary operator. Once acting on the basis state | 0⟩_ we obtain the state |ϕ⟩, i.e.,

    U  | 0⟩_ = |ϕ⟩.

Let |ϕ⟩ = (x_1,x_2,...,x_n𝒩)/C where

    C = √(∑_k=1^n x_k^2 )

is the normalization factor. In this paper, we work in the real regime, i.e., x_i ∈ℝ for all i = 1,2,..., n.  If we break such a vector into  parts, and denote each part as _j (j=1,2,..., ), e.g., _j = (x_(j-1)n +1,x_(j-1)n+2, ..., x_(j-1)n +n). Note that given such notation, the normalization factor C is essentially equal to C = ∑_i=1^ |_i|^2 where |.|^2 refers to the usual l_2 norm of a vector. It is then straightforward to decompose |ϕ⟩ as: |ϕ⟩ = 1/C∑_j=1^𝒩|j⟩⊗_j. Before proceeding, we remark on the unitary U that generates the desired state |ϕ⟩. We are interested in  points {_i}_i=1^. If we choose these  points classically, which means that we know their coordinates plus the norms respectively, then we can use the well-known amplitude encoding method <cit.> to load these entries into a quantum state, resulting in a quantum circuit U of depth 𝒪(log(n)). On the other hand, if we choose U to be a random unitary circuit, then a constraint is imposed on the coordinates of |ϕ⟩, which implies C=1. In this case, the set of points {_i}_i=1^ must have their norms summing up to 1, which may reduce the number of points considered in a given domain in one go. 

Now we append another ancillary system having dimension  initialized in | 0⟩_, to obtain the state | 0⟩_⊗|ϕ⟩ = 1/C∑_j=1^𝒩| 0⟩_⊗|j⟩⊗_j. Using CNOT gates to copy the second register to the first one, i.e.,

    1/C∑_j=1^𝒩| 0⟩_⊗|j⟩⊗_j ⟶1/C∑_j=1^𝒩|j⟩_⊗|j⟩⊗_j.

If we trace out the first register (with subscript ), we obtain the state (1/C^2)∑_j=1^|j⟩⟨j| ⊗ _j_j^T (we recall that it should be _j_j^†, but since we work in the real regime, ^T and ^† are identical). Lemma <ref> allows us to prepare the (exact) block encoding of (1/C^2) ∑_j=1^|j⟩⟨j| ⊗ _j_j^T (in complexity 𝒪(log (n) ), whose matrix representation is as follows:

    1/C^2∑_j=1^|j⟩⟨j| ⊗ _j_j^T = 1/C^2[  _1_1^T       ·       ·       ·;       · _2 _2^T       ·       ·;       ·       ·       ·       ·        ;       ·       ·       ·    __^T ].

If C is greater than 1, the factor C^2 from the above representation can be removed using the amplification technique (basically uniform singular value amplification from <cit.>), with a further complexity 𝒪(C^2). In reality, if we prefer to choose  points {_i}_i=1^ uniformly within the hypersphere, then we can expect that C^2 = ∑_i=1^ |_i|^2 ≤, which means that the complexity of the above step can be 𝒪() (for C ≥ 1). 

For C smaller than 1, one cannot use the amplification method <cit.>.  Therefore, the factor C cannot be removed by amplification. (Of course, this can be avoided using different or more points.) For now, we continue the construction with C being greater or equal to 1 (which means it is removed from the above equation). Subsequently, we will return to the case C being smaller than 1 and show that the structure of a homogeneous polynomial allows us to factor out some power of C, resulting in a different expression for the Hessian evaluated at given points. Hence, the final complexity is different (from the case C> 1) by a factor of power of C. 

Define |_i⟩≡_i/|_i|, e.g., the normalized vector. We have that _i_i^T = |_i|^2 |_i⟩⟨_i|. Our goal is to transform the above operator into its square root, e.g., for each i, we aim to transform |_i|^2 |_i⟩⟨_i|⟶ |_i||_i⟩⟨_i|. In order to achieve such a goal, we recall two results from <cit.> and <cit.>:


    Given δ, ϵ∈ (0, 1/2], c∈ (0,1] and let f(x) = 0.5 x^c. There exists an even/odd polynomial of degree 𝒪( 1/δlog( 1/ϵ ) ) such that
    
    ||P-f||_ [δ,1] ≤ϵ,  ||P||_[-1,1]≤ 1.



[<cit.> Theorem 56]

Suppose that U is an
(α, a, ϵ)-encoding of a Hermitian matrix A. (See Definition 43 of <cit.> for the definition).
If P ∈ℝ[x] is a degree-d polynomial satisfying that


  * for all x ∈ [-1,1]: |P(x)| ≤1/2,

then, there is a quantum circuit Ũ, which is an (1,a+2,4d √(ϵ/α))-encoding of P(A/α) and
consists of d applications of U and U^† gates, a single application of controlled-U, and 𝒪((a+1)d)
other one- and two-qubit gates.

Choosing c= 1/2, by using polynomial from Lemma <ref> as an approximation to 0.5√(x) plus Lemma <ref>, we can obtain the following (ϵ-approximated) transformation (remind that for all i, _i_i^T ≡ |_i|^2 |_i⟩⟨_i|):

    [  _1_1^T       ·       ·       ·;       · _2 _2^T       ·       ·;       ·       ·       ·       ·        ;       ·       ·       ·    __^T ]⟶[ 0.5 |_1|  |_1⟩⟨_1|                  ·                  ·                  ·;                  ·  0.5 |_2| |_2⟩⟨_2|                  ·                  ·;                  ·                  ·                  ·                  ·                   ;                  ·                  ·                  ·     0.5 |_| |_⟩⟨_| ]

Note that we can not remove the factor 0.5 by amplification, as the norm needs to be less than or equal to 1/2. Moreover, Lemma <ref> admits an approximation of a given positive power function on the interval [δ,1]. To apply such a result to our context, we need to make sure that our interval is appropriate, which means that δ needs to be less than or equal to min_i { |_i|^2 }_i=1^. In order to find such a minimum, we can use Lemma <ref> to the left-hand side of the above equation in the case that the points are unknown. However, if we classically pick these points, we know their norms, including the minimum.


In Sec. <ref>, we have mentioned a very simple way to prepare the block encoding of the identity matrix of any dimension. Suppose we have a block encoding of _n(p-1), then Lemma <ref> allows us to prepare the block encoding of ∑_i=1^|i⟩⟨i|⊗ |_i||_i⟩⟨_i|⊗_n^⊗ p-1. We note that with 2log(n) SWAP gates, we can swap the order of |_i||_i⟩⟨_i| and any _n among p-1 such ones. Therefore, it takes p-1 further steps to achieve all of them, e.g., block encoding of operators of the form 

    ∑_i=1^|i⟩⟨i|⊗_n ⊗ |_i| |_i⟩⟨_i|⊗⋯⊗_n,  ∑_i=1^|i⟩⟨i|⊗_n ⊗_n ⊗ |_i||_i⟩⟨_i|⊗⋯⊗_n,  ..., ∑_i=1^|i⟩⟨i|⊗_n ⊗_n ⊗⋯⊗ |_i||_i⟩⟨_i|⊗_n.


Lemma <ref> yields the block encoding of their products, and it is easy to note that their product is ∑_i=1^|i⟩⟨i|⊗ ( |_i||_i⟩⟨_i|)^⊗ p-1⊗_n, which is exactly ⊕_i=1^𝒩 (|_i||_i⟩⟨_i|)^⊗ p-1⊗_n by a simple tensor algebraic property. The complexity of this step is 𝒪(plog(n) 1/|_min|), where |_min| = min_i { |_i|^2 }_i=1^. To summarize what we have so far, we state the following:

 Assuming that C = √(∑_i=1^ |_i|^2)≥ 1. An ϵ-approximated block encoding of ⊕_i=1^𝒩 (|_i||_i⟩⟨_i|)^⊗ p-1⊗_n can be prepared in time complexity 𝒪( p 1/|_min|log(n) log(1/ϵ)  ). 


We are now ready to construct the “multi-point” Hessian, which is straightforward by combining the methods outlined in <ref>, Lemmas <ref> and <ref> and the following simple property of matrix multiplication:

    [ A_1   0;   0 A_2 ]·[ B_1   0;   0 B_2 ] 
        =  [ A_1B_1      0;      0 A_2B_2 ],

which holds for any higher dimension N, i.e.,

    ⊕_i=1^N A_i  ⊕_i=1^N B_i = ⊕_i=1^N A_i B_i.

More specifically, from Lemma <ref>, we have  the block encoding of ⊕_i=1^ (M_H+M_D)/(2sp^2). Using Lemma <ref> to obtain the block encoding of 

    ⊕_i=1^𝒩( 0.5|_i||_i⟩⟨_i|) ^⊗ p-1⊗_n   · ⊕_i=1^(M_H+M_D)/2sp^2· ⊕_i=1^𝒩( 0.5 |_i| |_i⟩⟨_i|)^⊗ p-1⊗_n,

which is exactly 

    ⊕_i=1^( 0.5|_i||_i⟩⟨_i|)^⊗ p-1⊗_n  ·(M_H+M_D)/ 2sp^2·( 0.5|_i||_i⟩⟨_i|)^⊗ p-1⊗_n    = 1/4^p-1⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗(_i)/2sp^2.

Due to the fact that, for all i, |_i||_i⟩⟨_i| = |_i⟩_i^T (basically we absorb the norm to the ⟨_i| to obtain _i^T), then we use property (<ref>) to obtain the block encoding of ⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗(_i)/(2sp^2). We remark a subtlety that, for any i, 0.5|_i||_i⟩ is essentially proportional to |_i||_i⟩, which means that one can `absorb' that factor 0.5 into the calculation of Hessian. This is a particular property of homogeneous polynomial (see Sec. <ref>), whereby a rescale of the given input ⟶λ (for λ∈ℝ) would result in a rescaling of Hessian, i.e., (λ) = λ^2(p-1)(). 

At the beginning of this section, we mentioned two cases: the normalization factor C ≥ 1 or C < 1. For C <1, we cannot remove the factor C in Eqn (<ref>); the aforementioned property of homogeneous polynomial then allows us to handle the case C ≤ 1 in a simple manner, as we treat the point _i/C as a scaled point _i ⟶_i/C, which means that eventually there will be a factor C^2p-2 absorbed, e.g, in the above formulation, we would have the following operator: 

    1/(4C^2)^p-1⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗(_i)/2sp^2.

    
One might wonder that for the case C ≥ 1 if we did not use amplification to remove the factor C, then we would end up having the same form as above. It does not seem to be an issue for homogeneous polynomials, as discussed above, due to the difference in the Hessians being just a scaling factor. However, there are two reasons. First, our method is aimed at dealing with polynomials of arbitrary type, as we will generalize subsequently. This means that the homogeneous property will not hold; the homogeneous polynomial is just a base on which we can build and achieve the generalization conveniently. Second, as mentioned previously, for C < 1, we can not use the amplification method to remove the factor C. Besides that, the subsequent strategy that we will use to dissect convexity (in Section <ref>) is tracking the value of some operators that have the form as the above. If there is an extra factor C (being greater than 1), we will need to choose the error tolerance to be smaller (by a factor of C^2p-2) to reveal the correct eigenvalue of the desired Hessian, which would result in a substantial running time. 


With the above operator (for both cases C ≥ 1 and C < 1), we are able to generalize it to polynomials of arbitrary type. The following section shows how to achieve our goal based on what we have obtained so far.



  §.§.§ Generalization to Polynomial of Arbitrary Kind


Previously, we used the particular form of f, which is a homogeneous polynomial of even degree. Now we generalize our method to deal with arbitrary polynomials, or more specifically, monomials as also mentioned in <cit.>, which include homogeneous polynomials of odd degree and inhomogeneous polynomials. According to <cit.>, an inhomogeneous function can be given below by inserting an extra factor into a homogeneous one:

    f() = ∑_q=1^P-1 (c_q^T )  ∏_k=1^q-1 (^T B_k q ).

We recognize that the term ∏_k=1^q-1 (^T B_kq) is in fact an alternative expression for a homogeneous polynomial of even degree 2(q-1), and c_q^T is a -dependent coefficient that adds inhomogeneities. Therefore, all terms in the above summation share a similar form, and the function f(x) is simply adding them. Since the derivative of a sum of functions is the sum of the derivative of the constituting functions, we consider the following part separately g()  =  (c^T )  ∏ (^T B ) where we already ignore the subscript and treat with full generalization, e.g., with arbitrary order of the polynomial. The partial derivative:

    ∂ g/∂ x_m = c^T ∂∏ (^T B )  /∂ x_m  + ∏ (^T B )  ∂ c^T /∂ x_m.

Now we take further partial derivative:

    ∂^2 g/∂ x_n ∂ x_m = c^T ∂^2 ∏ (^T B )  /∂ x_n ∂ x_m  + ∂∏ (^T B )  /∂ x_m ∂ c^T /∂ x_n + ∂∏ (^T B )  /∂ x_n ∂ c^T /∂ x_m.

Denote ∏ (^T B )  ≡ h(). Since the term ∏ (^T B ) is a homogeneous polynomial of even degree, we know its gradient (composed of partial derivatives, see Eqn. (<ref>)) and Hessian matrix (composed of partial derivatives of second order, see Eqn. (<ref>)), therefore, the above equation implies that the Hessian of g is generally expressed as:

    (g()) = (c^T )  (h() ) +  ∇ h()  c^T + c  ∇ h() ^T.

Since h() is a homogeneous polynomial of even degree, its gradient and Hessian can be computed using known technique <cit.> (see further Sec. <ref>, more specifically Eqn (<ref>) and Eqn (<ref>)). The only difference we need to consider is the contribution from c, which accounts for the inhomogeneities part. How to deal with the extra term c^T depends greatly on what kind of access we have to c. We now outline a solution in the case where c is generated by some unitary U_C, e.g., U_C | 0⟩ = |c⟩≡ c (assuming further that |c|=1). From the unitary U_C, using Lemma <ref> we have the block encoding of cc^T.  


For now, we assume to work in the regime where the normalization factor C ≥ 1. The first goal is to produce the block encoding of some operator that includes ∇ h() c^T as a tensor component (as in Sec. <ref>). Since h() is the regular homogeneous part, the gradient operator can be computed according to Sec. <ref>, i.e., there exists a procedure similar to what was outlined in Sec. <ref> and Sec. <ref> that produces the block encoding of the following operator  

    ⊕_i=1^1/4^p-1(|_i⟩⟨_i|)^⊗ p-1⊗1/sp^2∇ h(_i) _i^T.

More specifically, from equation (<ref>), if we ignore the M_H operator, we then obtain a similar property for D():

    (^T) ^⊗ p-1⊗_n · M_D · (^T) ^⊗ p-1⊗_n    = (^T) ^⊗ p-1⊗ D().

Consequently, we can use the result of Section <ref> and <ref> to obtain the block encoding of the operator ⊕_i=1^( 0.5 |_i||_i⟩⟨_i|)^p-1⊗_n. Then we use the block encoding of operator M_D/sp^2 plus lemma <ref> to obtain the block encoding of their products, which is basically 

    1/4^p-1⊕_i=1^  (|_i⟩⟨_i|)^⊗ p-1⊗1/sp^2 D(_i).

From Sec. <ref>, we have the operator ⊕_i=1^_i_i^T (the construction above Lemma 3). If we use Lemma <ref> to construct the block encoding of ⊕_i=1^_i _i^T ⊗_n^⊗ p-1. Then, with 2log(n) SWAP gates, we can swap the first and last register, e.g., obtaining ⊕_i=1^_n^⊗ p-1⊗_i_i^T. Then we use Lemma <ref> to obtain the block encoding of:

    1/4^p-1⊕_i=1^  (|_i⟩⟨_i|)^⊗ p-1⊗1/sp^2 D(_i)   ·⊕_i=1^_n^⊗ p-1⊗_i_i^T    = 1/4^p-1⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗1/sp^2  D(_i)_i_i^T 
       =  ⊕_i=1^1/4^p-1(|_i⟩⟨_i|)^⊗ p-1⊗1/sp^2∇ h(_i) _i^T.

Note that we used the property of gradient operator of homogeneous even degree function D(_i)_i = ∇ h(_i). 


From the block encoding of cc^T, it is trivial to produce the block encoding of ⊕_i=1^_n^⊗ p-1⊗ cc^T. Denote c^T _i = _i^T c ≡β_i. We then use Lemma <ref> to produce the block encoding of the multiplied operator:

    ⊕_i=1^ ( |_i⟩⟨_i|)^⊗ p-1⊗∇ h(_i) _i^T/4^p-1 sp^2·⊕_i=1^_n^⊗ p-1⊗ cc^T    = ⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗β_i/4^p-1 sp^2∇ h(_i) c^T.

In order to remove the factor β_i from the above formulation, we use the following procedure. 

In Appendix <ref>, we prove the following thing:

   
   Given the block encoding of cc^T, then it is possible to construct the block encoding of the diagonal matrix ℬ with entries ℬ_ij = β_iβ_j δ_ij where β_i = _i^T c.

Recall that for the homogeneous part h(), procedure in Sec. <ref> obtains the following operator:

    1/4^p-1⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗1/2sp^2(_i).

Thus, we use Lemma <ref> and Lemma <ref> to obtain the block encoding of:

    1/4^p-1⊕_i=1^|_i⟩⟨_i|^⊗ p-1 ⊗β_i^2 (_i) /2sp^2.

Recall that we have the block encoding of the operator 

    ⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗β_i/4^p-1 sp^2∇ h(_i) c^T.

Using Lemma <ref> to insert an extra factor 1/2 into the above operator, i.e., we obtain: 

    ⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗β_i/4^p-1 2 sp^2∇ h(_i) c^T.

We remark that, the transpose of the block encoding of the above operator is the block encoding of: 

    ⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗β_i/4^p-1 2 sp^2 c ∇  h(_i)^T


We then use Lemma <ref> to obtain the block encoding of a sum of two operators above, 

    ℙ   = 1/3(  1/4^p-1⊕_i=1^|_i⟩⟨_i|^⊗ p-1 ⊗β_i^2 (_i) /2sp^2 +  ⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗β_i/4^p-1 2 sp^2∇ h(_i) c^T +  
       ⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗β_i/4^p-1 2 sp^2 c ∇  h(_i)^T) 
       = 1/31/4^p-1sp^2(  ⊕_i=1^β_i |_i⟩⟨_i|^⊗ p-1 ⊗_inho(_i)    ),

where _inho(_i) refers generally to the Hessian evaluated at _i of the given inhomogeneous function (see derivation in Eqn. <ref>). Our goal is to remove the factor β_i (for all i) in the above operator. Recall that from Lemma <ref> we have block encoding of ℬ that contains β_i^2 (for i=1,2,...,) on the diagonal. Note that we also have block encoding of ℬ⊗_n  ≡⊕_i=1^β_i^2 _n (which is trivial to obtain using block encoding of identity matrix _n plus Lemma <ref>). We then can use the following polynomial approximation of the negative power function from <cit.>:

    Let δ, ϵ∈ (0,1/2], c>0 and let f(x) = δ^c/2 x^-c, then there exists a (could be even or odd) polynomial P such that || P - f(x) ||_δ,1≤ϵ, ||P||_-1,1≤ 1/2. The degree of the polynomial P is 𝒪(max [1,c]/δlog(1/ϵ)  ).

To apply the above lemma to the operator ⊕_i=1^β_i^2 _n, we need to know the lower bound, i.e., the minimum of {β_i^2}_i=1^, denoted as β_min which can be estimated using Lemma <ref>. Then, we  use the above Lemma (with c = 1/2) with Lemma <ref> to ϵ approximately transform the operator 

    ⊕_i=1^β_i^2 _n ⟶⊕_i=1^√(β_min)/2β_i_n.

The complexity of this step is 𝒪( 1/β_minlog(1/ϵ ) ). Then, we can use Lemma <ref> to obtain the approximation of the block encoding of 

    ⊕_i=1^√(β_min)/2β_i_n ·ℙ  =  √(β_min)/2· 3 · 4^p-1sp^2  (  ⊕_i=1^|_i⟩⟨_i|^⊗ p-1 ⊗_inho(_i)    ).

Finally, in order to remove the factor √(β_min)/6, we can use  the amplification method <cit.> with further complexity 𝒪(6/√(β_min)  ) = 𝒪(1) as it doesn't scale as any input parameter. In the end, we obtain the following operator:

    1/4^p-1 sp^2(  ⊕_i=1^|_i⟩⟨_i|^⊗ p-1 ⊗_inho(_i)    ),

which has a similar form to Eqn. (<ref>). The above construction was discussed for the case C ≥ 1. For C < 1, all the steps are essentially the same (as in Section <ref>), except that there will be an appearance of factor C at the end, i.e., we would obtain the following operator: 

    1/(4C^2)^p-1 sp^2(  ⊕_i=1^|_i⟩⟨_i|^⊗ p-1 ⊗_inho(_i)    ).


Before moving further, we make a simple remark that for either case, C ≥ 1 and C < 1, the operators of interest (e.g., the above one) only differ by a factor of C^2p-2. Now, we are ready to tackle our main objective, which is to dissect the convexity of a given function of arbitrary type. 



 §.§ Testing Positive-semidefiniteness


Remind that we are interested in the spectrum, or more specifically, the minimum eigenvalue of Hessian  at some given point . From Sec. <ref> we have that |||| ≤ p^2 at any point , which means the eigenvalues of /p^2 lie within (-1,1). Given such a range of eigenvalues of , the shifted matrix (_n - /p^2)/2 would have a spectrum lying in the range (0,1). The reason why we change the attention to (_n - /p^2)/2 is because we can apply the improved quantum power method proposed in <cit.> to find its maximum eigenvalue. 


Given the block encoding of some positive-semidefinite matrix A whose eigenvalues are ∈ (0,1), then its largest eigenvalue can be estimated up to additive accuracy δ in time 

    𝒪( T_A/δ( log(1/δ) + log(n)/2)  ),

where T_A is the complexity of producing the block encoding of A.


To see how it applies to our main problem, let λ_min denotes the minimum eigenvalue of /p^2. If λ_min < 0 (being negative) then (1 - λ_min)/2 > 1/2, and that (1 - λ_min)/2 is the maximum eigenvalue of (_n - /p^2)/2. Therefore, we can track sign of λ_min from estimating minimum eigenvalue of (_n - /p^2)/2. 


From the previous section we have obtained the block encoding of 1/4^p-1⊕_i=1^|_i⟩⟨_i|^⊗ p-1 ⊗(_i)/2sp^2. The last recipe that we need to complete our algorithm is the following. First, suppose {A_i}_i=1^M is a set of Hermitian operators, then the eigenspace of ⊕_i=1^M A_i is simply the union of the eigenspace of all {A_i}_i=1^M <cit.>. Second, the eigenspace of ⊗_i=1^M A_i is the tensor product of the eigenspace of each {A_i}, which means that the eigenvalues of ⊗_i=1^M A_i is the multiplication of eigenvalues of contributing matrices. From the second property, we can claim that for any i, the eigenvalues of (1/4^p-1) |_i⟩⟨_i|^⊗ p-1 ⊗(_i)/(2sp^2) is the eigenvalues of (1/4^p-1) (_i)/(2sp^2). The reason is that, each operator |_i⟩⟨_i| is a projector and hence, its only eigenvalue is 1. Therefore, from the first property, we have that 

    maximum eigenvalue of ⊕_i=1^1/21/2s 4^p-1( _n -  |_i⟩⟨_i|^⊗ p-1 ⊗(_i)/p^2) = max__i{_n - (_i)/p^2/4s 4^p-1}_i=1^.


One may wonder why the right-hand side is important. We remind that if the function is convex in the given domain 𝒟, then its Hessian matrix does not have non-negative eigenvalues at all points in 𝒟. Therefore, if max__i{(_i)  } is non-negative, then the Hessian is semi-definite at all  points of consideration. As we pointed out from the beginning of this section, the condition of min__i{(_i)  } being non-negative is equivalent to min__i{  (_n - (_i)) / 2 } being no less than 1/2. 

Fortunately, given the block encoding of (1/4^p-1)⊕_i=1^|_i⟩⟨_i|^⊗ p-1 ⊗(_i)/2sp^2, it is simple to use lemma <ref> (plus a trivial block encoding of identity matrix of dimension n) to obtain the block encoding of their summation, which is exactly: 

    ⊕_i=1^1/4 sp^2 4^p-1( _n +  |_i⟩⟨_i|^⊗ p-1 ⊗(_i) ).

Then, we can find its minimum eigenvalue using recent work of <cit.>, as we mentioned at the beginning of this subsection, e.g., Lemma <ref>. 

We finally remark that we are estimating the minimum eigenvalue of the operator 

    ⊕_i=1^1/4s 4^p-1( _n -  |_i⟩⟨_i|^⊗ p-1 ⊗(_i)/p^2),

which contains the extra factor 4s 4^p-1. What we actually want is the maximum eigenvalue of the operator 

    ⊕_i=1^1/2( _n -  |_i⟩⟨_i|^⊗ p-1 ⊗(_i) /p^2),

which corresponds exactly with 

    max__i{_n - (_i)/2p^2}_i=1^.

Therefore, we need to use Lemma <ref> with an adjusted multiplicative error, e.g., by choosing

    δ⟶δ/2s 4^p-1.


Summary of the quantum algorithm procedure. For convenience, we provide key points of our framework (which means we leave out technical details) plus corresponding complexity at each step.


∙  We begin with a set of points {_i}_i=1^ (each _i ∈ℝ^n and |_i| ≤ 1) of interest. Define C^2 = ∑_i=1^ |_i|^2 as normalization factor. 


∙  Load the above points into quantum state |ϕ⟩ = 1/C∑_i=1^|i⟩⊗_i, then construct the block encoding of the following operator

    1/C^2[  _1_1^T       ·       ·       ·;       · _2 _2^T       ·       ·;       ·       ·       ·       ·        ;       ·       ·       ·    __^T ].

The complexity of the above step is 𝒪( log(n ).


∙  Break into two cases C ≥ 1 and C < 1. First, consider C ≥ 1, then use amplification <cit.> to remove the factor C^2. The complexity for amplification step is 𝒪(C^2 log(n)). Then we 
construct the following operator 

    [ 0.5 |_1|  |_1⟩⟨_1|                  ·                  ·                  ·;                  ·  0.5 |_2| |_2⟩⟨_2|                  ·                  ·;                  ·                  ·                  ·                  ·                   ;                  ·                  ·                  ·     0.5 |_| |_⟩⟨_| ].

The complexity of this step is 𝒪( C^2/_minlog(n)   ) where _min≡min{ |_i|^2 }_i=1^.


∙  Using the above operator, we construct the following operator 

    1/2^p-1[ ( |_1| |_1 ⟩⟨_1 |  )^⊗ p-1⊗_n                                                                                          ;                               ( |_2| |_2 ⟩⟨_2 |  )^⊗ p-1⊗_n                                                            ;                                                                                         ⋯                              ;                                                                                               ( |_| |_⟩⟨_ |  )^⊗ p-1⊗_n ].

The complexity upon this step is 𝒪( p C^2/_minlog(n)  ). 


∙ Use oracle access to A (defined in section <ref>) to construct the block encoding of (M_H + M_D)/2sp^2. The complexity of this step is 𝒪( p^2 log(n)  ). 


∙  Employing mathematical property of homogeneous polynomial (see <ref>), we construct the block encoding of 

    1/4^p-1⊕_i=1^ (|_i⟩⟨_i|)^⊗ p-1⊗(_i)/2sp^2

that contains our Hessian  (for homogeneous polynomial). 


∙  Generalize the above construction to a polynomial of arbitrary type. The complexity of this step (including everything from the beginning) is

    𝒪( C^2/_minlog(n)  + p^2 log(n) ).


∙  Construct the block encoding of the operator 

    ⊕_i=1^1/4 s 4^p-1( _n -  |_i⟩⟨_i|^⊗ p-1 ⊗(_i/p^2) ).


∙  Find minimum eigenvalue of the above operator by using Lemma <ref> (with accuracy being scaled ϵ/4sp^2 4^p-1, and infer the convexity from such eigenvalue. More specifically, if the minimum eigenvalue is less than 1/2, then the function is not convex. Otherwise, it is convex. The complexity of this step is:

    𝒪( 4^p-1 sp^2 ( log(n) + log( 4^p-1 sp^2 ) ) ( p C^2/_min (log(n) + p^2 log(n) ) ).


∙  Finally consider the case where the normalization factor C < 1. Repeat the same procedure, which results in almost the same complexity, i.e.,

    𝒪( (4C^2)^p-1 sp^2 ( log(n) + log( 4^p-1 sp^2 ) ) ( p C^2/_min (log(n) + p^2 log(n) ) ).


We remark that in the above summary, we did not take into account the error term ϵ that appears in multiple steps, such as amplification, encoding matrices M_H,M_D, and taking the square root of the operator. To dissect the convexity, we can ignore the error-dependence factors because we only care about the sign of the minimum eigenvalue instead of a real estimation. Therefore, one can set the error to be some constant. Furthermore, the complexity dependence on the error is polylogarithmic, which is efficient. The generalization to polynomials of arbitrary type is carried out in Section <ref>. A subtle detail about the above running time is that the scaling depends on |_min|, which is the minimum norm of the length of chosen points, as well as C, which is the normalization factor. Let's say we choose  points to be distributed uniformly in the domain of interest (hypersphere), then we can assume that |_min| ∼𝒪(1) (being some constant) and hence 1 ≤ C^2 ≤. So the factor C^2/|_min| ∈𝒪(), which is linear in . For C being smaller than 1, which means that we choose  points having very small norm (at least, for most of them). Then we can assume that |_min| is smaller than 1/, which means that 1/|_min| is greater than , but no greater than 𝒪(). All in all, for both cases where C ≥ 1 or C < 1, the factor C^2/|_min|  is 𝒪(), which is linear in .  


Now we are ready to state our main result formally:


    Let an objective function f (of arbitrary type), with certain assumptions defined as in Sections <ref> and <ref>. Over the domain 𝒟⊆ (-1,1)^n, let 𝒩 be the chosen sample points, denoted as {_i }_i=1^ for convexity testing. Let C^2 = ∑_i=1^ |_i|^2 (≤))  and _min = min_i { |_i|^2 }_i=1^. If C ≥ 1, then the quantum algorithm outlined above can reveal the positive-semidefiniteness of corresponding Hessian in time 
    
    𝒪( 4^p-1  sp^2   ( log(n) + log(4^p-1 2sp^2))  ( p  (log(n)) +  p^2 log(n) ),

    where α is some bounded constant. In particular, if C < 1, then the running time is: 
    
    𝒪( (4C^2)^p-1 sp^2   (log(n) + log(4^p-1 2sp^2))  ( p  (log(n)) +  p^2 log(n) ).




Potential advantage

To see the advantage, we take a look at how classical computers can solve the above problem. The Hessian matrix is computed as

    () = _1,2,...,p-1( (M_H + M_D)   ( ρ_x^⊗ p-1⊗_n )   ),

which involves matrix operation of size n^p. Each operator M_H and M_D is computed from matrix A which takes further time 𝒪( pn^p + p^2 n^p ). The Hessian is then evaluated at  points to find the sign of minimum eigenvalue, which results in  the final complexity 𝒪( ( pn^p + p^2 n^p)).

The above running time clearly suggests that the quantum algorithm can test convexity superpolynomially more efficiently than its classical counterpart with respect to the number of variables n, meanwhile a bit slower in relative to the number of sample points . Therefore, the quantum method is very efficient when dealing with high-dimensional cases as well as high-degree polynomials.


An interesting question from the above classical procedure is: what if we follow the same routine using a quantum computer for all points? It means that one can run the quantum framework to check the spectrum of Hessian at each point , and repeat the procedure for  points. In such a case, the quantum procedure would be very similar to our prior construction, except that we do not care about other - 1 points, more specifically, if we take a look at Eqn. <ref> and pay attention only to the top-left corner, e.g., the operator _1 _1^T (the factor C is treated in a similar manner as the above construction), which means that we treat this point _1 as a point of interest. Then, the same procedure can be carried out to first build the squared operator:

    [ _1_1^T      ·      ·      ·;      ·      ·      ·      ·;      ·      ·      ·      ·       ;      ·      ·      ·      · ]⟶[ 0.5 |_1|  |_1⟩⟨_1|                  ·                  ·                  ·;                  ·                  ·                  ·                  ·;                  ·                  ·                  ·                  ·                   ;                  ·                  ·                  ·                  · ].

Then, from the right-handed side operator, one proceeds to build the block encoding of ( 0.5 |_1||_1⟩⟨_1|)^⊗ p-1⊗_n (using the routine below Equation (<ref>)). Then, one uses the same properties of Hessian (see Eqn. (<ref>)) to obtain the following:

    ( 0.5 |_1||_1⟩⟨_1|)^⊗ p-1⊗_n  · (M_H+M_D)/2sp^2·( 0.5 |_1||_1⟩⟨_1|)^⊗ p-1⊗_n = 1/4^p-1|_1⟩⟨_1|^p-1⊗(_1)/2sp^2.

One can see that the above formula is essentially the top-left block corner of Eqn. (<ref>), which is quite obvious. Then, one can use the same Lemma <ref> to reveal the minimum value, which indicates the positivity of Hessian  at given point _1. One repeats the process for  different points, which results in asymptotically the same complexity. The difference is that, as we need to store  different eigenvalues (of Hessian at  points), the memory usage is required to be as much as . Meanwhile, the quantum process that we have outlined during this work requires 𝒪(log()) qubits to handle the same task, which is more effective in terms of memory usage but sharing the same running time. 


§ REMARKS AND DISCUSSION



In this section, we discuss our algorithm in a larger context, showing the potential application of our result in multiple directions. 



  §.§.§ From Hessian to Curvature and Geometric Structure of Manifold

As we mentioned, the Hessian of a function at a given point encodes the local geometric structures of such a function, e.g., convexity. While the sign of eigenvalues of Hessian can reveal the convexity, the magnitude of eigenvalues can actually reveal how curved it is (see Fig. <ref> for a simple example in 3D). Figure <ref> features two simple surfaces in 3D, which are a plane and a sphere. For the plane, it is easy to compute the Hessian of z(x,y) = -(a/c)x - (b/c)y + d/c in this case, which is simply a 2 × 2 matrix of 0 entries, which implies that the surface is not curved at all, i.e., exactly match the plane shape. On the other hand, for a sphere, z =√(r^2 - x^2 -y^2) so the Hessian matrix is not as straightforward as a plane, but as we can see, the function z is convex on the lower half of plane x-y, but not on the upper half. 


In fact, studying the underlying structure of a manifold using functions defined on the manifold is a fundamental aspect of differential geometry. Thus, from this angle, our work suggests a potential aid of a quantum computer in the study of the geometric structure of a manifold, given that a manifold is locally Euclidean and arbitrary smooth function of some variables can be approximated by certain polynomials in some domains.




  §.§.§ Variational Quantum Algorithm

A popular topic of recent progress in quantum computation is the variational quantum algorithm (VQA) <cit.>. One of the reasons that make variational quantum algorithms so appealing is that typically, they require low-depth circuits, which is very suitable in the near-term era. VQA has shown its success in many areas such as combinatorial optimization <cit.>, supervised learning <cit.>, etc. A common strategy in the context of VQA eventually boils down to the minimization problem: 

    min_θ f(θ) = ⟨0| U(θ)^† O U(θ) |0⟩,

where U(θ) refers to some variational circuit, e.g., a circuit composed of rotational gates with adjustable angles. A common method to optimize the above quantity is gradient descent, where the above observables are usually defined as a cost function, and its gradient is computed classically. Then, the parameters are updated iteratively by tuning corresponding rotational gates. Optimizing a quantum circuit is apparently not easy, as a phenomenon called the barren plateau can occur <cit.>, which prevents the efficient training of quantum circuits. In particular, training general variational quantum circuits is even NP-hard <cit.>.

Theoretically, the cost landscape of the above function is also a factor that affects the minimization. If the domain is convex, then any initial randomization falling into that domain can lead to a minimum, which implies that the optimization is efficient. While the exact expression for f(θ) might not necessarily be a polynomial, we remark that over some domains that are sufficiently small, arbitrary functions can be approximated by some polynomials, e.g., Taylor series. Therefore, a visible strategy that we see from our work is that we can consider multiple domains and check the convexity. Given that the task of dissecting convexity can be done superpolynomially fast in the quantum realm, it turns out to be quite useful in this case as we can scan the optimization landscape efficiently with respect to the number of parameters. The challenge then only arises as to how to represent f(θ) with proper polynomials, which allows computable parameters in a form that we assume in our work. At first glance, this challenge seems impossible for a random circuit U, but might not be so for a structured U that has been employed in several contexts, such as the so-called QAOA ansatz <cit.>. We leave this question as a motivation for future exploration of the application of our convexity tester. 



  §.§.§ Improving Quantum Newton's Method

As a striking corollary of our method is the generalization to arbitrary polynomial as well as major improvement on the resource requirement, upon the quantum Newton's method developed in <cit.>. In <cit.>, we already made progress on quantum gradient descent, and in fact, the method outlined in this work shares certain similarities with what in <cit.>, but the scope is different. 

To begin with, we recall some aspects of quantum Newton's method, which is basically a modified gradient descent method. We first follow the same notation and assumption from <cit.>. In such a problem, we are given an objective function f: ℝ^n ⟶ℝ, which is a homogeneous polynomial of even degree as we defined in Section <ref>. The goal is to find the point ∈ℝ^n at which f() is minimum. A standard method for this kind of optimization problem is gradient descent, which is an iterative method in that one begins with a random solution _0 and performs the update iteratively as follows:

    _t+1 = _t - η∇ f(_t).

As mentioned in Section <ref>, such particular form of f admits an analytical form for the gradient, as ∇ f() = D(), which was used in <cit.> (and improved in <cit.>) to construct the quantum process carrying out the above iteration procedure. As also mentioned in the same work <cit.>, Newton's method modifies directly upon the gradient descent method by taking account of the curvature of f, i.e., the update rule is as follows:

    _t+1 = _t -  η ^-1(_t) ∇ f(_t).

Roughly speaking, at a given time step t-th, the method in <cit.> takes multiple copies of _t, use oracle access to simulate exp(-iM_D t) and exp(-i M_H t), to construct the gradient and Hessian (via relation in Equation <ref> and <ref>). Then they perform the subtraction of vectors by using extra ancilla and Hadamard gates, to obtain _t+1. If one wish to obtain a normalized version of temporal solution _t+1, |⟩_t+1, then by measuring the ancilla and post-select on ancilla being |0⟩, one achieves the goal, as in <cit.>. Here, we lift such requirement and consider the problem in a more general manner, that is obtaining a temporal solution written in a form similar to “density operator”, _t _t^T for any given time step t.    


From Section <ref>, <ref> (particularly equation (<ref>)) and most importantly section <ref>, we have the block encoding of (^T)^⊗ p-1⊗()/(2sp^2) where  refers to the Hessian of polynomial of arbitrary kind. What is lacked is the gradient of polynomial of arbitrary kind, which was not constructed before and neither in the previous work <cit.>, so here we first fill this gap. Recall that from section <ref>, we have that for a general polynomial:

    f() = ∑_q=1^P-1 (c_q^T )  ∏_k=1^q-1 (^T B_k q ).

Since the derivative of a sum is equal to sum of derivative of each term within the summation, then for simplicity, as similar to what we did in section <ref>, we consider each constituent of the above summation, which has the form g()  =  (c^T )  ∏ (^T B ) ≡ (c^T ) h() (where we have defined ∏ (^T B ) ≡ h(). Then by chain rule, it is simple to see that:

    ∂ g/∂ x_m   = (c^T) ∂ h()/∂ x_m  + h() ∂ c^T /∂ x_m
       = (c^T) ∂ h()/∂ x_m + h() ∂∑_i=1^n c_j x_j /∂ x_m
       = (c^T) ∂ h()/∂ x_m  + h() c_m

where x_m is the m-th variable of  and c_m is the m-th entry of c (note that 1 ≤ m ≤ n and n is th dimension of ). As the gradient of a function is composed of partial derivative of such function with respect to all variables, we have that:

    ∇ g() = (c^T ) ∇ h() + h() c

Since h() is a homogeneous even degree polynomial, its gradient admits an explicit analytical expression. Recall that from section <ref> and property in equation <ref>, we have that:

    (^T) ^⊗ p-1⊗_n  ( M_D)  (^T) ^⊗ p-1⊗_n    = (^T) ^⊗ p-1⊗ D(h())

where D(h()) refers directly to the fact that it is the gradient operator of h() evaluated at the point . 

Recall that via oracle access to A, we have ϵ-approximated block encoding of M_D/(sp). Then we can use lemma <ref> to construct the block encoding of 

    (^T) ^⊗ p-1⊗_n  M_D/sp (^T) ^⊗ p-1⊗_n = (^T) ^⊗ p-1⊗D(h())/sp

We also have block encoding of cc^T (by assumption, see section <ref>). Lemma <ref> allows us to construct the block encoding of (^T)  (cc^T) = (^T c) c^T. Then using lemma <ref> using block encoding of , we can obtain block encoding of ⊗ (^T c) c^T. Then lemma <ref> allows us to construct the block encoding of 

    (^T) ^⊗ p-1⊗D(h()) (^T c) c^T /sp = (^T) ^⊗ p-1⊗ (^T c) ∇ h() c^T/sp

Using lemma <ref> again with block encoding of ⊗ (c^T) c^T (which is the transpose of ( c^T)  c^T), we obtain the block encoding of 

    (^T) ^⊗ p-1⊗ (^T c) ∇ h() c^T  (c^T) c^T  /sp  =    (^T) ^⊗ p-1⊗ ( c^T)^2 ∇ h() ^T  /sp

where we have use c^T c = 1 and ^T c = c^T due to the real regime that we work on. Now we handle the term h() c. Since h() is homogeneous even degree, it has the familiar form:

    h() = 1/2⟨|^⊗ p A |⟩^⊗ p

As we have oracle access to A, we can construct the block encoding of A/2s. From block encoding of ^T, it is trivial to obtain the block encoding of (^T)^⊗ p using lemma <ref>. We have that:

    (^T)^⊗ pA/2s (^T)^⊗ p = (^T)^⊗ p-1⊗h() /s^T

Now we use the block encoding of ⊗ cc^T plus lemma <ref> to construct the block encoding of

    ( ⊗ cc^T)  ( (^T)^⊗ p-1⊗h() /s^T   ) = (^T)^⊗ p-1⊗ (c^T) h() c^T   /s

We use lemma <ref> to add a scaling of 1/p to the above term, e.g, obtaining the block encoding of (^T)^⊗ p-1⊗ (c^T) h() c^T   /ps. Then we use lemma <ref> to construct the block encoding of:

    1/2(  (^T) ^⊗ p-1⊗ ( c^T)^2 ∇ h() ^T  /sp  +(^T)^⊗ p-1⊗ (c^T) h() c^T   /ps) 
       = 1/2(  (^T)^⊗ p-1⊗^T c/sp ( (^T c) ∇ h()^T + h()c ^T  )   )  
       =  1/2(  (^T)^⊗ p-1⊗^T c/sp ∇ g()^T )

Before moving further, we remind that we have the (ϵ-approximated) block encoding of the following operator:

    (^T )^⊗ p-1⊗()/2sp^2  and  1/2(  (^T)^⊗ p-1⊗^T c/sp ∇ g()^T )


Now we have enough recipe to deal with quantum Newton's method. Let us recall a lemma from <cit.> (note that in their context, xx^T is exactly ^T in our case):

     Given the block encoding of (x x^T)^⊗ p-1⊗D(x)/ps (with complexity T), then it is possible to obtain the block encoding of D(x)/p in time 𝒪( γ^2(p-1) s T_D ), where γ is some bounded constant.

The above form coincides with what we have right above. Therefore, using the same procedure, we can obtain the block encoding of ()/p^2 in similar time 

    𝒪(γ^2p-2 s  ( p^2log(n) + plog^2.5(1/ϵ))).
 
and we also obtain the block encoding of (^T c) ∇ g()^T/p in time

    𝒪(  γ^2(p-1)  s ( p^2log(n) + plog^2.5(1/ϵ))).    )

The extra factor x^T c can be estimated up to small accuracy by using routine in appendix <ref>, setting = 1 and using amplitude estimation. The concrete procedure will be elaborated in appendix <ref>. Then it is removed from the above operator using amplification method, with further 𝒪( 1/(x^T c)) = 𝒪(1) complexity. Therefore, we obtain the ϵ-approximated block encoding of ∇ g()^T/p in complexity 

    𝒪(  γ^2(p-1)  s ( p^2log(n) + plog^2.5(1/ϵ)))


The minimum eigenvalue (in magnitude) of ()/p^2 can be found exactly via the convexity testing framework (see section <ref> and <ref>), which is required to perform inversion of () in order to execute quantum Newton's method. The inversion of ()/p^2, or more precisely, the transformation from ()/p^2 ⟶^-1()/Γ, where Γ is roughly κ, where κ is roughly a reciprocal of minimum eigenvalue of () (according to <cit.>), is then carried out using popular methods <cit.>, with running time

    𝒪(γ^2p-2 2s  ( p^2log(n) + plog^2.5(1/ϵ) )  κ  polylog(κ/ϵ) ),

Therefore, it is sufficient to complete the improved quantum Newton's method. More specifically, at t-th time step, we are presented with _t _t^T. In particular, _t _t^T is related to _t+1_t+1^T via the following relation, that was used in <cit.>:

    _t+1 (_t+1)^T    =  (_t  - η^-1(_t) ∇ g(_t))( _t  - η^-1(_t) ∇ g(_t)  )^T 
       =   (_t  - η^-1(_t) ∇ g(_t)) (  _t^T - η(  ^-1() ∇ g(_t))^T     ) 
       = (_t _t^T) - η_t (  ^-1(_t) ∇ g(_t)  )^T  - η^-1(_t) ∇ g(_t) _t^T + 
       η^2 (^-1(_t) ∇ g(_t)  )  · (^-1(_t) ∇ g(_t)  )^T

The block encoding of _t _t^T is apparently presented, more concretely, in quantum Newtons' method, it is from the previous t-1-th step. Previously, we have constructed the (ϵ-approximated ) block encoding of H^-1/Γ and of ∇ g(_t) _t^T. The block encoding of ∇1/p g(_t) _t^T naturally yields the block encoding of (∇1/p g(_t) _t^T)^T = 1/p_t ∇ g(_t)^T. Then we can use lemma <ref> to construct the block encoding of 
   
    1/p_t  ∇ g(_t)^T  (^-1(_t)/Γ)^T

We then can use lemma <ref> to insert the factor η, that we transform the block encoding of the above operator into 

    η/p_t  ∇ g(_t)^T  (^-1(_t)/Γ)^T

The unitary transpose of the block encoding of above operator is exactly  the block encoding of

    η/p^-1(_t)/Γ ∇ g(_t)  _t^T

We can use lemma <ref> to construct the block encoding of their product, e.g.,

    ( η/p^-1(_t)/Γ ∇ g(_t)  _t^T ) ·(  _t  ∇ g(_t)^T  (^-1(_t)/Γ)^T η/p) 
     
       = |_t|^2 η^2/p^21/Γ^2 (^-1(_t) ∇ g(_t)  )  · (^-1(_t) ∇ g(_t)  )^T

In order to remove the factor |_t|^2, we just need to note that since we have the block encoding of _t _t^T, we can use lemma <ref> to efficiently find its maximum eigenvalue, which is exactly |_t|^2. Then we can use the amplification method to remove such factor, resulting in further complexity 𝒪( 1/|_t|^2 ). 

Due to the term p^2 Γ^2 (which is known as roughly one over minimum eigenvalue of (_t) , we need to use lemma <ref>) to change the block encoding of _t _t^T into 1/p^2Γ^2_t _t^T,  and of (η/p) _t  ∇ g(_t)^T  (^-1(_t)/p Γ)^T into  η_t  ∇ g(_t)^T  (^-1(_t)/p Γ^2)^T (same thing for its transpose). Then we employ lemma <ref> to construct the block encoding of 

    1/4p^2 Γ^2((_t _t^T) - η_t (  ^-1(_t) ∇ g(_t)  )^T  - η^-1(_t) ∇ g(_t) _t^T +  η^2 (^-1(_t) ∇ g(_t)  )  · (^-1(_t) ∇ g(_t)  )^T   )

which is exactly the block encoding of _t+1_t+1^T/(4pΓ^2) where we remind that Γ^2 is ∼ one over minimum eigenvalue of (_t), or ||(_t)^-1||, which is known because we estimated it at first. This factor can be removed by amplification, resulting in further 𝒪(Γ^2) complexity. Therefore, the total running time of a single step of our quantum Newton's method is

    𝒪(  𝒪(γ^2p-2 2s p^2  ( p^2log(n) + plog^2.5(1/ϵ) )  κ  polylog(κ/ϵ) )  )


The improvement compared to <cit.> turns out to be substantial as the number of copies requirement for ^T is significantly reduced. At a certain step of Newton's method, the work in <cit.> requires 𝒪(p^5/ϵ^3) copies of the temporal solution to perform an update of the solution, with a total time complexity 𝒪(p^8 log(n) κ /ϵ^4  ). Meanwhile, the running time of our work depends polylogarithmically in terms of error tolerance, and that the number of “copies” of (block encoding of) ^T is p, which is a major advantage compared to <cit.>. In particular, the dependence on p reduces by a power-of-4. The most important thing is that our framework generalizes to polynomials of all kinds, not limited to homogeneous, even-degree ones. We finally emphasize that according to the analysis given in <cit.>, the factor γ from the above running time can be chosen to be arbitrarily small by choosing the parameter η properly, e.g., choosing η≤ 1/(2p ||||^-1) guarantees that γ≤ 4 (see Section 4, Result 4 of <cit.> for detailed derivation). 


An important aspect of our improved Newton's method is that the inversion of Hessian () is more accurate. We recall that the method <cit.> requires a reasonable lower bound on the eigenvalues of (). However, in this case, the Hessian depends on , and it changes per each iteration. Therefore, in order to apply the method in <cit.>, we need to estimate the smallest eigenvalue in magnitude, which can be done efficiently, e.g., logarithmically with respect to the size of the matrix, using the method in <cit.>. We note that this extra step has a running time smaller than the inversion itself, which implies that, asymptotically, it does not increase the overall complexity. In <cit.>, the authors assumed the inversion is executed on some well-conditioned subspace of , with the cutoff threshold chosen to be some constant Λ_H^-1. Where and how to obtain a reasonable value for Λ_H^-1 is unclear to us; as we mentioned, the spectrum of Hessian  is not a constant over different iterations. Therefore, we believe that the extra step of estimating such a threshold is an important improvement upon <cit.>.  

As another application of our convexity testing framework, in particular, one can see that our framework can be applied to enhance the performance of (improved) quantum Newton's method as well as gradient descent method. The reason is straightforward to see, as the gradient descent and quantum Newton's method are aimed to optimize a function by shifting toward the minima from some initialization. Therefore, by scanning the landscape, one can see if the given function is convex within such region. If the function is convex, then a minima is guaranteed, and hence, the algorithm is considered successful. 




§ CONCLUSION

In this work, we have investigated the potential of quantum computers in the context of functional analysis, specifically testing the convexity of a given objective function. The problem of testing convexity is converted to the problem of determining the non-negativity of a matrix, e.g., all eigenvalues are non-negative. By employing a useful algebraic property of homogeneous polynomials of even degrees, we wrote down an explicit form of the so-called Hessian and built upon it to generalize the Hessian to arbitrary polynomials. We then combine it with a powerful quantum singular value transformation framework, plus a quantum power method to construct a quantum algorithm that allows us to test the positive-semidefiniteness of such Hessian matrix at multiple points from a given domain. The procedure has running time polylogarithmically respective to the dimension n and linear with the number of sample points , which is superpolynomially faster than the best-known classical approach relative to the number of variables, meanwhile keeping the same complexity on the number of sample points. We also point out two examples where we envision a potential application of convexity testing, including studying the geometric structure of manifold, testing training landscape of variational quantum algorithms as well as optimization landscape for gradient descent/Newton's method for optimization. In particular, as a striking corollary of our result, we provide a major improvement upon the work of <cit.> in multiple aspects: running time (with respect to error tolerance), generality (arbitrary polynomial type) and the subtle detail regarding the inversion of the Hessian. As a whole, our work has added one more interesting example to the field of quantum computation, revealing that the area of functional analysis, and many more areas, is a rich avenue deserving of further research from computational aspects. An interesting open direction that we believe is worth looking at is that the model we work on in this case is somewhat explicit, e.g., the function is a polynomial of computable coefficients given via an oracle. In some traditional problems, the input function is typically given as a blackbox, such as Grover's search problem <cit.>, and we wish to reveal hidden properties. Suppose, instead, we are given a blackbox function that computes some analytical function; then, how do we extract the hidden properties, such as the convexity of the function in some domain? We note that there are two prior relevant works such as <cit.>, where the author considered the gradient estimation, and <cit.>, where the authors considered a sampling problem from a blackbox oracle that computes the value of given function within some domain. While there are seemingly overlaps, we do not see a direct solution to our case, thereby we leave the challenge for future investigation. 




§ ACKNOWLEDGEMENT

We thank Hiroki Sukeno and Shuyu Zhang for carefully reading and detailed feedback on the manuscript. This work was supported by the U.S. Department of Energy, Office of Science, Advanced
Scientific Computing Research under Award Number DE-SC-0012704. 
We also acknowledge the support of a Seed Grant from
Stony Brook University’s Office of the Vice President for Research and the Center for Distributed Quantum Processing.




unsrt



§ PRELIMINARIES


Here, we summarize the main recipes of our work. We keep their statements brief but precise for simplicity, with their proofs/ constructions referred to in their original works.

 <cit.>
 
Let A be some Hermitian matrix of size N × N whose matrix norm |A| < 1. Let a unitary U have the following form:

    U = [ A ·; · ·;   ].

Then U is said to be an exact block encoding of matrix A. Equivalently, we can write:

    U = |0⟩⟨0|⊗ A + ⋯,

where | 0⟩ refers to the ancilla system required for the block encoding purpose. In the case where the U has the form 

    U  =  |0⟩⟨0|⊗A + ⋯,

where || A - A || ≤ϵ (with ||.|| being the matrix norm), then U is said to be an ϵ-approximated block encoding of A.


The above definition has multiple simple corollaries. First, an arbitrary unitary U block encodes itself. Suppose A is block encoded by some matrix U, then A can be block encoded in a larger matrix by simply adding ancillas (which have dimension m). Note that _m ⊗ U contains A in the top-left corner, which is a block encoding of A again by definition. Further, it is almost trivial to block encode the identity matrix of any dimension. For instance, we consider σ_z ⊗_m (for any m), which contains _m in the top-left corner.



Let ρ = _A |Φ⟩⟨Φ|, where ρ∈ℍ_B, |Φ⟩∈ℍ_A ⊗ℍ_B. Given unitary U that generates |Φ⟩ from | 0⟩_A ⊗| 0⟩_B, then there exists an efficient procedure that constructs an exact unitary block encoding of ρ.


The proof of the above lemma is given in <cit.> (see their Lemma 45). 




    Given the unitary block encoding of two matrices A_1 and A_2, an efficient procedure exists that constructs a unitary block encoding of A_1 A_2.


The proof of the above lemma is also given in <cit.>.  




    Given the unitary block encoding {U_i}_i=1^m of multiple operators {M_i}_i=1^m (assumed to be exact encoding), then, there is a procedure that produces the unitary block encoding operator of ⊗_i=1^m M_i, which requires a single use of each {U_i}_i=1^m and 𝒪(1) SWAP gates. 

The above lemma is a result in <cit.>. 


    Given the oracle access to s-sparse matrix A of dimension n× n, then an ϵ-approximated unitary block encoding of A/s can be prepared with gate/time complexity 𝒪(log n + log^2.5(1/ϵ)).

This is also presented in <cit.>.  One can also find similar construction in Ref. <cit.>. 

    Given unitary block encoding of multiple operators {M_i}_i=1^m. Then, there is a procedure that produces a unitary block encoding operator of ∑_i=1^m ± M_i/m in complexity 𝒪(m).
    




    Given a block encoding of some matrix A (as in <ref>), then the block encoding of A/p, where p > 1, can be prepared with an extra 𝒪(1) cost.




§ PROOF OF LEMMA <REF>


Remind that we are trying to prove the following

   Given the block encoding of ⊕_i=1^ cc^T = _ cc^T and ⊕_i=1^_i_i^T, then it is able to constructs the block encoding of the diagonal matrix ℬ with entries ℬ_ij = β_i^2  δ_ij where β_i = _i^T c

From the block encoding of ⊕_i=1^ cc^T = _ cc^T as in the above lemma, plus the block encoding of ⊕_i=1^_i _i^T, we can use lemma <ref> to construct the block encoding of their products, e.g, we would obtain the block encoding of ⊕_i=1^β_i c _i^T, denoted as U_β. Given the definition of block encoding (<ref>), we have that:

    U_β| 0⟩_u |Φ⟩ = | 0⟩_u (⊕_i=1^β_i c _i^T) |Φ⟩ + |Garbage⟩

where |Garbage⟩ satisfies: | 0⟩⟨ 0|⊗·|Garbage⟩ = 0. Given that ⊕_i=1^β_i c _i^T = ∑_i=1^|i⟩⟨i|⊗β_i c _i^T and if we choose |Φ⟩ =  |j⟩⊗|c⟩, we have that 

    ⊕_i=1^β_i c _i^T ·|Φ⟩ = δ_ijβ_i β_j  c

Recall further that we also have a unitary U_C such that U_C| 0⟩ = |c⟩≡ c. Therefore, the state |Φ⟩ could be created from | 0⟩⊗|j⟩ by simply applying U_C to the first register, e.g, to obtain |c⟩⊗|j⟩ and use SWAP gate to swap them, e.g, we obtain |j⟩⊗|c⟩. Denote such process as U_s.

To match the same dimension as U_β, we simply add extra register (with corresponding dimension). We have that:

    ⟨ 0|_u ⟨ 0|⟨i|  (_u ⊗ U_C^†)   U_β (_u ⊗ U_s) | 0⟩_u | 0⟩|j⟩ = δ_ijβ_i β_j

which basically is a block encoding of a diagonal matrix that contains β_i^2. Using the result from <cit.> to remove the power factor 2, we are left with block encoding of δ_ijβ_i β_j. 




