TECHNICAL FIELD OF THE INVENTION
The invention relates to calculating an observable sampled out of a quantum state by using a non-quantum computer instead of a quantum computer.
BACKGROUND OF THE INVENTION
At present time, existing quantum computers are limited to handling a very low number of qubits, typically less than twenty-five qubits, and their operation is limited by short decoherence time and instruction faults. In this context, comparing algorithms to be implemented on future quantum computers needs developing new tools that can be implemented using non-quantum computers as available now, also called classical computers. Issues are to test, debug and analyse execution of an algorithm as it will be executed by a quantum computer. To this end, quantum computer operation is emulated using a classical computer. Then, for handling by the classical computer, a n-qubit vector is represented as a 2n-length vector, also denoted 2n-vector, which is comprised of 2n complex coordinates where n is a non-zero integer.
Known methods for emulating the operation of the quantum computer consist in having the classical computer calculating matrix and vector products as linear combinations of their complex coefficients and coordinates. But for such methods, increasing n-value by one multiplies by a factor of about two the size of the memory that is required for the classical computer, and also the computation time. Put another way, required memory and calculation time increase exponentially with the number n of the qubits in the quantum states.
Starting from this situation, one object of the present invention consists in providing a new method of simulating the operation of a quantum computer using a classical computer, which is less demanding in memory size and calculation time.
SUMMARY OF THE INVENTION
For meeting this object or others, a first aspect of the present invention proposes a method for calculating, using a non-quantum computer, a value <0|U†HU|0> of an observable H sampled out of a quantum state U|0>, when this quantum state U|0> results from a unitary operator U that corresponds to a quantum circuit and operates with an initial quantum state |0>. Executing the invention method with the non-quantum computer thus emulates the operation of the quantum circuit which corresponds to the unitary operator U, for calculating the observable H. This invention method comprises the following steps:
- /1/ expressing the unitary operator U as an ordered decomposition product of at least one Pauli rotation RPj(θj), multiplied by the identity operator I2n or by a Clifford operator C different from this identity operator I2n, each Pauli rotation equalling a sum of the identity operator I multiplied by cosine of half of a rotation angle θj of this Pauli rotation and a first Pauli operator Pj multiplied by sine of the half of the rotation angle θj and by the opposite of the imaginary unit;
- /2/ expressing the observable H as a linear combination of second Pauli operators Pj′, or if step /1/ involves the Clifford operator C that is different from the identity operator I2n: base-transforming the observable H using the Clifford operator C and expressing the base-transformed observable, namely C†HC, as a linear combination of second Pauli operators;
- /3/ expressing each first Pauli operator Pj as a first sum of transfer operators that each select one coordinate of a quantum state expressed in a calculation base and transfer the selected coordinate onto a selected base state, the selected coordinates and the selected base states, and also respective phases of the transfer operators in the first sum being determined separately for each first Pauli operator Pj1;
/4/ calculating a final state |ψfinal> by applying each Pauli rotation RPj(θj) as expressed in step /3/, in a chained manner from the initial quantum state |0> according to the ordered decomposition product of the unitary operator U obtained in step /1/;
- /5/ expressing each second Pauli operator Pj′ as a second sum of transfer operators that each select one coordinate of a quantum state expressed in the calculation base and transfer the selected coordinate onto a selected base state, the selected coordinates and the selected base states, and also respective phases of the transfer operators in the second sum being determined separately for each second Pauli operator Pj′;
- /6/ for each second Pauli operator Pj′ and using the second sum obtained in step /5/ for this second Pauli operator Pj′, calculating a contribution value <ψfinal|Pj′|ψfinal> as a result of right- and left-combination of the second Pauli operator Pj′ with the final state |ψfinal> calculated in step /4/; and
- /7/ inputting the contribution values <ψfinal|Pj′|ψfinal> calculated in step /6/ into the linear combination obtained in step /2/ for the observable H or base-transformed observable C†HC, so as to obtain as a calculation result the value <0|U†HU|0> of the observable H sampled out of the quantum state U|0>.
Thanks to expressing each Pauli rotation using a Pauli operator, and then expressing each Pauli operator that is used in an elementary form, the invention method is over-optimized for reducing the necessary calculation resource and calculation time. In particular, steps /3/ and /5/ transform each Pauli operator into a sum of coordinate-switching operations, corresponding to the so-called transfer operators, which can be executed in a very simple and cheap manner by the non-quantum computer.
In preferred implementations of the invention, for each first Pauli operator Pj or second Pauli operator Pj′ involved in steps /3/ and /5/, the transfer operators and their respective phases in the first or second sum may be obtained from a x-vector and a z-vector both of length n and determined for this first Pauli operator Pj or second Pauli operator Pj′, where n is the qubit number of the quantum state and an integer index q is numbering the qubits of the quantum state according to the calculation base. The x-vector and z-vector are determined in the following way:
- the qth coordinate of the x-vector equals unity if the first or second Pauli operator acts as the X-Pauli 2×2 operator or Y-Pauli 2×2 operator onto the qth qubit of any quantum state expressed in the calculation base, otherwise the qth coordinate of the x-vector equals zero, and
- the qth coordinate of the z-vector equals unity if the first or second Pauli operator acts as the Z-Pauli 2×2 operator or Y-Pauli 2×2 operator onto the qth qubit of any quantum state expressed in the calculation base, otherwise the qth coordinate of the z-vector equals zero.
Generally for the invention, the observable H may be useful in a variational algorithm, thus encompassing a broad spectrum of applications. For such applications, the invention provides a reduction in the simulation time.
Generally again, the observable H may be useful in applications pertaining to quantum chemistry, combinatorial optimization and machine learning.
For example, the invention method may be implemented for executing algorithms such as Variation Quantum Eigensolver in quantum chemistry. It may be implemented in particular for estimating the ground energy of complex physical or chemical systems or materials, such as molecules, and thus lead to chemical and electromagnetic behaviours of these systems. Indeed, algorithms used for such physical or chemical applications often rely on execution of a sequence of Pauli rotations combined with sampling of an observable which is specified in the Pauli basis. The observable is then the Hamiltonian of the physical or chemical system considered. Such applications include designing molecules that are efficient as medicines, dyeing agents, liquid crystals, etc.
The most generic approach in combinatorial optimization, in particular for hard optimization problems, is that of Quantum Approximate Optimisation Algorithms, known as QAOA. It is also based on Pauli rotations followed by sampling of a Pauli-sparse Hamiltonian. The QAOA algorithms can thus be emulated efficiently using the invention method.
Finally, the invention method may be implemented to simulate simple noise models, in particular for noise as existing when executing a quantum program on an imperfect hardware. Indeed, depolarizing noise may be simulated by applying Pauli rotations.
In addition, a second aspect of the invention proposes a computer program product that comprises instructions for executing the method of the first invention aspect when this computer program product is implemented by a non-quantum computer. Thus, in some embodiments, non-transitory computer-readable storage device(s) may include instructions that, when executed by one or more non-quantum processors operably coupled to the non-transitory computer-readable storage device(s), the one or more non-quantum processors, individually or collectively, perform a method as disclosed herein.
These and other features of the invention will be now described with reference to the appended figures, which relate to preferred but not-limiting implementations of the invention.
BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 shows a classical computer implementing the invention method;
FIG. 2 is a block diagram that shows steps of an implementation of the invention method;
FIG. 3 displays a possible first algorithm that may be useful for executing the invention method, when applying a Pauli rotation to a quantum state;
FIG. 4 displays a possible second algorithm that may also be useful for executing the invention method, when sampling a Pauli operator out of a quantum state;
FIG. 5a is a diagram that compares calculation times between the invention method and a prior-art method, as a function of a qubit number; and
FIG. 5b is another diagram that compares calculation times between the invention method and the prior-art method of FIG. 5a, as a function of a number of involved Pauli rotations.
DETAILED DESCRIPTION OF THE INVENTION
In the present description, classical computer denotes a non-quantum computer, as opposed to an actual quantum computer.
The object of the invention is calculating the value <0|U†HU|0> by using a classical computer for emulating the operation of a quantum computer that would be used for calculating the same value <0|U†HU|0>. The following notations are used throughout hereafter:
- |0>: initial quantum state to be used for the invention, which may be one of the quantum base states used for calculation or any other quantum state which is preferably a stabilizer state;
- n: number of qubits combined through tensor product in each quantum state considered herein. Therefore, each quantum state is described classically by an ordered series of 2n complex coordinates and called n-qubit;
- H: observable, defined as a 2n×2n operator, for example a Hamiltonian in an application of the invention to quantum chemistry;
- U: unitary operator to be applied to the initial quantum state |0> before assessing the value of the observable H. The value to be calculated is therefore <0|U†HU|0>, where U† denotes the complex conjugate transposition of the operator U. The unitary operator U corresponds to a quantum circuit to be emulated using the classical computer;
- C: any Clifford operator involved when expressing the unitary operator U as a product of Pauli rotations. It is assumed that the Clifford group in well-known, composed of so-called Clifford operators. In the particular case of a Clifford operator C, C† is equal to the reverse operator of C, so that C†HC is equivalent to converting the observable H through a base transformation;
- P: any Pauli operator acting onto quantum states comprised on n qubits combined through tensor product. In a known manner, each Pauli operator P acts as the identity 2×2 operator I2 or as any of the X-, Y- and Z-Pauli 2×2 operators for the qth qubit of any n-qubit, q being an integer index comprised between 1 and n, values 1 and n included;
- RP(θ): the Pauli rotation acting onto quantum states comprised on n qubits combined through tensor product, and having the Pauli operator P as an axis and rotation angle θ;
- or I2n: identity operator sized for n qubits combined through tensor product, I2 then denoting the identity 2×2 operator;
- |ψ>: any quantum state comprised of n qubits combined through tensor product, and thus expressed as an ordered series of 2n complex coordinates in the calculation base;
- |k>: base state comprised of n base qubits combined through tensor product; and
- i: the imaginary unit of complex numbers.
In a well-known manner, any Pauli rotation RP(θ) can be written in the following manner (Eq. 1):
R
P(θ)=cos(θ/2)−i sin(θ/2)P
The minus sign before the imaginary unit i in this equation Eq. 1 is provided by convention and may be converted into plus sign if opposite convention is used.
Again in a well-known manner, any Pauli operator P can be written in the following manner (Eq. 2):
P=i
−v
z
·v
x
Z
v
z
X
v
x
- where X denotes the X-Pauli 2×2 operator, Z denotes the Z-Pauli 2×2 operator, vx and vz are n-length vectors that depend on the Pauli operator P, and that are constructed in the following manner:
- for vector vz, also called z-vector in the general part of the present description: the qth coordinate of vector vz equals unity if the Pauli operator P acts as the Z-Pauli 2×2 operator or Y-Pauli 2×2 operator onto the qth qubit of any quantum state |ψ> expressed in the calculation base, otherwise the qth coordinate of vz equals zero corresponding to the Pauli operator P acting as I2 onto the qth qubit; and
- for vector vx, also called x-vector in the general part of the present description: the qth coordinate of vector vx equals unity if the Pauli operator P acts as the X-Pauli 2×2 operator or Y-Pauli 2×2 operator onto the qth qubit of any quantum state |ψ> expressed in the calculation base, otherwise the qth coordinate of vx equals zero corresponding to the Pauli operator P acting as I2 onto the qth qubit.
In the above equation Eq. 2, Zvz denotes the operator acting on the quantum state |ψ> that is obtained by tensor product of the Z-Pauli 2×2 and I2 operators, with using the Z-Pauli 2×2 for the qth qubit if the qth coordinate of vector vz equals one, otherwise using the I2 operator for the qth qubit. Similarly, Xvx denotes the operator acting on the quantum state |ψ> that is obtained by tensor product of the X-Pauli 2×2 and I2 operators, with using the X-Pauli 2×2 for the qth qubit if the qth coordinate of vector vx equals one, otherwise using the I2 operator for the qth qubit. The phase factor i−vz·vx comes from the identity Y=−iZX, and vz·vx denotes the dot product of both vectors vz and vx.
Then, one can show that the expression of Pauli operator P can be transformed from the equation Eq. 2 into the following equation (Eq. 3):
In equation Eq. 3, k denotes the n-length vector that is comprised of bit values as coordinates each equal to the corresponding qubit figure in |k> expressed in the tensor-product calculation base, ⊕ is the well-known XOR operator acting on bits, and k·vz denotes the dot product of both vectors k and vz.
Intermediate development lines for obtaining equation Eq. 3 from equation Eq. 2 are:
When applied to the quantum state |ψ>, equation Eq. 3 becomes (Eq. 3′):
The Pauli-transformed state P|ψ> can be calculated in this way by implementing a number of complex operations that is proportional to 2n. |k><k⊕vx| has been called transfer operator in the general part of this description. It selects one coordinate of the quantum state |ψ> and transfers it onto the base state |k>.
In FIG. 1, reference number 1 denotes a classical computer used for implementing the invention method. References 1a and 1b denote inputs thereof for entering the unitary operator U and the observable H, respectively, and reference 1c denotes the output that delivers the value <0|U†HU|0>.
Referring now to FIG. 2, the invention method starts with step /1/ that consists in performing a product-decomposition of the unitary operator U using Pauli rotations. This product-decomposition involves an ordered series of Pauli rotations (RPj(θj)), each having the Pauli operator Pj as its axis and rotation angle θj. The angles θj are non-Clifford, i.e. each θj is not a multiple of π/2. Thus, the unitary operator U is expressed as set in the following equation (Eq. 4):
U=CR
P
m
(θm) . . . RP1(θ1)
In equation Eq. 4, C is a Clifford operator which depends on the unitary operator U. In particular, the Clifford operator C may be , but may also be any other Clifford operator. As it is well-known, identity operator I2n is a particular Clifford operator. The product of the Pauli rotations does not commutate although it is written in contracted form within box /1/ in FIG. 2. In addition, each one of the Pauli rotations RPj(θj) can be expressed in accordance with above equation Eq. 1, using each time Pauli operator Pj and rotation angle θj as P and θ, respectively. The Pauli operators Pj have been called first Pauli operators in the general part of the present description. The calculation time that is necessary for obtaining the decomposition product of the unitary operator U as displayed in the above equation Eq. 4 is proportional to n2m, where m is the number of Pauli rotations in the decomposition product and n is again the number of qubits that are combined through tensor product.
Step /2/ consists in performing a linear decomposition of the observable C†HC using Pauli operators, with the Clifford operator C as determined in step /1/. When the Clifford operator C is different from , C†HC has been called base-transformed observable. In this way, the observable or base-transformed observable is expressed as a linear combination of Pauli operators Pj′ in accordance with the following equation (Eq. 5):
H or C†HC=Σj′aj′Pj′
where aj′ are complex coefficients of the linear combination. The Pauli operators Pj′ have been called second Pauli operators in the general part of the present description.
Step /3/ consists in applying the above equation Eq. 3 to each one of the first Pauli operators Pj. The sums that appear in this equation Eq. 3 when applied to the first Pauli operators Pj have been called first sums in the general part of the present description.
Then, the expression of each first Pauli operator Pj as resulting from step /3/can be inserted in the above equation Eq. 1. The so-obtained expression for each Pauli rotation RPj(θj) is then used for applying this latter to any quantum state |ψ>. This may be achieved by implementing Algorithm 1 as displayed in FIG. 3 in the form of a pseudo-code. The inputs of Algorithm 1 are thus the quantum state |ψ> expressed by its 2n complex coordinates, and the Pauli rotation RP(θ) determined by its axis P and its rotation angle θ. The rotated quantum state RP(θ)|ψ> is then the output |ψ′> of Algorithm 1 also expressed by its 2n complex coordinates <k|ψ′>. The successive lines of the pseudo-code of Algorithm 1 as displayed in FIG. 3 have the following meanings:
- line 1: initializing of the output quantum state |ψ′>
- line 2: determination of the n-length vector vz by implementing a subroutine to the Pauli operator P, for determining its components that equal Z- or Y-Pauli 2×2 operator. Such ComputeMaskZ subroutine is well-known and basic to implement
- line 3: determination of the n-length vector vx by implementing another subroutine to the Pauli operator P, for determining its components that equal X- or Y-Pauli 2×2 operator. Such ComputeMaskX subroutine is also well-known and basic to implement
- line 4: calculation of the phase factor p equalling to the phase factor i−vz·vx as appearing in the above equation Eq. 2. This is achieved by implementing a third subroutine ComputeYPhase that uses in a basic manner the results of both subroutines ComputeMaskZ and ComputeMaskX applied to the Pauli operator P
- lines 5 to 7: iteration loop for calculating each coordinate of the rotated quantum state RP(θ)|ψ>=|ψ′>, by using the combined equations Eq. 1 and Eq. 3′.
Algorithm 1 is simple and does not require any synchronisation between sub-results being made available.
In step /4/, Algorithm 1 is implemented repetitively in a chained manner so as to calculate the result of the product RPm(θm) . . . RP2(θ2)RP1(θ1) as appearing in equation Eq. 4, when applied to the initial quantum state |0>: it is first applied with P=P1, θ=θ1 and |ψ>=|0> for obtaining RP1(θ1)|0>, then applied with P=P2, θ=θ2 and |ψ>=RP1(θ1)|0> for obtaining RP2(θ2)RP1(θ1)|0>, and then with P=P3, θ=θ3 and |ψ>=RP2(θ2)RP1(θ1)|0> for obtaining RP3(θ3)RP2(θ2)RP1(θ1)|0>, etc until RPm(θm) . . . RP2(θ2)RP1(θ1)|0> is obtained. The resulting quantum state has been called final state and denoted |ψfinal>.
Step /5/ consists in applying the above equation Eq. 3 to each one of the second Pauli operators Pj′. The sums that appear in this equation Eq. 3 when applied to the second Pauli operators Pj′ have been called second sums in the general part of the present description.
Step /6/ consists in calculating a contribution value for the result of <0|U†HU|0>, that originates from each second Pauli operator Pj′ as involved in the above equation Eq. 5. This contribution value is <ψfinal|Pj′|ψfinal>, and the expression obtained in step /5/ for the Pauli operator Pj′ is used to this end. This corresponds to applying the following equation (Eq. 6) to the Pauli operator Pj′:
in which the meanings already explained apply again. The contribution value <ψfinal|Pj′|ψfinal> can be calculated in this way with a computation time that is proportional to 2n.
This calculation may be achieved by implementing Algorithm 2 as displayed in FIG. 4 in the form of another pseudo-code. The inputs of Algorithm 2 are thus the final state |ψfinal> expressed by its 2n complex coordinates as resulting from step /4/, and the Pauli operator P successively equals all second Pauli operators Pj′. For each execution, the output of Algorithm 2 is the contribution value <ψfinal|Pj′|ψfinal>. The successive lines of the pseudo-code of Algorithm 2 as displayed in FIG. 4 have the following meanings:
- line 1: initializing of an addition result r
- lines 2 to 4: determination of the n-length vectors vz and vx and also the phase factor p in a manner similar to that used in Algorithm 1, by applying the subroutines ComputeMaskZ, ComputeMaskX and ComputeYPhase to the Pauli operator Pj′
- lines 5 to 7: iteration loop for calculating after one another each term of the sum that appears in the above equation Eq. 6. Once calculated, each term is added to the previous value of the addition result r
- line 8: multiplication by the phase factor p so as to obtain the contribution value <ψfinal|Pj′|ψfinal>.
Algorithm 2 is simple again and does not require either any synchronisation between sub-results.
Algorithm 2 is applied to each second Pauli operator Pj′ that is involved in the above equation Eq. 5.
Finally, step /7/ consists in combining the contribution values as resulting from step /6/ in accordance with the above equation Eq. 5. The combination result is the desired value of <0|U†HU|0>.
The diagram of FIG. 5a displays the calculation times necessary for obtaining the value of <0|U†HU|0> for a unitary operator U which involves fifty Pauli rotations in its decomposition product, when using a classical computer which includes twenty threads. Vertical axis indicates the calculation time values expressed in seconds (s), in logarithmic scale. Horizontal axis indicates the values of the qubit number n. Both curves correspond to the following calculations:
- curve labelled “linalg”: the calculation is carried out using standard matrix and vector products as implemented before the present invention; and
- curve labelled “paulisim”: the calculation is carried out using the present invention, including implementations of Algorithm 1 and Algorithm 2.
All other calculation conditions and parameters are identical between both curves. It is clear that the invention provides significant reduction in the calculation time.
The diagram of FIG. 5b corresponds again to implementing standard matrix and vector products (curve labelled “linalg”) and implementing the present invention including Algorithm 1 and Algorithm 2 (curve labelled “paulisim”), but with fixed qubit number n equal to twenty and varying Pauli rotation number m as appearing in the decomposition product of the unitary operator U. All other calculation conditions and parameters are identical again between both curves. Vertical axis indicates again the calculation time values expressed in seconds (s), but in linear scale. Horizontal axis indicates the values of the Pauli rotation number m. Reduction in the calculation time is also clear in this way.