Ansätze for noisy variational quantum Eigensolvers

Universidade do Minho Escola de Engenharia Mafalda Francisco Ramôa da Costa Alves Ansätze for Noisy Variational Quantum Eigensolvers November, 2021 Universidade do Minho Escola de Engenharia Mafalda Francisco Ramôa da Costa Alves Ansätze for Noisy Variational Quantum Eigensolvers Master Thesis Engineering Physics Work developed under the supervision of: Ernesto Galvão Mikhail Vasilevskiy November, 2021 COPYRIGHT AND TERMS OF USE OF THIS WORK BY A THIRD PARTY This is academic work that can be used by third parties as long as internationally accepted rules and good practices regarding copyright and related rights are respected. Accordingly, this work may be used under the license provided below. If the user needs permission to make use of the work under conditions not provided for in the indicated licensing, they should contact the author through the RepositoriUM of Universidade do Minho. License granted to the users of this work Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International CC BY-NC-SA 4.0 https://creativecommons.org/licenses/by-nc-sa/4.0/deed.en iv Acknowledgements I would like to offer my special thanks to Raffaele Santagati for suggesting this project, and for his guidance and support throughout the course of it. I thank also my co-supervisor Ernesto Galvão, for his assistance at every stage of the project, as well as for his invaluable advice; and Mikhail Vasilevskiy, for his support. I thank my family and friends for their encouragement. In particular, I’m deeply thankful to my sister, colleague and friend Alexandra. This year would not have been the same without her continuous support both inside and outside the academic sphere. I want to express my gratitude also to all the teachers and professors that, throughout the years, have taught me so many different things in so many different ways. The genuine concern for students and dedication for teaching I have encountered throughout my life have certainly impacted it along the years. Additionally, I wish to thank the Programme New Talents in Quantum Technologies of the Gulbenkian Foundation (Portugal) for the financial support, and for the opportunity of being involved in this experience. The network of students under the support of the Gulbenkian Foundation is unmatched in the enthusiasm for learning and growing. Finally, I wish to acknowledge the use of IBMQ’s systems [42], Google Colaboratory’s cloud service [30], and the NOVAThesis L A T EX template [56]. v STATEMENT OF INTEGRITY I hereby declare having conducted this academic work with integrity. I confirm that I have not used plagiarism or any form of undue use of information or falsification of results along the process leading to its elaboration. I further declare that I have fully acknowledged the Code of Ethical Conduct of the Universidade do Minho. vi Abstract Ansätze for Noisy Variational Quantum Eigensolvers Simulating quantum mechanical systems is one of the main applications envisioned for quantum computers. In contrast with the first algorithms created for this purpose, that were devised to be implemented in a fault-tolerant quantum computer, the Variational Quantum Eigensolver (VQE) aims to adjust to the constraints of Noisy Intermediate-Scale Quantum (NISQ) devices. There is hope that the class of Variational Quantum Algorithms (VQAs), to which VQE belongs, will be the first to achieve quantum advantage. The choice of ansatz can dictate the success (or lack thereof) of a VQA: too deep ansätze can hinder near-term viability, or lead to trainability issues that render the algorithm inefficient. In this context, this dissertation aimed to analyse different ansätze for quantum chemistry, examining their noise-resilience and viability in state-of-the-art quantum computers. In particular, dynamic ansätze were explored, and their performance compared against predetermined ansätze, with a focus on susceptibility to noise. Multiple variants of VQE, namely Unitary Coupled Cluster Singles and Doubles (UCCSD)-VQE (predetermined) and Adaptive Derivative-Assembled Pseudo-Trotter (ADAPT)-VQE (dynamic), were implemented both in simulators and cloud quantum computers. Using noise models, the impact of several noise sources on convergence was assessed. Additionally, the importance of the operator pool in ADAPT-VQE was analysed, and strategies to manipulate the ansatz beyond the ADAPT-VQE algorithm were explored. Several conclusions could be drawn from this work. Adapting the ansatz to the problem and system was concluded to be fundamental in avoiding trainability issues, decreasing the circuit depth required for a given accuracy, and improving noise-resilience (against circuit depth dependent and independent sources alike). Dynamic ansätze were shown to be capable of enduring significantly larger error rates than predetermined alternatives, and were thus proved to be better suited for NISQ devices. For 𝐻2, as far as ground state energy calculations are concerned, ADAPT-VQE was shown to tolerate a 20 times lower shot count, 150 times larger error rates in state preparations and measurements, and 850 times lower coherence times than UCCSD-VQE. The difference is expected to increase with the size of the system. Additionally, it was observed that there is still a margin for improving upon ADAPT-VQE. Further manipulation of the ansatz was shown to be capable of producing yet shallower circuits for the same accuracy. Using an idea previously proposed in the literature, a more conservative selection criterion was tested. Additionally, removing operators on the fly based on available data was attempted as a new possibility. Both approaches were shown to be capable of improving upon the ADAPT-VQE ansatz, resulting in an up to 35-fold decrease in the error for a similar circuit depth within the first 10 iterations of ADAPT-VQE. Keywords: adaptive ansätze, quantum chemistry, quantum computing, variational quantum algorithms, variational quantum eigensolver vii Resumo Ansätze para Eigensolvers Variacionais Quânticos Ruidosos Simular sistemas quânticos é uma das mais relevantes potenciais aplicações dos computadores quânticos. Os primeiros algoritmos criados para tal visam implementação em computadores quânticos tolerantes a falhas; em contraste, o Eigensolver Variacional Quântico (VQE) procura adaptar-se às limitações dos dispositivos quânticos ruidosos de escala intermédia (NISQ). Há esperança que a classe de algoritmos variacionais quânticos (VQAs), à qual o VQE pertence, seja a primeira a alcançar a vantagem quântica. O ansatz pode ditar o sucesso (ou falta de) de um VQA: a sua profundidade pode impedir a viabilidade a curto prazo, ou conduzir a problemas de treinabilidade que levam à ineficiência. Neste contexto, esta dissertação pretendeu analisar diferentes ansätze para química quântica, examinando a sua resiliência ao ruído e viabilidade em dispositivos contemporâneos. Em particular, exploraram-se ansätze dinâmicos, cujo desempenho se comparou ao de ansätze predeterminados com foco na suscetibilidade ao ruído. Diferentes variantes do VQE, nomeadamente UCCSD-VQE (predeterminado) e ADAPT-VQE (dinâmico), foram implementadas em simuladores e em computadores quânticos. Com recurso a modelos de ruído, o impacto de diversas fontes de ruído na convergência foi analisado. Inspecionou-se a importância do conjunto de operadores no ADAPT-VQE, e exploraram-se estratégias adicionais para manipular o ansatz . Várias conclusões puderam ser retiradas deste projeto. Concluiu-se que adaptar o ansatz ao problema e ao sistema é fundamental para evitar problemas de treinabilidade, diminuir a profundidade dos circuitos, e melhorar a resiliência ao ruído (seja a fonte dependente da profundidade do circuito ou não). Mostrouse que ansätze dinâmicos toleram taxas de erro superiores aos predeterminados, sendo portanto mais adequados para dispositivos NISQ. Quanto aos cálculos da energia do estado fundamental, mostrou-se que, para a molécula 𝐻2, o ADAPT-VQE tolera um número de repetições do circuito 20 vezes inferior, taxas de erro na preparação e medição de estados 150 vezes superiores, e tempos de coerência 850 vezes inferiores do que o UCCSD-VQE. É expectável que a diferença aumente com o tamanho do sistema. Adicionalmente, observou-se que há ainda uma margem para melhorar além do ADAPT-VQE. Mostrouse que outras formas de manipulação do ansatz são capazes de produzir circuitos ainda menos profundos para uma mesma precisão. Usando uma ideia previamente sugerida na literatura, tentou-se usar um critério de seleção mais conservador. Alternativamente, remover operadores com base em dados disponíveis ao longo da execução foi testado como uma nova possibilidade. Mostrou-se que ambas as estratégias são capazes de superar o ansatz criado seguindo o protocolo do ADAPT-VQE, resultando num erro até 35 vezes menor para profundidades de circuito comparáveis dentro de 10 iterações do ADAPT-VQE. Palavras-chave: algoritmos variacionais quânticos, ansätze adaptativos, computação quântica, eigensolver variacional quântico, química quântica viii Contents List of Figures xii Acronyms xx 1 Introduction 1 1.1 Context and Motivation .............................. 1 1.2 State of the Art .................................. 3 1.3 Objectives .................................... 7 1.4 Outline of the Document ............................. 8 2 Theoretical Framework 10 2.1 Quantum Computing ............................... 10 2.1.1 Basic Concepts ............................. 10 2.1.2 Measuring Pauli Strings ......................... 13 2.1.3 Quantum Simulation and Trotterization .................. 16 2.1.4 Current Limitations ............................ 19 2.2 Quantum Chemistry ............................... 21 2.2.1 The Electronic Problem .......................... 21 2.2.2 Finding the Ground State: Classical Approaches .............. 25 2.3 Variational Quantum Algorithms .......................... 33 2.3.1 Motivation ................................ 33 2.3.2 Outline .................................. 34 2.3.3 Variational Quantum Eigensolver ..................... 39 3 Static Ansätze for VQE 48 3.1 Problem Agnostic Ansätze ............................. 48 3.1.1 Application to Helium Hydride ...................... 48 3.2 Problem Tailored Ansätze ............................. 61 3.2.1 Application to Molecular Hydrogen .................... 61 4 Problem Tailored, Dynamically Created Ansätze: ADAPT-VQE 72 4.1 The ADAPT-VQE Algorithm ............................. 72 4.1.1 Fermionic-ADAPT-VQE .......................... 72 ix LIST OF FIGURES 27 Plots showing the impact of different types of noise on the performance of ADAPT-VQE for 𝐻2 at an interatomic distance of 0.74Å. For thermal relaxation, a simple model assuming average gate times and a constant ratio between the coherence times T1 and T2 was used. Figure 27a shows the evolution of the error in the energy as a function of these times. For SPAM noise, the ratio between the probability of measuring 0 after preparing 1 and the probability of measuring 1 after preparing 0 was fixed (and greater than one, for the model to be realistic). Figure 27b shows the evolution of the error as a function of these probabilities. Figure 27c illustrates the impact of sampling noise by plotting the error as a function of the number of shots. Twenty runs were used for each point; the median was used as the estimator (as it filters out aberrant runs), and the error bars show the interquartile ranges. All noise models were created using Qiskit [3]. The chosen optimizer was COBYLA. ........... 88 28 Evolution of the optimization of the first iteration of ADAPT-VQE (figure 28a) and of the UCCSDVQE algorithm, for 𝐻2at an interatomic distance of 0.74Å. IBMQ’s Belem processor was used in both cases. ..................................... 89 29 Comparison of the effect of different types of noise on ADAPT-VQE and UCCSD-VQE for 𝐻2 at an interatomic distance of 0.74Å. The same noise models described in figures 21 and 27 were used. ...................................... 91 30 Plots comparing the costs associated with UCCSD-VQE and ADAPT-VQE for a similar final error, in a setting with no noise other than sampling noise. For each interatomic distance, 10 runs of the algorithms were performed. The median of the errors is plotted in figure 30a. Figures 30b, 30c, and 30d present respectively the average number of energy evaluations, measured observables, and total shots per run. The number of shots per Pauli string was set to 213 and 218 for ADAPT-VQE and UCCSD-VQE respectively. These values were chosen so that the error was roughly matched between the two algorithms (the average error was 0.00038 a.u. and 0.00039 a.u. by the same order). .................. 93 31 Impact of the Jordan-Wigner string on the convergence of the qubit-ADAPT-VQE algorithm. The graphs compare the evolution of the error in the energy along the 20 first iterations of the algorithm, with and without the Z-strings in the pool operators. The molecule in study is the 𝐻4molecule at an interatomic distance of 1.45Å. ................. 95 32 Evolution of the error in the ADAPT-VQE energy along the evolution of the algorithm for the Spin-Complemented Generalized Singles and Doubles (SCGSD) pool, the SGSD pool, and the qubit pool. The molecule in study is the 𝐻4molecule. ............... 96 33 Circuit for implementing the operator in 5.8 using the ladder-of-CNOTs method explained in subsection 2.1.3. The circuit was drawn in the IBM Quantum Composer [41]. . . . . . . 107 34 Circuit for implementing the operator in 5.8 using the more efficient method introduced in [98]. The circuit was drawn in the IBM Quantum Composer [41]. ............ 107 xvi LIST OF FIGURES 35 Evolution of the 20 first iterations of the ADAPT-VQE algorithm for three different pools: the XXYX One Pool, the XXYX Two Pool, and the original Qubit Pool (without the Jordan-Wigner string). The upper plot shows the evolution of the error in the energy. The lower plot shows the evolution of the total probability amplitude of computational basis states representing Slater determinants with an altered particle number or 𝑍spin projection. The molecule in consideration is the 𝐻4molecule at an interatomic distance of 1.45Å. ......... 108 36 Evolution of the 20 first iterations of the ADAPT-VQE algorithm for different pools. The upper plot shows the evolution of the error in the energy. The lower plot shows the number of Slater determinants in the superposition state with altered particle number, 𝑆𝑧, or neither. In 36c and 36d, altered 𝑆𝑍and particle number correspond to the same colour because the curves coincide. When the green curve is not visible, it lies under the orange one. The molecule in consideration is 𝐻4at an interatomic distance of 1.45Å. ................ 110 37 Comparison of the convergence of the ADAPT-VQE algorithm using several different One Pools (figure 37a) and Two Pools (figure 37b) against the full Qubit Pool (without JordanWigner strings). In figure 37b, a curve corresponding to a One Pool was also plotted for reference. The error in the energy is plotted as a function of the iteration number. The molecule in consideration is 𝐻4at an interatomic distance of 1.45Å. .......... 112 38 Comparison of the convergence of the ADAPT-VQE algorithm using the Generalized Singles and Doubles (GSD) Pool and the Eight Pool. The error in the energy is plotted as a function of the iteration number. The molecule in consideration is 𝐻4at an interatomic distance of 1.45Å. ........................................ 113 39 Impact of the choice of pool (between One, Two, Four, and Eight Pools) in the ADAPT-VQE algorithm, for two different molecules. The error in the energy is plotted against the iteration number. In figure 39b, the XXYX Four Pool and the Eight Pool correspond to the same color because the curves coincide. These curves also overlap in some iterations in figure 39a. 114 40 Norm of the selected operator and change in energy in the first ten iterations of ADAPT-VQE. The molecule in consideration is the 𝐿𝑖𝐻 molecule at an interatomic distance of 1.45Å. The SGSD pool was used. ................................. 117 41 Absolute value of the performance ratio in the first ten iterations of ADAPT-VQE, for the same pool and molecule as figure 40. The performance ratio at a given iteration corresponds to the operator added then. ............................... 119 42 Relevant data structures at iteration 3 of the ADAPT-VQE algorithm, with an additional term removal feature based on blocking and unblocking operators. The molecule in consideration is 𝐿𝑖𝐻, and the pool the SGSD pool. The units of energy, omitted to avoid overcrowding of the figures, are Hartree throughout the whole chapter. Because the variational parameters are dimensionless, the gradient shares the energy units. ................ 122 43 Relevant data structures at iteration 4. ......................... 123 xvii LIST OF FIGURES 44 Relevant data structures at iteration 7. ......................... 124 45 Relevant data structures at iteration 8. ......................... 125 46 Relevant data structures at iteration 9. ......................... 125 47 Relevant data structures at iteration 11. ........................ 126 48 Relevant data structures at iteration 12. ........................ 126 49 Relevant data structures at iteration 14. ........................ 127 50 Relevant data structures at iteration 3 of the ADAPT-VQE algorithm, with an additional term removal feature based on a performance penalty. The molecule in consideration is 𝐿𝑖𝐻, and the pool the SGSD pool. ................................ 131 51 Relevant data structures at iteration 4. ......................... 132 52 Evolution of the relevant data, iteration 5 through 13. On top, we have data on the evolution of the algorithm and on the selected operators. On the bottom, we have data specifically regarding operator 25, that was removed at iteration 4. ................. 132 53 Relevant data structures at the end of iteration 11. ................... 135 54 Plot comparing the original ADAPT-VQE and the version with term removal, for LiH at an interatomic distance of 1.45Å and using the SGSD pool. The error in the energy (top) and the number of operators in the ansatz (bottom) are plotted against the iteration number. The shaded blue area on the upper plot marks the error values within chemical accuracy (less than 1kcal/mol). ................................... 137 55 Plots comparing the energy and energy error in the original ADAPT-VQE and the version with term removal, for LiH at an interatomic distance of 1.45Å and using the SGSD pool. The shaded blue areas on the plots mark the values within chemical accuracy both for the energy (less than 1kcal/mol away from the FCI energy) and the error (less than 1kcal/mol). 139 56 Plot comparing the original ADAPT-VQE and the version with term removal, for OH−at an interatomic distance of 1.45Å and using the SGSD pool. The error in the energy (top) and the number of operators in the ansatz (bottom) are plotted against the iteration number. The shaded blue area on the upper plot marks the error values within chemical accuracy (less than 1kcal/mol). ................................... 140 57 Plots comparing the energy and energy error in the original ADAPT-VQE and the version with term removal, for OH−at an interatomic distance of 1.45Å and using the SGSD pool. The shaded blue areas on the plots mark the values within chemical accuracy both for the energy (less than 1kcal/mol away from the FCI energy) and the error (less than 1kcal/mol). . . 141 58 Plot comparing the original ADAPT-VQE and the version with term removal, for OH−at an interatomic distance of 1.45Å and using the Eight Pool. The error in the energy (top) and the number of operators in the ansatz (bottom) are plotted against the iteration number. The shaded blue area on the upper plot marks the error values within chemical accuracy (less than 1kcal/mol). ................................... 142 xviii LIST OF FIGURES 59 Plots comparing the energy and energy error in the original ADAPT-VQE and the version with term removal, for OH−at an interatomic distance of 1.45Å and using the Eight Pool. The shaded blue areas on the plots mark the values within chemical accuracy both for the energy (less than 1kcal/mol away from the FCI energy) and the error (less than 1kcal/mol). . . 143 60 Plot showing the minimum error reached for each ansatz size, for multiple ansatz growth strategies. The molecule in case is LiH at an interatomic distance of 1.45Å. The SGSD pool was used. In the points where the green curve, corresponding to conservative ansatz growth, is not visible, it lies directly under the term removal curve. ............... 144 61 Plot showing the minimum error reached for each ansatz size, for multiple ansatz growth strategies. The molecule in case is 𝑂𝐻−at an interatomic distance of 1.45Å. The Eight Pool was used. The term removal and conservative ansatz growth strategies correspond to the same color because the curves coincide. ....................... 146 62 Plot showing the minimum error reached for each ansatz size, for multiple ansatz growth strategies. The molecule in case is 𝐻4at an interatomic distance of 1.45Å. The original ADAPT-VQE and the version with term removal correspond to the same color because the curves coincide. ................................... 147 xix Acronyms ADAPT Adaptive Derivative-Assembled Pseudo-Trotter vii,xv,xvi,xvii,xviii,xix,6,7,8,9,36,72, 73,74,75,76,82,83,94,96,97,98,99,112,113,115,116,117,118,119,122,123, 130,131,135,136,137,138,139,140,141,142,143,145,147,148,150,151,152 ADAPT Adaptive Derivative-Assembled Problem-Tailored xv,xvi,xvii,7,8,9,72,75,76,79,80,81, 82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,108, 109,110,111,112,113,114,150 BCH Baker–Campbell–Hausdorff 18,32 CC Coupled Cluster 5,31,32,33,61,72 CCSD Coupled Cluster Singles and Doubles 31 CI Configuration Interaction 30,31 CISD Configuration Interaction Singles and Doubles 31 CNOT Controlled NOT xiv,xvi,8,15,17,51,61,64,67,68,69,75,83,90,96,97,106,107, 112 COBYLA Constrained Optimization By Linear Approximation xiii,xiv,xvi,54,55,57,58,65,88 FCI Full Configuration Interaction xiv,xv,xviii,xix,30,48,61,62,63,69,82,84,85,116,139, 141,143 GSD Generalized Singles and Doubles xvii,101,112,113 NISQ Noisy Intermediate-Scale Quantum vii,3,4,5,7,8,19,20,34,75,86,96,150,152 QEB Qubit-Excitation-Based 102,143 q-UCCD Quantum Unitary Coupled Cluster Doubles 44 QAA Quantum Adiabatic Algorithm 34 QAOA Quantum Approximate Optimization Algorithm 34,38,39 QPE Quantum Phase Estimation 4,33,34 QUBO Quadratic Unconstrained Binary Optimization 34 xx ACRONYMS RDM Reduced Density Matrix 7 SAGSD Spin-Adapted Generalized Singles and Doubles 76 SCF Self-Consistent Field 27,28 SCGSD Spin-Complemented Generalized Singles and Doubles xvi,72,78,79,96,97 SGSD Singlet Generalized Singles and Doubles xv,xvi,xvii,xviii,xix,76,78,79,80,83,96,97, 116,117,122,131,136,137,139,140,141,144,148 SPAM State Preparation and Measurement xv,xvi,8,9,19,20,21,33,57,61,70,71,88,89, 90,91,150 SVD Singular Value Decomposition 49,50,52 UCC Unitary Coupled Cluster 5,32,44 UCCGSD Unitary Coupled Cluster Generalized Singles and Doubles 44 UCCSD Unitary Coupled Cluster Singles and Doubles vii,xiv,xv,xvi,5,6,7,8,9,32,33,36,44, 61,62,63,64,65,66,67,68,69,70,71,72,75,76,82,83,84,85,86,89,90,91, 92,93,94,150 UpCCGSD Unitary Paired Coupled Cluster Generalized Singles and Doubles 44 VCC Variational Coupled Cluster 32 VQA Variational Quantum Algorithm vii,3,4,8,34,38,54,61 VQE Variational Quantum Eigensolver vii,xiii,xiv,xv,xvi,xvii,xviii,xix,4,5,6,7,8,9,33,36, 38,39,43,44,46,48,49,52,53,54,55,56,57,58,59,61,62,63,64,65,66,67, 68,69,70,71,72,73,74,75,76,79,80,81,82,83,84,85,86,87,88,89,90,91,92, 93,94,95,96,97,98,99,100,101,102,108,109,110,111,112,113,114,115,116, 117,118,119,122,123,130,131,135,136,137,138,139,140,141,142,143,145, 147,148,150,151,152 xxi 1 Introduction 1.1 Context and Motivation Richard Feynman raised the question, What kind of computer are we going to use to simulate physics? , calling attention to the fact that the physical world is quantum mechanical, so quantum physics is what we are concerned with [25]. So arose the answer: a new kind of computer. The simulation of a quantum system being a difficult task to perform classically is but a symptom of the fact that bits, as basic units of information in classical computing, are unfit to describe nature under the laws of quantum mechanics. Quantum computing was thus born with the purpose of using quantum mechanics itself to perform computation, replacing bits by their quantum counterparts: qubits . Making use of quantum phenomena (namely entanglement; superposition; interference) empowers quantum computers with the ability to do operations without a classical analogue. The scaling of the description of highly entangled states using classical data inspired hope that quantum computers would be able to perform tasks that are intractable even to classical supercomputers. The memory required for storing a general state of a quantum system in a classical computer scales exponentially on the size of the system; in contrast, the scaling is linear on a quantum computer, as seems intuitive in view of it being such a system itself. In fact, an interesting reason to believe that a quantum computer may bring advantage is that there is no known classical algorithm that can simulate such a computer efficiently; representing the state of a few hundred qubits may already require more bits than the number of atoms in the visible universe [69]. However, this fact does not constitute in itself proof that quantum computers are more powerful than the classical ones. It is a misconception that quantum computers are exponentially parallel : even though superposition implies that qubits can be in a linear superposition of states that are exponentially many on the number of qubits, the information held by this superposition isn’t necessarily accessible. Measuring the state of the qubits will cause the collapse of the wave function: transforming quantum information into classical information destroys it irreversibly, unless the measurement reveals no information about the quantum state it is inflicted upon [64]. Famously, by Holevo’s theorem, the accessible classical information in n qubits is bounded above by n bits. The fact that the qubits carry more information is to no avail, unless one devises a clever way of extracting information in an efficient manner. As such, the design of techniques that allow employing 1 CHAPTER 1. INTRODUCTION quantum resources advantageously is all but straightforward; and the quest for quantum algorithms has developed into a large field of research. In 1980, Benioff published the first article on quantum computation, a pioneering work describing a quantum mechanical model of Turing machines [12]. Since then, there has been significant progress. Notable quantum algorithms were introduced in the following decade, such as Deutsch–Jozsa in 1992 [21] and Simon’s in 1994 [78]. Albeit not having much practical value, they hold their merit as some of the first examples of how a quantum computer could offer an exponential speedup over classical computers. Advantage of practical interest ensued when Shor, inspired by Simon’s proposal, introduced his famous factoring algorithm [76]. Shor’s algorithm enables a (hypothetical) quantum computer to factor integers in a time polynomial on their size - an exponential improvement as compared to the best-known classical algorithm for the task. While this remarkable speedup was not unprecedented, Shor’s algorithm owes its distinct success to the possible applications of the problem it was created to solve. An efficient factoring algorithm defies the security of the vastly used RSA cryptosystem, which brought attention to the possible implications of quantum computers. If resource requirements in what concerns quantum hardware were met, this algorithm would allow breaking internet security: it manages to solve in a matter of minutes a problem that would take years for the best-known classical algorithm to solve, and upon the difficulty of which the RSA algorithm relies [70]. This was a significant landmark, as well as a stronger motivation for developing the hardware that would pave the way for physically realising these theoretical proposals. Albeit offering an unarguably less astounding quadratic speedup, Grover’s algorithm [36] for searching in an unstructured database has also drawn formidable attention. Aside from covering a wide array of applications, its strength is in offering a provable speedup: evidently, the number of evaluations required to solve the problem classically has to scale linearly (half of elements in the search space will have to be checked on average). By offering the uncanny possibility of finding a marked element after a number of operations that is sub-linear on the size of the search space, Grover’s algorithm also contributed to a build up of the public reputation of quantum computing as a curious and strange field. While they seem compelling, these theoretical results are yet to be proven by experimental realization. Current hardware limitations hinder the implementation of large-scale instances of algorithms such as the aforementioned ones. There has been notable progress since the concept of quantum computing was introduced. The first implementation of an algorithm on a physical quantum computer took place in 1998 [18], and several proof-of-concept demonstrations of small-scale instances of different algorithms followed. In 2012, Shor’s algorithm was used to factor N=21 for the first time [58]. However, practical experiments up to now have been limited to toy problems and small, classically tractable demonstrations. Building a quantum computer capable of outperforming a classical one in tasks that may be of practical use seems to be strenuous at best, impossible at worse. The difficulty of building the physical hardware is indissociable from the very same nature on top of which the concept of quantum computing rests. Quantum mechanical states are very sensitive, causing 2 1.2. STATE OF THE ART qubits to require close to perfect isolation from the world to store and process information reliably. Building a quantum computer thus implies granting the proper conditions, i.e. conditions that reduce the coupling with the surrounding environment to the point that useful computations can be carried. Often, this involves temperatures near absolute zero and shielding from radiation; this is aggravated by the fact that the difficulty greatly increases with the system size, preventing large-scale quantum computers to have been built so far [70]. There is a dichotomy between the scaling of the most compelling quantum algorithms, and the scaling of the difficulty of building a quantum computer capable of actually bringing them to fruition. In the words of two prominent French physicists vocally skeptical of quantum computing [37], ...the large-scale quantum machine, though it may be the computer scientist’s dream, is the experimenter’s nightmare . After Shor’s algorithm was discovered, concerns were raised regarding the viability of quantum computing: algorithms built under the assumption of an ideal quantum computer might be rendered useless by decoherence and other types of noise [52,86]. In light of the no-cloning theorem [23,66,96], that forbids quantum data from being copied, error correction didn’t seem as trivial as in the classical case. In 1995, Peter Shor himself recovered some hope with a viable way of realising error correction in quantum computing [77]. Quantum error correcting codes are of the greatest importance: if error correction is important for computing in general, it is pivotal for quantum computing in particular. Unfortunately, the overhead in terms of qubit requirements and number of necessary gates is severe [31]. An error corrected, fault-tolerant quantum computer does not seem to belong to a near future. John Preskill coined the term Noisy Intermediate-Scale Quantum (NISQ) to describe the current era [69]. The name encompasses qubit limitations as well as noise and decoherence. After the promising algorithms of the 1990s shed some light on the potential of a fault-tolerant quantum computer, there has been a growing effort to build such a device - but that possibility seems at best distant. In an attempt to make as good use of NISQ computers as possible, the attention shifted in part from the asymptotic scaling of algorithms to their near-term viability. The focus has widened, and thus, in addition to creating conceptually advantageous algorithms lacking certainty that they are ever to be implemented, we now observe a growing dedication to devising algorithms that seem implementable but aren’t necessarily known to be advantageous. Variational quantum algorithms have emerged. 1.2 State of the Art VQAs are hybrid quantum-classical algorithms that resemble machine-learning techniques typical to classical computing. The usual outline involves using a classical computer to find the variational parameters that minimize a cost function, which is evaluated by the quantum computer. An interesting review of the structure, applications and challenges of variational quantum algorithms can be found in reference [17]. 3 2 Theoretical Framework 2.1 Quantum Computing The purpose of this section is to offer a simple overview of the concepts in quantum computing that are necessary for the following chapters. A more in-depth explanation on these topics can be found on the book Quantum Computation and Quantum Information, reference [64], upon which the following was based. 2.1.1 Basic Concepts 2.1.1.1 Pure Single Qubit States The basic unit of information for classical computing, the bit , undertakes two possible states, usually labeled 0and 1. Similarly, the quantum bit , or qubit , is a two-level system. The difference resides in the quantum mechanical nature of the qubit: while the classical bit is either in state 0or 1, the quantum bit can be in a superposition of these two states. A qubit can be physically implemented by a two-level quantum-mechanical system. Irrespective of the specific implementation, the state of an isolated qubit can be represented abstractly in Dirac notation as |𝜓i=𝛼|0i+𝛽|1i.(2.1) In the expression above |0iand |1iare the computational basis states ; they form an orthonormal basis that spans the state space of the system, a two-dimensional Hilbert space ℋ(a complex vector space equipped with inner product). By convention, |0i= 1 0!,|1i= 0 1!.(2.2) The numbers 𝛼and 𝛽are the respective complex amplitudes, related to the respective probability of occurrence. As discussed previously, even though they characterize the state of the qubit, it is not possible to discover the values of these numbers with access to a single copy of the state: when a measurement is done on the computational basis, the outcome will be a computational basis state. The squared absolute values of the amplitudes correspond to the probability of, upon measurement, obtaining state |0ior |1irespectively: 𝛼𝛼∗=|𝛼|2=Probability of |0i, 10 2.1. QUANTUM COMPUTING 𝛽𝛽∗=|𝛽|2=Probability of |1i. Since the probabilities must sum to one, we have a normalization condition that reads h𝜓|𝜓i=|𝛼|2+ |𝛽|2=1.(2.3) To describe the state of a perfectly isolated qubit, it thus suffices to specify the corresponding unit vector in ℋ. This can be done by specifying the values of 𝛼in 𝛽in equation 2.1 subject to condition 2.3. The normalization removes one degree of freedom and allows the state to be parameterized by three real numbers as |𝜓i=𝑒𝑖𝛾 cos 𝜃 2|0i+𝑒𝑖𝜑 sin 𝜃 2|1i. The parameter 𝛾is an irrelevant degree of freedom that does not impact any observable. Since the global phase associated with it is not detectable and lacks physical significance, it can be promptly ignored and removed from our parameterization altogether, leaving |𝜓i=cos 𝜃 2|0i+𝑒𝑖𝜑 sin 𝜃 2|1i.(2.4) With two degrees of freedom, this parameterized expression can now be visualized as lying on the surface of a sphere. The Bloch sphere depicted in figure 1is a popular and convenient geometric representation of the state of a (pure) qubit. Figure 1: Representation of the state of a qubit in the Bloch sphere. The angles 𝜃and 𝜑correspond to the parameters from equation 2.4. A qubit can be acted on by operators that are represented by two-by-two matrices. There are three matrices of particular importance, the Pauli matrices : 11 CHAPTER 2. THEORETICAL FRAMEWORK 𝑋≡𝜎𝑥≡ 0 1 1 0!, 𝑌 ≡𝜎𝑦≡ 0−𝑖 𝑖0!, 𝑍 ≡𝜎𝑧≡ 1 0 0−1!.(2.5) Frequently one considers the identity matrix 𝐼as a Pauli matrix itself, oftentimes denoted 𝜎0. All Pauli matrices are unitary (𝜎𝑖𝜎† 𝑖=𝜎† 𝑖𝜎𝑖=𝐼) and Hermitian (𝜎𝑖=𝜎† 𝑖). These four matrices are orthogonal to each other, and together they span the space of the observables of a qubit. Observables are named so because they are associated with observable physical quantities . Each projective measurement is associated with an observable, whose eigenvalues correspond to the possible outcomes of the measurement. These outcomes are always real-valued, which is related to the fact that all observables are Hermitian. From definitions 2.2 and 2.5 we can see that the computational basis states are eigenstates of 𝑍. Computational basis measurements are projective measurements associated with the observable 𝑍, and their possible outcomes correspond to the eigenvalues of the operator. These are the measurements typically available in quantum computers. 2.1.1.2 Pure Multiple Qubit States All these ideas generalize to larger systems: the tensor product can be used to obtain a larger vector space in which the state of the multiple qubits lives. The description grows rapidly: to specify the state of n qubits, we need 2𝑛complex amplitudes, that correspond to 2𝑛+1−2real-valued degrees of freedom after removing those corresponding to the normalization condition and the global phase. When we have more than one qubit, measurementoutcomes can be correlated beyond what is possible in classical systems. This is because it is not always the case that an n qubit state can be written in the form |𝜓0i⊗|𝜓1i⊗... ⊗|𝜓𝑖i⊗... ⊗|𝜓𝑛−2i⊗|𝜓𝑛−1i,(2.6) where the 𝜓𝑖are the states of the individual qubits. When this is in fact possible, the state is called separable ; otherwise, we say that it is entangled . Entanglement plays a big role in quantum algorithms: it is an exclusively quantum-mechanical resource and a bedrock for quantum speedup. For two qubits, the maximally entangled states are the famous Bell states : 𝜓+=1 √2|01i+|10i,|𝜓−i=1 √2|01i−|10i, 𝜙+=1 √2|00i+|11i,|𝜙−i=1 √2|00i−|11i. (2.7) It can be readily seen that the outcome of a (computational basis) measurement done on one of the qubits fully determines the outcome of a measurement on the other one. 12 2.1. QUANTUM COMPUTING 2.1.1.3 Mixed States Most of what was discussed so far only holds for pure states - those that are known exactly, and can be written in the form of equation 2.4. During a computation, one would ideally deal with pure states exclusively. Unfortunately, this is never the case: there are unknown errors occurring during the state preparation and the computation itself. The qubits are not fully isolated, as there is always some coupling with the environment. This results in the system being in a mixed state , a statistical distribution over pure states. This mixture simply translates lack of knowledge about the system. States that are not pure cannot be written as in equation 2.4 and cannot be represented in the surface of the Bloch sphere. The density matrix 𝜌provides a way to describe a system in a more general state as 𝜌=Õ 𝑖 𝑝𝑖|𝜓𝑖ih𝜓𝑖|.(2.8) Here the 𝑝𝑖are the probabilities of occurrence of the respective |𝜓𝑖i, the possible pure state configurations. The density matrix offers a way to test the purity of the state: tr(𝜌2)=1−→ state is pure; tr(𝜌2)<1−→ state is mixed. (2.9) There is also a way of generalizing the geometrical representation of the states: the Bloch ball . We rewrite the density matrix representing the state as 𝜌=𝐼+®𝑟· ®𝜎 2,(2.10) where ®𝜎=(𝑋,𝑌, 𝑍). The three-dimensional vector ®𝑟, with norm lower than or equal to unity, is the Bloch vector whose coordinates are written in 2.11. Õ 𝑖 𝑝𝑖𝑥𝑖,Õ 𝑖 𝑝𝑖𝑦𝑖,Õ 𝑖 𝑝𝑖𝑧𝑖!(2.11) The sum is over the possibilities of pure states, with 𝑝𝑖their probabilities and (𝑥𝑖,𝑦𝑖, 𝑧𝑖)their coordinates on the Bloch sphere. It is easy to verify that, for the special case of pure states, the norm of the Bloch vector ®𝑟is 1 and the coordinates are as before. For mixed states, the norm is in [0,1[: they are inside the Bloch sphere, in the Bloch ball . The lower the purity, the lower the norm of ®𝑟. The totally mixed state lies in the center of the Bloch sphere, representing a complete lack of knowledge about the system. 2.1.2 Measuring Pauli Strings As it was mentioned before, one can perform a 𝑍measurement on a qubit. This yields either of the operator’s eigenvalues, with the respective probabilities depending on the state of the qubit. Repeated measurements allow obtaining the expectation value of this operator in this specific state, namely 13 CHAPTER 2. THEORETICAL FRAMEWORK h𝜓|𝑍|𝜓i. It is common for problems to require measurements in bases other than the computational basis; this can be done with basis rotations . Measuring a generic observable amounts to performing a computational basis measurement, preceded by the unitary that rotates from the eigenbasis of the desired observable to the computational basis. Evidently, in order to be appropriate for an observable 𝐴, the unitary 𝑈needs to obey the condition h𝜓|𝑈†𝑍𝑈 |𝜓i=h𝜓|𝐴|𝜓i for all states |𝜓i, so that the requirement can be written simply as 𝑈†𝑍𝑈 =𝐴. (2.12) This generalizes for observables acting on multiple qubits. Implementing the basis rotation unitary for a generic observable may prove to hinder efficiency. Fortunately, Pauli measurements in particular are very simple to do. (a) Circuit to measure 𝑍. (b) Circuit to measure 𝑋. (c) Circuit to measure 𝑌. Figure 2: Examples of circuits that perform single-qubit Pauli measurements, drawn in the IBM Quantum Composer [41]. Circuits for all three types of single-qubit Pauli measurements (𝑍,𝑋,𝑌) are drawn in figure 2. The circuit in 2a implements the typical computational basis measurement: no basis rotation is necessary. For 𝑋we can use the unitary 𝑈=𝐻as the basis rotation - this is the popular Hadamard gate. The circuit is represented in figure 2b. Finally, figure 2c implements a measurement of 𝑌, with the basis rotation being a rotation around the x axis 𝑈=Rx(𝜋/2). 14 2.1. QUANTUM COMPUTING These choices of unitaries are not unique. As an example, using 𝑈=𝐻𝑆†for measuring 𝑌is also a popular option. One can choose the rotation that is the most convenient, given the operations that can natively be implemented in the hardware. As well as single qubit Pauli measurements, one can easily measure Pauli strings . Pauli strings are observables of the form Ì𝜎𝑖𝑖∈ {0, 𝑧, 𝑥,𝑦}.(2.13) There are multiple ways to perform such measurements. The simplest and most evident one amounts to performing the individual Pauli measurements and multiplying the outcomes. However, it is possible to include two qubit gates in addition to the single qubit rotations as a means of reducing the number of measurements that needs to be performed. As an example, two versions of circuits to perform a 𝑍⊗𝑍measurement are illustrated in figure 3. (a) Using two single-qubit measurements. (b) Using one single-qubit measurement. Figure 3: Examples of circuits that perform a 𝑍⊗𝑍measurement, drawn in the IBM Quantum Composer [41]. The circuit in 3a measures both qubits in the computational basis; the result is then obtained by multiplying the outcomes. The final value will be in {−1,1}. The actual value depends on the parity of the computational basis state that was measured: we obtain 1 if the parity was even, meaning that the qubits were measured on the same state (|00ior |11i). The outcome is -1 if the parity was odd, meaning that they were measured on opposite states (|01ior |10i). In contrast, the circuit in 3b performs only one measurement. A CNOT gate is used to compute and store the parity of the state of both qubits in the state of the first one alone. Essentially, for a two-qubit observable 𝐴, we can use an unitary that obeys 𝑈†(𝑍⊗𝑍)𝑈=𝐴 and measure both qubis, or use an unitary that obeys 𝑈†(𝑍⊗𝐼)𝑈=𝐴 and measure only one of them. 15 CHAPTER 2. THEORETICAL FRAMEWORK With this, it has become clear that any Pauli string, i.e. operator of the form of equation 2.13, can be measured in a quantum computer. This is of great importance, because these Pauli strings span the space of Hermitian operators, and all observables are Hermitian. As a consequence, all observables can be written as a linear combination of Pauli strings; and given the linearity of quantum mechanics, all expectation values of observables can be obtained as a weighed sum of expectation values of Pauli strings. If we want to measure an observable ˆ 𝑂, we can start by decomposing it in the form of equation 2.14. ˆ 𝑂=Õ 𝑖 ℎ𝑖ˆ 𝑃𝑖(2.14) Here the ℎ𝑖are real coefficients and the ˆ 𝑃𝑖are Pauli strings (i.e. observables of the form of 2.13). It is then straightforward to evaluate the desired expectation value: it amounts to measuring the expectation value of each Pauli string appearing in the decomposition of the observable, and plugging it into formula 2.15. hˆ 𝑂i=Õ 𝑖 ℎ𝑖hˆ 𝑃𝑖i(2.15) This procedure can be followed for any observable; however, it is only tractable if the observable can be decomposed into a sum of polynomially-many Pauli strings. 2.1.3 Quantum Simulation and Trotterization The behaviour of an isolated quantum-mechanical system is dictated by the Schrödinger equation in 2.16. 𝑖ℏ𝑑 𝑑𝑡 |𝜓i=𝐻|𝜓i(2.16) When the Hamiltonian operator for the system, 𝐻, is time-independent, this implies that the time evolution of the wave function has the form presented in 2.17. |𝜓(𝑡)i=𝑒−𝑖𝐻𝑡 |𝜓(0)i(2.17) As it was mentioned before, the simulation of physical systems is one of the most promising applications of quantum computers, as it is generally a hard task for classical computers. As such, it is important to know how to simulate this time evolution on a quantum computer. A general Hamiltonian is not easy to exponentiate - so the question becomes, how do we approximate the time evolution with sufficient accuracy? Typically, one deals with special classes of Hamiltonians that allow easy and efficient exponentiation. A Hamiltonian can be decomposed into a sum of several, simpler terms: 𝐻= 𝐿 Õ 𝑘=1 𝐻𝑘.(2.18) 16 2.1. QUANTUM COMPUTING And the 𝐻𝑘may allow easier implementation on a quantum computer - for example, each one may act only on a small subsystem, or they may have a special form. In fact, as was discussed in subsection 2.1.2, we know that any Hermitian operator (such as the Hamiltonian) can be decomposed into a linear combination of Pauli strings. Conveniently, when the 𝐻𝑘 are Pauli strings with real coefficients, it is easy to construct the circuit that applies 𝑒−𝑖𝐻𝑘𝑡. (a) Circuit for simulating the time evolution of a system under the Hamiltonian 𝐻𝑘=𝑍⊗ 𝑍⊗𝑍⊗𝑍. (b) Circuit for simulating the time evolution of a system under the Hamiltonian 𝐻𝑘=𝑋⊗ 𝑌⊗𝑍⊗𝑋. Figure 4: Examples of circuits that apply time evolution under simple Hamiltonians 𝐻𝑘, consisting of Pauli strings. They apply 𝑒−𝑖𝐻𝑘𝑡using a parameterized 𝑍rotation. The circuits were drawn in the IBM Quantum Composer [41]. Figure 4shows two circuits that can be used to implement the exponentiation of two specific 𝐻𝑘. The circuit in figure 4a is the simplest one: it implements 𝑒−𝑖𝑡𝑍 ⊗𝑍⊗𝑍⊗𝑍. This amounts to applying a local phase 𝑒−𝑖𝑡 to the computational basis states that have an even parity, and 𝑒𝑖𝑡 to those that have an odd parity; the drawn circuit does precisely that. The first ladder of CNOT gates computes the parity of the state and stores it on the last qubit. Then, all that’s left is to apply 𝑒−𝑖𝑡𝑍 to the qubit that holds the parity, which can be done by rotating it around the 𝑍axis by an angle of 2𝑡(i.e. applying the gate Rz(2𝑡)). Finally, a ladder of CNOT gates is applied to uncompute the parity before the computation proceeds. 17 CHAPTER 2. THEORETICAL FRAMEWORK In figure 4b, we have the circuit that implements 𝑒−𝑖𝑡𝑋 ⊗𝑌⊗𝑍⊗𝑋. This can be done with minor adjustments to the previous circuit. Appendix Ashows that, when 𝐻𝑘is a Pauli string, one can write any operator of the form 𝑒𝑖𝑡𝐻𝑘as 𝑒−𝑖𝑡 Ë𝑖𝑃𝑖=Ì 𝑖 𝑈† 𝑖𝑒−𝑖𝑡 Ë𝑖𝑍𝑖Ì 𝑖 𝑈𝑖,(2.19) where the 𝑈𝑖are single qubit unitaries obeying 𝑈† 𝑖𝑍𝑖𝑈𝑖=𝑃𝑖. For any qubit 𝑖, acted on by 𝑃𝑖in 𝐻𝑘. This is the same condition as 2.12: as we did to make Pauli measurements, we need an unitary that rotates from the eigenbasis of the operator in question to the computational basis. The form of equation 2.19 is very convenient: it tells us that we can use the circuit for 𝑒−𝑖𝑡 Ë𝑖𝑍𝑖(figure 4a) for implementing any operator of the form 𝑒−𝑖𝑡 Ë𝑖𝑃𝑖(where 𝑃𝑖∈ {𝐼𝑖, 𝑍𝑖, 𝑋𝑖, 𝑌𝑖}), as long as we apply the appropriate Ë𝑖𝑈𝑖before and Ë𝑖𝑈† 𝑖after. The circuit in figure 4b does precisely that, for the specific case of the Pauli string Ë𝑖𝑃𝑖=𝑋⊗𝑌⊗ 𝑍⊗𝑋. The employed basis rotations were the same as those from section 2.1.2. So now we know how to simulate the time evolution under simple Hamiltonians 𝐻𝑘consisting of a single Pauli string. However, the end goal is simulating the full Hamiltonian 𝐻from 2.18. This amounts to applying 𝑒−𝑖𝑡 Í𝐿 𝑘=1𝐻𝑘, but 𝑒−𝑖𝑡 Í𝐿 𝑘=1𝐻𝑘≠ 𝐿 Ö 𝑘=1 𝑒−𝑖𝑡𝐻𝑘 in general. In fact, for two operators 𝐻𝑖,𝐻𝑗, we can calculate 𝑒𝐻𝑖𝑒𝐻𝑗as 𝑒𝐻𝑖𝑒𝐻𝑗=𝑒𝐻𝑖+𝐻𝑗+1 2[𝐻𝑖,𝐻𝑗]+ 1 12 [𝐻𝑖,[𝐻𝑖,𝐻𝑗]]− 1 12 [𝐻𝑗,[𝐻𝑖,𝐻𝑗]]+....(2.20) Equation 2.20 is known as the Baker–Campbell–Hausdorff (BCH) formula. We can see that the product of the exponential series is the exponential series of the sums only in the special case of commuting operators ([𝐻𝑖, 𝐻𝑗]=0). Fortunately, we can use the Trotter formula : lim 𝑛→∞(𝑒𝑖𝐻𝑖𝑡/𝑛𝑒𝑖𝐻𝑗𝑡/𝑛)𝑛=𝑒𝑖(𝐻𝑖+𝐻𝑗)𝑡.(2.21) Evidently, it is not feasible to implement the circuit for 𝑛→ ∞. However, we can make 𝑛finite to obtain a circuit that approximates the desired time evolution. 𝑒𝑖(𝐻𝑖+𝐻𝑗)𝑡≈ (𝑒𝑖𝐻𝑖𝑡/𝑛𝑒𝑖𝐻𝑗𝑡/𝑛)𝑛(2.22) 18 2.1. QUANTUM COMPUTING Here, 𝑛is the number of repetitions of the basic unit of our Trotter circuit, and 𝑡/𝑛is the time step . We can interpret equation 2.22 as approximating the evolution of a system under (𝐻𝑖+𝐻𝑗)for a time 𝑡with 𝑛repetitions of the evolution under 𝐻𝑗followed by the evolution under 𝐻𝑖, both always for a time step of 𝑡/𝑛. We switch between evolving the system under the two elements in the sum repeatedly, with the duration of this evolution being always less than the actual evolution time. The larger the number of repetitions 𝑛, the shorter the time step and the better the approximation: in the limit of 𝑛→ ∞, one would get the exact result. The error in the approximation will depend on the evolution time; as an example, a Trotter approximation with a single repetition will have an error O(𝑡2) [64]. 2.1.4 Current Limitations As discussed in the previous chapter, it will likely be a while before a fault-tolerant quantum computer (FTCQ) is available. Currently, one has to deal with noisy intermediate-scale quantum (NISQ) devices. These are not error corrected and come with a multitude of limitations, namely: • Coherent noise. • Incoherent noise. •SPAM errors. A brief analysis of each follows. 2.1.4.1 Coherent Noise Coherent errors are unitary and don’t undergo fast changes (relative to the gate time); this type of error has various causes, such as systematic control noise, global external fields, cross-talk, and unwanted qubit-qubit interactions [33]. For a single qubit, a coherent source of noise could be performing a rotation by an angle of 𝜃+ 𝛿𝜃 instead of just 𝜃due to miscalibration. For two qubits, a well-known example is the unwanted 𝑍𝑍 interaction in superconducting quantum computers. It can significantly lower the fidelity of two-qubit gates; as such, it can have a great impact on the performance of the device, and consequently pose a challenge on the scalability of the architecture [100]. Coherent errors interfere constructively, and may take a significant toll on the computation. In the worst case scenario, the error rate of coherent noise can scale quadratically both in the number of qubits and in the length of the circuit [74,44]. For that reason, it can accelerate the error accumulation rate and be a challenge in the context of error correction codes. Regardless, there is hope that error correction can be used to remove undesirable coherence from the noise channel, thus avoiding the possibility that constructive interference leads to a quadratic scaling of the error rate. For the toric code and a particular 19 CHAPTER 2. THEORETICAL FRAMEWORK Typically, the methods for finding an approximate solution to the ground state of the electronic system aim to minimize h𝜓|ˆ 𝐻|𝜓i. They differ only in the variational form of |𝜓i, a concept identical to that of an ansatz . It is a guess that fixes the structure of the wave function, but only specifies it up to some tunable variational parameters . Plugging some test parameters into the variational form, we obtain a fully specified wave function in which we can calculate the expectation value of ˆ 𝐻. The task is then reduced to finding the best parameters: given the variational principle, this is just a minimization problem. A more thorough overview of the variational principle can be found in reference [34]. 2.2.2.1 Self-Consistent Mean Field Methods As much as the Born-Oppenheimer approximation greatly simplifies our eigenproblem, solving the Schrödinger equation with the electronic Hamiltonian (equation 2.30) exactly is still out of reach. The problematic term is now the electron-electron repulsion (expression 2.27), which couples the coordinates of all the electrons. Finding the exact solution implies dealing simultaneously with the coordinates of them all; decoupling these variables would allow solving the equation. Such is made possible by the mean-field approximation . The essence of the mean-field approximation resides in the simplification of the effect of the electrons in each other: instead of contemplating the sum of all the individual repulsion terms, we assume that each electron experiences an average potential created by all the other electrons. Each one of them is considered to be moving in the mean field created by the others. This gives rise to the operator in equation 2.32. ˆ 𝑓(𝑖)=−1 2∇2 𝑖− 𝑀 Õ 𝐴=1 𝑍𝑎 𝑟𝑖𝐴 +𝑣(𝑖)(2.32) This operator is just the electronic Hamiltonian in equation 2.30, separated for each electron 𝑖. The impact of the remaining electrons is contained within 𝑣(𝑖), an average potential. At this point we have managed to separate the Schrödinger equation into single fermion equations: instead of one complicated many-body problem, we face multiple simple one-body problems: ˆ 𝑓(𝑖)𝜒(®𝑥𝑖)=𝜖 𝜒(®𝑥𝑖).(2.33) This is another eigenvalue equation, thus its form is identical to that of equation 2.23. But it is much simpler: we now have ˆ 𝑓(𝑖), a one-body operator, instead of the many-body operator ˆ 𝐻. Instead of the many-body wave function |𝜓i, depending on the coordinates of all of the electrons, we have single-body wave functions (spin-orbitals) 𝜒(®𝑥𝑖), each of which depends solely on the coordinates of one of them. If we have N electrons, we have N equations of the form of 2.33. We want to find the eigenvalues 𝜖and the eigenfunctions 𝜒(®𝑥𝑖)of each operator ˆ 𝑓(𝑖). However, this set of equations is non-linear: for the 𝑖th electron this operator depends (through 𝑣(𝑖)) on the wave functions of the other electrons 𝜒(®𝑥𝑗), which are the solutions of all the other N-1 equations. Thus, the solution must be obtained iteratively . 26 2.2. QUANTUM CHEMISTRY The iterative method for solving the equation 2.33 is the Self-Consistent Field (SCF) method. This is because the convergence condition is the self-consistency of the field. The variational parameters are contained within the single-particle wave functions, the 𝜒(®𝑥𝑖). We start with an initial guess for the values of these parameters. With a concrete set of single-particle wave functions, we are then in place to calculate the average potential that each particle 𝑖feels (𝑣(𝑖)). Plugging the values in the expressions of the operator ˆ 𝑓(𝑖), we get N single-particle equations of the form of 2.33. These are much simpler to solve than the full Schrödinger equation. The obtained solutions are new single-particle wave functions 𝜒(®𝑥𝑖). From these we can calculate new average potentials 𝑣(𝑖), and solve for the eigenfunctions of the operator ˆ 𝑓(𝑖)once again. The procedure is repeated until self-consistency is achieved, i.e. the fields 𝑣(𝑖)stop changing and the trial functions 𝜒(®𝑥𝑖)used to calculate those fields are also the solutions to the equations written with them. At this point, the trial functions used to construct the operator ˆ 𝑓(𝑖)are simultaneously its eigenvalues: the field is self-consistent. 2.2.2.2 Variational Form: Hartree and Hartree-Fock Methods The previous discussion purposefully glossed over what the many-body wave functions look like, and how they are related to the one-body wave functions. While the essence of SCF methods is similar, they differ significantly in this one key aspect: the ansatz , i.e. the variational form we choose for the wave functions. It seems like the simplest approach would be to consider that the many-body wave functions are just products of one-body wave functions: 𝜓(®𝑥1,®𝑥2, ..., ®𝑥𝑁−1,®𝑥𝑁)=𝜒𝛼(®𝑥1)𝜒𝛽(®𝑥2)...𝜒𝜋(®𝑥𝑁−1)𝜒𝜌(®𝑥𝑁).(2.34) This is what Hartree proposed upon introducing an SCF method. However, the Hartree method does not provide satisfactory results, because it does not account for a fundamental principle: the antisymmetry principle . The antisymmetry principle is a crucial axiom in quantum mechanics. It states that the many-body fermionic wave function must be antisymmetric under exchange of any two particles, which we can write as 𝜓(®𝑥1, ..., ®𝑥𝑖, ..., ®𝑥𝑗, ..., ®𝑥𝑁)=−𝜓(®𝑥1, ..., ®𝑥𝑗, ..., ®𝑥𝑖, ..., ®𝑥𝑁).(2.35) Evidently, the Pauli exclusion principle is enforced by the antisymmetry condition 2.35. The condition is rooted in the indistinguishability of particles. Exchanging indistinguishable, or identical , particles must result in a physically equivalent state. Such only allows variability up to a global phase. There are two main types of indistinguishable particles: fermions and bosons, with complex phase factors 𝑒𝑖𝜋 =−1and 𝑒𝑖2𝜋=1respectively. These factors give rise to the antisymmetric nature of fermionic wave function, and the symmetric nature of the bosonic wave function. 27 CHAPTER 2. THEORETICAL FRAMEWORK In either case, the simplest possible way of combining the individual wave functions, the Hartree product (definition 2.34), is inadequate. In the scenario we are concerned with (the fermionic one), this can be dealt with by using antisymmetrized products instead. 𝜓(®𝑥1,®𝑥2, ..., ®𝑥𝑁−1,®𝑥𝑁)=  𝜒1(®𝑥1)𝜒2(®𝑥1)... 𝜒𝑁−1(®𝑥1)𝜒𝑁(®𝑥1) 𝜒1(®𝑥2)𝜒2(®𝑥2)... 𝜒𝑁−1(®𝑥2)𝜒𝑁(®𝑥2) . . .. . ..... . .. . . 𝜒1(®𝑥𝑁−1)𝜒2(®𝑥𝑁−1)... 𝜒𝑁−1(®𝑥𝑁−1)𝜒𝑁(®𝑥𝑁−1) 𝜒1(®𝑥𝑁)𝜒2(®𝑥𝑁)... 𝜒𝑁−1(®𝑥𝑁)𝜒𝑁(®𝑥𝑁)  (2.36) Wave functions of the form of 2.36 are typically called Slater determinants . Matrix determinants are antisymmetric under exchange of any two rows or columns, making them a natural choice here. For bosons, one would use matrix permanents . Slater determinant ansätze were used in the Hartree-Fock method, an updated version of the Hartree method that accounted for the antisymmetry principle. As it was mentioned before, the choice of ansatz is what shapes a specific SCF method. The form of the operator ˆ 𝑓(𝑖)as written in 2.32, and its place in the algorithm, are common to several SCF methods. The difference is in the potential 𝑣(𝑖): it will be specific to the assumptions we make about the wave function, i.e. it depends on how our ansatz looks like. This operator is then used to write Schrödinger-like one-particle equations 2.33, which means that the chosen ansatz greatly impacts the formulation of our eigenproblems. In the Hartree-Fock method, the operator of the form of equation 2.32 is called the Fock operator : ˆ 𝑓(𝑖)=−1 2∇2 𝑖− 𝑀 Õ 𝐴=1 𝑍𝑎 𝑟𝑖𝐴 +𝑣𝐻𝐹 (𝑖)(2.37) The key difference is that the 𝑣𝐻𝐹 (𝑖)term here has a superscript, referring to the fact that the potential in this operator is the Hartree-Fock potential. The eigenvalue equations 2.33 defined through the Fock operator are the Hartree-Fock equations. As an illustration of how the potential changes depending on the ansatz, a shallow, simplified comparison of the Hartree and Hartree-Fock potentials follows. The method proposed by Hartree used the Hartree product as the variational form. This was associated with a simple potential, consisting of a local operator that accounted solely for the electron-electron repulsion. In contrast, the potential in Hartree-Fock method includes, in addition to this local operator, an exchange operator. This is a non-local operator that simply results from choosing Slater determinants as wave functions: it is an artifact associated with the imposition of the anti-symmetry principle. 28 2.2. QUANTUM CHEMISTRY 2.2.2.3 Limitations of the Hartree-Fock Method The Hartree-Fock method simplifies the problem in multiple ways. Demanding that the solution is a Slater determinant is in itself a restriction that precludes exactness, as is typical of variational methods: the optimal solution is only guaranteed to be the best within the variational form . What the Hartree-Fock method does is to find the best solution of Slater determinantal form. As was discussed before, the motion of the nuclei is neglected (Born-Oppenheimer approximation). Further, the method also assumes the mean-field approximation, which fails to properly account for electron correlation outside of antisymmetry. Finally, even if none of these approximations were used, the exact solution would still be beyond reach, because for the method to be tractable the basis set must be finite. For the Slater determinants to constitute a complete N-electron basis, they would have to be built from an infinite basis set of spin-orbitals. Regardless, the Hartree-Fock method is a cornerstone in quantum chemistry, and its approximate solution is often used as a starting point in more accurate methods, both in classical and quantum computing. 2.2.2.4 Post-Hartree-Fock Methods Given the shortcomings of the Hartree-Fock solution, a multitude of more accurate methods has been proposed to improve upon the Hartree-Fock approximation (at the expense of extra computational cost). The following abbreviated explanation of post-Hartree-Fock methods is based, not only on the main reference for this section ([84]), but also on references [71] and [38]. Several relevant simplifications assumed by the Hartree-Fock method were mentioned previously. Of these, one typically undone by posterior methods is the assumption that the solution is of the form of a Slater determinant. To correct for the electronic correlation ignored by the Hartree-Fock method, a typical approach is to expand the wave function as a linear combination of Slater determinants, rather than a single one of them. Since the Hartree-Fock solution is a good starting point, the wave function is often described in terms of excitation operators that will act on this state. It is useful to define the cluster operator 𝑇: 𝑇= 𝑁 Õ 𝑖=1 𝑇𝑖, 𝑇1=Õ 𝑖,𝑎 𝑡𝑖 𝑎𝑎† 𝑎𝑎𝑖 𝑇2=Õ 𝑖>𝑗,𝑎>𝑏 𝑡𝑖𝑗 𝑎𝑏𝑎† 𝑎𝑎† 𝑏𝑎𝑖𝑎𝑗 . . . (2.38) 29 CHAPTER 2. THEORETICAL FRAMEWORK The lower case𝑡are expansion coefficients. The indices𝑖,𝑗and𝑎,𝑏run over orbitals that are occupied and unoccupied (virtual) in the reference state, respectively. The operator 𝑎† 𝑘is a creation operator, which creates a fermion in orbital 𝑘;𝑎𝑘is an annihilation operator, which removes a fermion from orbital 𝑘. These operators belong to the second quantization formalism. The cluster operator consists of a sum of all possible excitation operators 𝑇𝑖: the one that generates single excitations from the reference state (𝑇1), the one that generates double excitations (𝑇2), and so forth. With this operator, it is possible to write the FCI variational form in 2.39. |𝐹𝐶𝐼i=(1+𝑇)|𝐻𝐹i(2.39) This is just the Hartree-Fock state acted on by an operator (1+𝑇); the expansion coefficients will be the variational parameters. The term 𝑇|𝐻𝐹iconsists of a linear combination of all possible excited determinants, since the cluster operator includes all possible excitation operators. The FCI state is then a linear combination of all Slater determinants in the N-electron space formed from the chosen set of spinorbitals. The optimal FCI energy 2.40 will be the exact solution within the N-electron subspace spanned by these determinants. Of course, other than the finite basis set, the solution is also only exact up to the Born-Oppenheimer approximation. 𝐸𝐹𝐶𝐼 =min ® 𝑡 h𝐹𝐶𝐼 |𝐻|𝐹𝐶𝐼i h𝐹𝐶𝐼 |𝐹𝐶𝐼i(2.40) In short, the inclusion of more Slater determinants increases the variational freedom and improves the solution. The name configuration interaction arises from the fact that each Slater determinant in the expansion is associated with a specific electronic configuration (configuration of spin-orbitals). Full refers to the fact that all excitation operators are included. Unfortunately, while FCI offers accurate results, it is only tractable for small molecules: the number of determinants that need to be included in the expansion grows with the factorial of the total number of spin orbitals. Consequently, even with minimal single electron basis sets, the difficulty of calculations increases rapidly as the size of the system increases. To make the calculations tractable, the cluster operator 𝑇must be truncated. 𝑇(𝑘)= 𝑘 Õ 𝑖=1 𝑇𝑖(2.41) The operator 𝑇(𝑘)defined in 2.41 is just the cluster operator from definition 2.38, truncated up to excitations of order 𝑘. When a truncated version of the operator is used, only a portion of the Slater determinants formed from the one-electron basis set will appear in the expansion; the associated method is simply called Configuration Interaction (CI). The variational form is written as in 2.42. |𝐶𝐼i=(1+𝑇(𝑘))|𝐻𝐹i(2.42) 30 2.2. QUANTUM CHEMISTRY For the specific case that 𝑘is two, the method is called Configuration Interaction Singles and Doubles (CISD), because it includes only single and double excitations from the reference state (variational form in 2.43). |𝐶𝐼𝑆𝐷i=(1+𝑇1+𝑇2)|𝐻𝐹i(2.43) Truncating the operator to make the computation tractable results in a few issues. Namely, if we have two independent, non-interacting subsystems (e.g. two infinitely separated molecules), the CI energy of the full system will not be the sum of the energies of the subsystems. Because it lacks this property, it is said that CI is not size-consistent . A size-consistent version can be created by means of exponentiation. The corresponding method is called CC; the variational form is written in 2.44. |𝐶𝐶i=𝑒𝑇|𝐻𝐹i(2.44) Once again, the cluster operator 𝑇is typically truncated at some order of excitation. For the usual choice that excitations up to order two are included, the method is called Coupled Cluster Singles and Doubles (CCSD) (variational form in 2.45). |𝐶𝐶𝑆𝐷i=𝑒𝑇1+𝑇2|𝐻𝐹i(2.45) Even when truncated, CC methods don’t suffer from the problem of size-inconsistency. Further, they have the interesting property that even when the cluster operator𝑇is truncated up to order 𝑘, higher-order excitations occur in the exponential series. However, they have some problems of their own. A relevant weakness in CC is that the similarity-transformed Hamiltonian ¯ 𝐻for conventional CC (definition 2.46) is not variational. ¯ 𝐻=𝑒−𝑇(𝑘)𝐻𝑒𝑇(𝑘)(2.46) This is due to the CC expectation value 𝐸𝐶𝐶 =h𝐻𝐹 |𝑒−𝑇(𝑘)𝐻𝑒𝑇(𝑘)|𝐻𝐹i(2.47) not being symmetric, since h𝐻𝐹 |𝑒−𝑇(𝑘)≠(𝑒𝑇(𝑘)|𝐻𝐹i)†, which is just a consequence of the operator 𝑒𝑇(𝑘)not being unitary. This means that we can’t use the variational principle 2.31, because the associated expectation value isn’t of the proper form. The asymmetry of the expectation value impedes the CC energy from being an upper bound to the ground energy. 31 CHAPTER 2. THEORETICAL FRAMEWORK To circumvent this, Variational Coupled Cluster (VCC) has been proposed. This method aims to minimize 𝐸𝑉𝐶𝐶 =h𝐻𝐹 |𝑒𝑇(𝑘)†𝐻𝑒𝑇(𝑘)|𝐻𝐹i h𝐻𝐹 |𝑒𝑇(𝑘)†𝑒𝑇(𝑘)|𝐻𝐹i .(2.48) Here, the denominator is not unity because of the non-unitarity of 𝑒𝑇(𝑘). The VCC method is obviously variational, because the numerator is of the form h𝜓|𝐻|𝜓i; further, much like traditional CC, it is sizeconsistent. Unfortunately, the cost of the VCC method scales exponentially with the system size, regardless of the truncation order 𝑘. Another approach to make a variational version of VCC is UCC, that replaces the operator 𝑒𝑇with another exponential operator, now unitary: |𝑈𝐶𝐶i=𝑒𝑇−𝑇†|𝐻𝐹i.(2.49) The Hermitian conjugate of the operator (𝑇−𝑇†)reads (𝑇−𝑇†)†=𝑇†−𝑇††=𝑇†−𝑇=−(𝑇−𝑇†), which means that this operator is anti-Hermitian, and the commutators between (𝑇−𝑇†)and (𝑇− 𝑇†)†are all zero. This grants unitarity to the UCC operator. 𝑒𝑇−𝑇†𝑒𝑇−𝑇††=𝑒𝑇−𝑇††𝑒𝑇−𝑇†=𝐼 Since the UCC operator is unitary, the method is variational, and the UCC energy is an upper bound to the ground state energy. Finding the UCC ground state then amounts to a minimization problem (expression 2.50). 𝐸𝑈𝐶𝐶 =min ® 𝑡h𝐻𝐹 |𝑒−(𝑇−𝑇†)𝐻𝑒𝑇−𝑇†|𝐻𝐹i(2.50) As usual, the cluster operator should be truncated for the computation to be tractable, so that we get an UCC operator of the form 𝑒𝑇(𝑘)−𝑇(𝑘)†. For the frequent case that only excitations up to order two are included (𝑘=2), the method is called UCCSD. This corresponds to the variational form in 2.51. |𝑈𝐶𝐶𝑆𝐷i=𝑒(𝑇1+𝑇2)−(𝑇† 1+𝑇† 2)|𝐻𝐹i(2.51) The UCCSD method is variational, size-consistent, and has more Slater determinants contributing to the wave function than CISD, even though the truncation order of 𝑇is the same. In spite of its benefits, the UCCSD ansatz is not used in classical computations because it can’t be efficiently evaluated. Regardless of the truncation order of the cluster operator, expression 2.50 is classically intractable, because the BCH expansion for 𝑒−(𝑇−𝑇†)𝐻𝑒𝑇−𝑇†is infinite. 32 2.3. VARIATIONAL QUANTUM ALGORITHMS What is interesting is that while originally the unitary variant of the CC operator was created with the purpose of obtaining a variational version of CC theory, this unitarity has the convenient side effect of allowing straightforward implementation of the variational form in quantum computers. In fact, the UCCSD ansatz was used in the article that introduced the VQE (reference [67]). Further, this ansatz has inspired many other proposals of variational forms for use in chemistry applications. 2.3 Variational Quantum Algorithms This section will present the general outline of variational quantum algorithms, as well as the context that motivated their ingress into the field of quantum computing and prompted their popularity. An overview of the subject can be found in the article [17], used as the main reference throughout the subsections to follow. 2.3.1 Motivation The main incentives for devising variational quantum algorithms are no other than the main limitations of quantum computers. The problems mentioned in subsection 2.1.4 - reduced qubit count and connectivity; coherent, incoherent and SPAM noise - currently inhibit quantum hardware from being large enough and good enough to implement full algorithms. The prospect of a fault-tolerant era in quantum computing seems distant, if promising. As it was mentioned before, in 2014 the VQE was proposed as the first algorithm inserting a quantum computer in an optimization loop and obtaining a solution via the variational principle. The problem it meant to solve was one familiar to us from the previous section: finding the eigenstates of a Hamiltonian representing a chemical system. Variational algorithms were far from the first approach for using quantum resources to solve the Schrödinger equation. In fact, the problem is strongly tied with the idea of using quantum computers for the simulation of quantum mechanical systems as proposed by Feynman. Because of the exponential growth of the classical resources required to describe a molecule as its size grows, the task at hand is an example of something that quantum computers were created to tackle. And by the time quantum computers were still in the conceptual realm, algorithms were certainly not developed with noise-resilience or qubit limitations in mind - rather provable speedups. If fault-tolerant quantum computers were available to us now, the problem could be solved by means of the QPE algorithm. QPE allows obtaining the eigenvalue of a unitary operator corresponding to an eigenstate given as input. By applying QPE to the (evidently unitary) time evolution operator 𝑒𝑖𝐻𝑡 , it is possible to find the energy eigenstates. The only requirement is that we have a way of preparing states with high overlap with the eigenstates, which can be met by resorting to adiabatic state preparation (ASP). This procedure of using the QPE routine in combination with ASP was detailed in [94]. 33 CHAPTER 2. THEORETICAL FRAMEWORK However, this approach does not seem feasible in the NISQ era: employing the QPE routine to solve relevant problems would require the application of millions to billions of quantum gates [67], implying a circuit depth beyond what can be endured by near-term quantum computers (given that they’re not error-corrected, and coherence times are still limiting). This is where variational hybrid quantum-classical algorithms come in. They resort to the help of a classical optimizer to spare the quantum computer from bearing the load of the algorithm in full. The goal is to provide a near-term solution for problems that are classically hard to solve but, hopefully, easy for quantum computers. Another good example is the Quantum Approximate Optimization Algorithm (QAOA), a variational algorithm aiming to solve combinatorial optimization problems introduced in [24]. Among the problems suitable for QAOA are the Quadratic Unconstrained Binary Optimization (QUBO) problems, of which MaxCut is a famous example with several real-life applications (ranging from physics to communication networks and circuit design). What is interesting is that an infinitely slow adiabatic evolution leading to the solution of the problem is a special case of the QAOA circuit, in the limit of infinitely many layers. If we allow the QAOA ansatz to be big enough, there is a performance guarantee based on the adiabatic theorem: the ansatz is sure to contain the solution. The QAOA both parallels and contrast with the Quantum Adiabatic Algorithm (QAA), which provides a more solid solution by relying in full in the adiabatic theorem. QAOA reduces the circuit depth under what would be necessary for a truly adiabatic evolution, dispenses with the assurance of the adiabatic theorem, and inserts the circuit in a classical optimization loop, converting some gate parameters to variational parameters. The hope is that the extra variational freedom may make up for what is lost by using a shallower circuit. While the parallel with QAA seems like a good reason to believe the QAOA optimization might be successful, if one employs a lower number of layers than required by the adiabatic theorem, the theorem asserts nothing. So against QAA, an algorithm with solid foundations but demanding the unattainable in the NISQ era, we have QAOA: a NISQ-friendly algorithm of tentative construction. The priority shift from performance guarantees offered under the assumption of a fault-tolerant scenario to near-term viability is the bedrock of variational quantum algorithms. 2.3.2 Outline As VQAs might differ significantly in some key aspects and cover a wide variety of applications, quantum algorithms belonging to this class exhibit great structural affinity among themselves. A high-level, greatly simplified scheme of the outline of VQAs is presented in figure 6. The basic idea of VQAs is to use a classical computer to optimize an expectation value measured in a quantum computer. In the following, the four steps labeled in the scheme will be explained in more detail, assuming the quantum circuit model. 34 2.3. VARIATIONAL QUANTUM ALGORITHMS Figure 6: Typical outline of a variational quantum algorithm. 1 - Preparation of the Reference State The quantum part of the algorithm starts with the preparation of a reference state 𝜓𝑟𝑒𝑓 . This is the first portion of the ansatz , and it will typically be a shallow (often constant-depth) circuit with fixed gates and no variational freedom. The role of this circuit is to start the computation in a state other than the typical computational basis state |0...0i. In the specific case of chemistry applications, this step might be of great importance in assuring that the computation starts in a state belonging to the correct symmetry subspace dictated by the problem. For example, it might be desired that the initial state has the correct particle number, and the correct total spin and 𝑍spin projection expectation values. Of course, it isn’t always necessarily the case that a reference state different from the all-zero state is desired; when such occurs, the reference state preparation step is simply omitted. 2 - Application of the Parameterized Unitary After the reference state has been prepared, the parameterized portion of the ansatz circuit follows. This is the part that contains the variational freedom. This is represented in the scheme by the parameterized unitary 𝑈(® 𝜃). This notation simply means that, rather than a fully-fixed circuit, a parameterized circuit is applied in this step. The parameters are considered to be organized into a parameter vector ® 𝜃. After application of the ansatz, we are left with the parameterized state 𝑈(® 𝜃)𝜓𝑟𝑒𝑓 , which will differ among iterations (via the dependence on the parameter vector, updated each round). The complete ansatz then consists of the combination of the circuit that prepares the reference state 35 CHAPTER 2. THEORETICAL FRAMEWORK 𝜎+ 𝑖|0i𝑖=|1i𝑖𝜎+ 𝑖|1i𝑖=0 𝜎− 𝑖|1i𝑖=|0i𝑖𝜎− 𝑖|0i𝑖=0 (2.56) For a single fermionic mode, identifying 𝜎+ 𝑖,𝜎− 𝑖with 𝑎† 𝑖,𝑎𝑖respectively would suffice. However, the antisymmetry requirement would not be satisfied, so this simpler mapping would not be fitting in a multi-mode scenario. That is the purpose of the phase that appears in 2.54. The operator 𝑒𝑖𝜋 Í𝑖−1 𝑘=1𝜎+ 𝑘𝜎𝑘is typically called the Jordan-Wigner string operator . It can be readily seen that it adds a phase 𝑒𝑖𝜋𝑁𝑘<𝑖to the operators acting on qubit 𝑖, where 𝑁𝑘<𝑖is the number of the occupied fermionic modes 𝑘for which 𝑘<𝑖. This is crucial for the anti-symmetry requirement, and it is where the previously mentioned notion of order comes in. We have to establish and be consistent with a sequence, i.e. choose an order ... <𝑖<𝑗<... relating all fermionic modes. It is easy to see what is the effect of the Jordan-Wigner string by considering the successive creation of fermions in two different unoccupied modes. When the first fermion is added to the mode with the lowest index, it affects the string operator of the operator that creates the second fermion, but when the first fermion is added to the mode with the highest index, the string operator of the operator that creates the second fermion is left unchanged. This causes the wave function in these cases to differ by a phase factor 𝑒𝑖𝜋 =−1, as antisymmetry requires: whenever we swap the order of creation of fermions in two modes, we get a flipped sign. It is convenient to write the Jordan-Wigner transform 2.54 using Pauli operators. This is easily done using 2.55 and noting that the effect of the Jordan-Wigner string on operators acting on mode 𝑖amounts to adding a minus sign when the parity of the occupation numbers of the modes under 𝑖is odd. This parity calculation can be done through a string of Pauli 𝑍operators. 𝑎† 𝑖→1 2 𝑖−1 Ö 𝑘=1 𝑍𝑘· (𝑋𝑖−𝑖𝑌𝑖) 𝑎𝑖→1 2 𝑖−1 Ö 𝑘=1 𝑍𝑘· (𝑋𝑖+𝑖𝑌𝑖) (2.57) The formulation of the Jordan-Wigner transform in 2.57 is very convenient for quantum computing applications. Finally, it is worth noting that this definition of the Jordan-Wigner mapping is not unique. Often the roles of the operators 𝜎+ 𝑖and 𝜎− 𝑖in 2.54 are reversed, causing the empty and occupied fermionic modes to be represented by qubit states |1iand |0irespectively (the opposite of what the definitions presented here lead to). 42 2.3. VARIATIONAL QUANTUM ALGORITHMS (a) OpenFermion [2] orbital ordering. (b) Qiskit [3] orbital ordering. Figure 7: Circuit for preparing the Hartree-Fock reference state, for the 𝐻2molecule. The circuits were drawn in the IBM Quantum Composer [41]. 2.3.3.2 Preparation of the Reference State The reference state used in VQE is typically the Hartree-Fock approximation to the ground state. Because the computational basis states usually represent the Slater determinants obtained as approximate eigenstates of the system in performing the self-consistent mean field Hartree-Fock calculations, the reference state preparation state simply amounts to a O(1)circuit with 𝑁𝑒parallel Pauli 𝑋gates, where 𝑁𝑒is the number of electrons. Preparing the reference state amounts to applying Pauli 𝑋gates on the qubits that, under the JordanWigner mapping (see 2.3.3.1), correspond to occupied spin-orbitals. The exact circuit depends on the convention for organizing the spin-orbitals. Figure 7contrasts the reference state preparation for the hydrogen molecule in the two software packages used for creating the circuits: OpenFermion [2], used in conjunction with CIRQ [22], and Qiskit [3]. In the basis set that was always used in the simulations (the STO-3G minimal basis set), 𝐻2has four spin-orbitals, so that the electronic states are represented by four qubits. OpenFermion alternates between 𝛼(up) and 𝛽(down) type orbitals. With this ordering, each spatial orbital is represented by two qubits that are adjacent on the register. Figure 7a shows how to create the state |0011i, assuming little endian ordering. Under OpenFermion’s convention, this corresponds to creating two fermions on the two spin-orbitals associated with the lowest-energy spatial orbital (the orbitals are assumed to be organized in ascending energy from right to left). In contrast, Qiskit uses block-wise orbital ordering: all 𝛼orbitals come first, followed by all 𝛽orbitals. This means that the first qubit of the register corresponds to the same spatial orbital as the one in the middle of the register, with the two representing opposite spins. Figure 7b shows how to create the state |1010i(assuming little endian ordering once again). In spite of the differences in convention, the states 43 CHAPTER 2. THEORETICAL FRAMEWORK prepared by circuits 7a and 7b represent the same fermionic state. Evidently, the ansatz and Hamiltonian averaging must be done in accordance with the chosen convention. 2.3.3.3 Choice of Ansatz As was discussed, the typical ansätze for VQE are chemistry-inspired. Even though they may not lead to the shallowest circuits, the fact that they’re tailored to chemistry problems in particular helps create trial states in the appropriate symmetry subspace. As a consequence, the optimization is helped by the fact that subspaces where the ground state is sure not to be are entirely avoided. Ansätze that don’t leverage previous knowledge about the system typically lead to barren plateaus. These trainability issues have been shown to be a consequence of the choice of ansatz that cannot be helped through the choice of optimizer [60,10,39]. UCC theory enjoys great popularity due to its natural genesis in the context of previous computational chemistry methods. UCC proposals are also fairly broad, and this class has been expanding as a consequence of research efforts to make the ansätze more compact and more accurate. The article [80] explored singlet and pair Quantum Unitary Coupled Cluster Doubles (q-UCCD) approaches (where the q stands for quantum ), which reduce the double excitations appearing in the ansatz to a smaller subset. Singlet q-UCCD splits the double excitations into singlet and triplet components, and makes use of symmetry / anti-symmetry to reduce the number of excitation operators. Pair q-UCCD limits the excitations to those which excite a pair of electrons with opposite spins from one spacial orbital to another. Both of these ansätze were tested in several systems with results of compelling accuracy, while reducing the number of two-qubit gates in the ansatz circuit by up to 75%. A different ansatz termed k-Unitary Paired Coupled Cluster Generalized Singles and Doubles (UpCCGSD) was introduced in [53]. k-UpCCGSD consists of k products of unitary paired coupled-cluster double excitations, as well as all generalized single excitations. Here generalized indicates that virtual-to-virtual and occupied-to-occupied excitations are also allowed (orbitals are considered virtual or occupied if they are so in the reference state). The K-UpCCGSD ansatz increases the variational freedom as compared to the simple pair Unitary Coupled Cluster Generalized Singles and Doubles (UCCGSD) (pUCCGSD) ansatz, while still enjoying a linear scaling of the circuit depth with the system size. The k-UpCCGSD was shown to be capable of achieving chemical accuracy with a better resource scaling than UCCSD and UCCGSD (the generalized version of the former), both of which require an O(𝑁4)circuit depth, where N is the number of orbitals. Changing the form of the ansatz is not the only option. Low-rank decomposition of the basic UCC operator was used in [62] to reduce the gate complexity of each Trotter step from O(𝑁4)to O(𝑁3). By exploiting the structure and symmetries of the cluster operator and truncating small terms, the depth of the circuits can be reduced without precluding chemical accuracy. 44 2.3. VARIATIONAL QUANTUM ALGORITHMS 2.3.3.4 Hamiltonian Averaging The part of the algorithm that consists of obtaining the expectation value of the cost Hamiltonian (in this case a physical Hamiltonian) is typically called Hamiltonian averaging or quantum expectation estimation . In subsection 2.1.2, a procedure to measure arbitrary observables was outlined; it consisted of decomposing them as a linear combination of Pauli strings (formula 2.14) and calculating the expectation value as a weighed average of the corresponding expectation values (formula 2.15). The expectation value of the Hamiltonian can be obtained exactly this way; the question is whether this can be done efficiently. Luckily, it is often the case that the Hamiltonian of a physical system can be decomposed as a linear combination of a polynomial number of Pauli terms. This includes the electronic structure Hamiltonian (formula 2.30), as well as the Hamiltonians of the Ising and Heisenberg models [67]. The problems addressed in this dissertation concern electronic Hamiltonians with a number of terms that scales like O(𝑁4)on the number of spin-orbitals/qubits 𝑁[75]. An example of such a Hamiltonian, for the very simple case of the hydrogen molecule, can be found in appendix B. While the scaling is polynomial, obtaining the expectation value can still imply a large number of measurements, especially in chemistry problems which typically require exceptional accuracy. One way this can be improved is by organizing the Pauli strings appearing in the decomposition of the operator into commuting groups, that can be measured simultaneously. As a simple two-qubit example, we can consider the following set of observables. 𝐼𝐼, 𝑍𝐼, 𝐼𝑍, 𝑍𝑍 The expectation values of all four terms can be obtained from the same circuit. There are multiple ways of arranging Pauli strings into commuting groups; for example, if we have the set of observables 𝑍𝐼, 𝐼𝑍, 𝑍𝑍, 𝑋𝐼, 𝐼𝑋, 𝑋𝑋, 𝑋𝑍, 𝑍𝑋, we can chose to split the Pauli strings into groups as {𝑍𝐼, 𝐼𝑍, 𝑍𝑍 };{𝑋𝐼, 𝐼𝑋, 𝑋𝑋 };{𝑍𝑋 };{𝑋𝑍 }, but also as {𝑍𝐼, 𝐼𝑋, 𝑍𝑋 };{𝑋𝐼, 𝐼𝑍, 𝑋𝑍};{𝑍𝑍 };{𝑋𝑋 }. These divisions only exploit qubit-wise commutativity. However, two Pauli strings also commute if they have non-commuting Pauli operators in an even number of indices. This general commutativity allows us to further reduce the circuit repetitions that are required to obtain the final expectation value. In our example, it brings them down to 3 8as compared to when no grouping is used; we can organize the operators as {𝑍𝐼, 𝐼𝑋, 𝑍𝑋 };{𝑋𝐼, 𝐼𝑍, 𝑋𝑍};{𝑍𝑍, 𝑋𝑋 }. 45 CHAPTER 2. THEORETICAL FRAMEWORK Even though this is evidently not the only possibility; we can also have {𝑍𝐼, 𝐼𝑍, 𝑍𝑍 };{𝑋𝐼, 𝐼𝑋, 𝑋𝑋 };{𝑍𝑋, 𝑋𝑍 }. In this example, both options reduce the necessary circuit repetitions by the same amount (assuming these repetitions are to be distributed equally). However, this is not always necessarily the case. What is more, the task of finding the optimal partitioning becomes harder the larger the Hamiltonian is, and taking into account sampling noise complicates the task further. Proposals of grouping strategies can be found in references [59,93,49,72,45,47,65,46,97,28,29,89,19,40]. The last article among these proposed an approach that allows for a cubic reduction in term groupings over previous proposals, with the drawback of requiring the execution of a linear-depth circuit before the measurement. A table comparing all these strategies can be found in this same reference. 2.3.3.5 Estimating Shot Requirements The number of terms in the Hamiltonian decomposition, and the number of commuting groups, are not the only factors weighing in the shot requirements of VQE. The impact of sampling noise in the error, which typically also increases with the size of the system, is another relevant aspect. The following estimation of shot requirements was based on reference [93]. Assuming that the Hamiltonian operator has been decomposed in the form of formula 2.14 and that for each Pauli string 𝑃𝑖we use 𝑀𝑖shots, the error coming from the finite number of shots in the expectation estimation can be written 𝜖2=Õ 𝑖 |ℎ𝑖|2𝑉𝑎𝑟 (ˆ 𝑃𝑖) 𝑀𝑖 . The variances depend, not only on the Pauli operators appearing in the Hamiltonian, but also on the particular state. We can bound them as 𝑉𝑎𝑟 (𝑂𝑖) ≤ 1, since any Pauli operator has eigenvalues ±1and so will their product. Without further knowledge about their variance, the best choice is to distribute the shots by the operators in proportion to the norm of their coefficients: 𝑀𝑖∝ |ℎ𝑖|. Assuming this distribution, one can estimate the necessary number of shots to hit a precision 𝜖as 𝑀≈(Í𝑖|ℎ𝑖|)2 𝜖2. The term (Í𝑖|ℎ𝑖|)2will depend strongly on the molecule in study: it will be larger for larger molecules, requiring a larger number of shots for the same precision if the system size is increased. Given that chemical accuracy is 1kcal/mol ≈1.59 ×10−3Hartree, achieving it requires 𝑀≈0.4×106× (Õ 𝑖|ℎ𝑖|)2 46 2.3. VARIATIONAL QUANTUM ALGORITHMS shots. This is already of the order of 107to 108for small molecules like helonium and lithium hydride; 𝐹𝑒2𝑆2, for instance, would require around 1013 shots. This number of shots needs to be performed every time one needs to evaluate the energy, which happens many times throughout the optimization. The shot requirements are heavier for larger molecules, since there are more Pauli strings in the Hamiltonian, larger coefficients |ℎ𝑖|, and more parameters to optimize (resulting in more iterations and more function evaluations per iteration being required). Further, this shot count doesn’t guarantee proper functioning of the optimizer: it ensures that the difference between the energy estimate and the true energy in a state are within chemical accuracy, but the final state isn’t necessarily close to the ground state. The error might still be enough to confuse the optimizer into outputting a wrong guess for the ground state (e.g. by tricking it into an early convergence). 47 3 Static Ansätze for VQE This chapter aims to present the results of VQE simulations with static ansätze. In the first section, an entirely problem-agnostic ansatz will be applied in searching for the ground state of the helium hydride ion. The impact of the optimizer, optimization hyperparameters, and noise on performance will be assessed. In the second section, a problem-inspired ansatz will be applied in searching for the ground state of the hydrogen molecule. The impact of several types of noise on performance will be analysed. For creating and executing quantum circuits, CIRQ [22] and Qiskit [3] were used in alternation; the CIRQ simulator was used as well as IBMQ’s simulators and real quantum computers. Which package and backend was used in obtaining which results will be clear along the chapter. When analysing the results, the final VQE energy and state will be compared with those obtained from performing exact diagonalization on the FCI Hamiltonian. The FCI solution is exact up to the BornOppenheimer approximation, relativistic effects, and the chosen basis set. 3.1 Problem Agnostic Ansätze As was explained before, problem agnostic ansätze come with several advantages: they don’t demand previous knowledge about the problem or system, nor a way of leveraging such knowledge. Additionally, taking the focus away from the problem allows for a more hardware-friendly circuit implementation, which may decrease the depth of the circuits. The main disadvantage of this type of ansatz is that randomly initialized parameterized quantum circuits often lead to barren plateaus in the optimization landscape [60]. This section will present simple a example of a problem agnostic ansatz, applied to the helium hydride ion. An ansatz spanning the whole Hilbert space was created and employed in the simulations; such an approach is feasible given the dimensionality of the problem. The section contains results from running the algorithm on the CIRQ simulator and on IBMQ’s backends (simulators or real quantum processors). 3.1.1 Application to Helium Hydride The purpose of the following is to implement VQE as in the article that originally introduced it [67], which applied the algorithm to the helium hydride ion using a photonic quantum processor. In addition to the CIRQ and Qiskit implementations, an independent noise-free version of VQE, based on matrix algebra, was implemented to verify correctness. The Hamiltonians for the different bond distances were taken from the article (which in turn obtained them via the PSI3 computational package [20]), as was the parameterization of an arbitrary 2-qubit state. 48 3.1. PROBLEM AGNOSTIC ANSÄTZE The ansatz employed here searches through the whole Hilbert space, much like the ansatz in the original article. However, while in the article the parameterization was rewritten into a convenient form that directly dictated the values of six phase shifters in the photonic chip, here the state preparation was abstracted of the specifics of any physical implementation and instead mapped to a generic circuit. 3.1.1.1 Preparing an Arbitrary 2 Qubit State An arbitrary (pure) two qubit state can be written as in equation 3.1. |𝜓i=𝛼|00i+𝛽|01i+𝛾|10i+𝛿|11i(3.1) The computational basis coordinates 𝛼,𝛽,𝛾and 𝛿are complex numbers that fully determine the state. Since optimizers typically deal with real variables, it is convenient to rewrite the state using real parameters. Further, the degrees of freedom corresponding to the normalization and the global phase should be removed. With this, the two qubit state can be described by six real variational parameters. A possible parameterization is presented in 3.2. |𝜓i=cos 𝜃0 2cos 𝜃1 2|00i+cos 𝜃0 2sin 𝜃1 2𝑒𝑖𝜔1|01i+ sin 𝜃0 2𝑒𝑖𝜔0cos 𝜃2 2|10i+sin 𝜃0 2𝑒𝑖𝜔0sin 𝜃2 2𝑒𝑖𝜔2|11i (3.2) Here, 𝜃𝑖∈ [0, 𝜋]and 𝜔𝑖∈ [0,2𝜋]. This was the parameterization used in the original VQE article [67], and it will be used here as well. We wish to transform this parameterized description of the state into a parameterized circuit that we can use as our ansatz. The concept of Schmidt decomposition, presented in theorem 2.7 of reference [64], is useful here. If we have a pure state of a composite system 𝐴𝐵,𝜓𝐴,𝐵, then the Schmidt decomposition allows us to write it as 𝜓𝐴,𝐵=Õ 𝑖 𝜆𝑖|𝑖𝐴i|𝑖𝐵i,(3.3) where the |𝑖𝐴i,|𝑖𝐵iare orthonormal states of the systems 𝐴,𝐵respectively. They form the Schmidt basis of the respective system. The 𝜆𝑖, called Schmidt coefficients , are real non-negative numbers that satisfy Í𝑖𝜆2 𝑖=1. We can find the Schmidt decomposition of a two qubit state written in the computational basis (formula 3.1) by using Singular Value Decomposition (SVD). Firstly, we must organize the coordinates into a two-bytwo matrix 𝑀. 𝑀= 𝛼 𝛽 𝛾 𝛿!(3.4) 49 CHAPTER 3. STATIC ANSÄTZE FOR VQE This is done by identifying the coordinate of computational basis state |𝑖𝑗iwith the matrix entry 𝑀𝑖,𝑗 : the left qubit, 𝐴, is associated with row indices, and the right qubit, 𝐵, is associated with column indices. Calculating the scalar product of the state of the composite system with a state of system 𝐴or 𝐵(|𝜓𝐴i or |𝜓𝐵i) corresponds to leftor right-multiplying matrix 𝑀by the bra h𝜓𝐴|or the ket |𝜓𝐵i, which results in a linear superposition of the rows or columns of matrix 𝑀, respectively. Calculating the scalar product of the state of the composite system 𝐴𝐵 with a computational basis state of the same system amounts to doing both simultaneously, i.e. 𝑖𝑗𝜓𝐴,𝐵≡h𝑖| 𝛼 𝛽 𝛾 𝛿!|𝑗i. We can now use SVD decomposition to write 𝑀in the more convenient way of 3.5. 𝑀=𝑈 𝜆00 0𝜆1!𝑉(3.5) The 𝜆𝑖are the singular values of matrix 𝑀, non-negative real numbers that we can easily identify with the Schmidt coefficients of formula 3.3.𝑈and 𝑉𝑇(note the transpose) are two-by-two unitary matrices, whose column vectors form the Schmidt basis for 𝐴and 𝐵respectively. It will be necessary to create a circuit to implement 𝑈and 𝑉𝑇. A generic single-qubit unitary can be written by means of three rotations, along with a global phase. This is done in 3.7. 𝑈=𝑒𝑖𝜃0𝑅𝑍(𝜃1)𝑅𝑌(𝜃2)𝑅𝑍(𝜃3)=𝑒𝑖𝜃0𝑒−𝑖𝜃1𝑍 2𝑒−𝑖𝜃2𝑌 2𝑒−𝑖𝜃3𝑍 2(3.6) =𝑒𝑖𝜃0©« 𝑒−𝑖(𝜃1+𝜃3) 2·cos𝜃2 2−𝑒−𝑖(𝜃1−𝜃3) 2·sin𝜃2 2 𝑒−𝑖(−𝜃1+𝜃3) 2·sin𝜃2 2𝑒−𝑖(−𝜃1−𝜃3) 2·cos𝜃2 2ª®¬(3.7) The rotation angles can easily be found from this; for example, they can be obtained using the expressions in 3.8, where ∠𝑈𝑖𝑗 denotes the phase 𝜙of the complex number 𝜌𝑒𝑖𝜙 in line i and column j of matrix 𝑈, and |𝑈𝑖 𝑗 |denotes its modulus 𝜌. 𝜃0=∠𝑈00 +∠𝑈11 2𝜃1=−∠𝑈00 +∠𝑈10 𝜃2=2 arccos |𝑈00|𝜃3=−∠𝑈00 +∠(−𝑈10) (3.8) It should be noted that 𝜃0corresponds to a global phase with no relevance. The only parameters that are important and should be found are 𝜃1,𝜃2and 𝜃3. A procedure to create a circuit that prepares a state from its four complex coordinates in the computational basis can now be easily described. Once 𝑈,𝑉,𝜆0and 𝜆1have been found by SVD decomposition (allowing us to rewrite the coordinate matrix as in 3.5), it can be divided into the three following steps: 50 3.1. PROBLEM AGNOSTIC ANSÄTZE 1. Create the state 𝜆0|00i+𝜆1|10i(from |00i). This can be done by rotating qubit 𝐴around the 𝑌 axis by an angle of 2 arccos |𝜆0|. Because 𝜆2 0+𝜆2 1=1, this is the same angle as 2 arcsin |𝜆1|. 2. Create the state 𝜆0|00i+𝜆1|11i(from 𝜆0|00i+𝜆1|10i). This can be done by applying a CNOT, with qubit 𝐴as the control and qubit 𝐵as the target. 3. Create the state 𝑈⊗𝑉𝑇(𝜆0|00i+𝜆1|11i)(from 𝜆0|00i+𝜆1|11i). This can be done by decomposing the single qubit unitaries 𝑈and 𝑉𝑇into single qubit rotations, and applying them to qubits 𝐴and 𝐵respectively. Steps 1 and 2 create the state 𝜆0|00i+𝜆1|11i, which is the desired state in the Schmidt basis . Step 3 rotates each qubit from its Schmidt basis to the computational basis, creating the final state. The corresponding unitaries are those formed by the Schmidt vectors of the respective qubits. The outlined procedure leads to the state preparation circuit presented in figure 8. Such an ansatz has 7 variational parameters: one from the rotation on the first qubit in step 1, and 3 from the decomposition of each of the unitaries in step 3. Figure 8: Circuit with 7 variational parameters for preparing an arbitrary 2-qubit state, drawn in the IBM Quantum Composer [41]. One thing that stands out is that there were six variational parameters in the parameterization 3.2. In fact, there are six degrees of freedom left in a two-qubit state once we remove those corresponding to the normalization and the global phase. This means that there is a redundant degree of freedom in our parameterized circuit. Luckily, this can easily be solved by noting that after the two first steps (i.e. after the CNOT gate), we have a symmetric state 𝜆0|00i+𝜆1|11i. This means that the first two parallel Z rotations, with parameters 𝑝2and 𝑝5in figure 8, can be merged into a single Z rotation with angle 𝑝2+𝑝5on either of the qubits, with exactly the same effect. The final 6-parameter ansatz is presented in figure 9, where the parameter labels were redefined to run from 𝑝1to 𝑝6. 51 CHAPTER 3. STATIC ANSÄTZE FOR VQE (a) Statevector simulator (b) QASM Simulator (c) IBMQ Lima Figure 13: Evolution of the VQE optimization for 𝐻𝑒𝐻+at an interatomic distance of 90pm, using three different IBM Quantum [42] backends: the statevector simulator (with no noise of any kind), the QASM simulator (with sampling noise) and the Lima device (a 5-qubit Falcon processor). In the last two cases, 1024 shots were used. The energy is plotted in blue and the overlap with the true ground state in red; the pale blue line marks the true ground energy. The same random starting point was used for all backends. The chosen optimizer was COBYLA. 58 3.1. PROBLEM AGNOSTIC ANSÄTZE (a) 256 shots (b) 1024 shots (c) 4096 shots (d) 8192 shots Figure 14: Evolution of the VQE optimization for 𝐻𝑒𝐻+at an interatomic distance of 90pm, for different shot counts and using the Lima device. The remaining conditions (starting point, optimizer) are the same as in figure 13. 59 CHAPTER 3. STATIC ANSÄTZE FOR VQE accuracy. 60 3.2. PROBLEM TAILORED ANSÄTZE 3.2 Problem Tailored Ansätze Leveraging problem-specific knowledge to design the ansatz allows narrowing the search space by excluding regions that are not expected to contain the solution. Such techniques reduce the expressibility of the ansatz without abdicating precision, because the solution is still contained within the variational form. Since problem-tailoring is of great importance in preventing barren plateaus, and this trainability issue renders variational quantum algorithms inefficient, a great part of the research surrounding VQAs is aimed at developing precise and shallow ansätze of this type. This section contains the results obtained from applying VQE to molecular hydrogen (𝐻2) both in simulators and real quantum computers., using a problem-tailored ansatz. The employed ansatz was UCCSD, introduced in section 2.2. As was explained then, this variational form arises in the context of conventional CC theory as a natural modification to the non-variational wave form of the latter, which in turn is the natural size-consistent successor to the Configuration Interaction variational form. While UCCSD is not classically tractable, the unitarity that forces variationality allows for convenient implementation in quantum circuits. Due to being unitary and having solid roots in older computational chemistry methods, the UCCSD ansatz has become a popular choice in VQAs addressing chemistry applications. 3.2.1 Application to Molecular Hydrogen The UCCSD-VQE algorithm was implemented in this project using CIRQ [22] and Qiskit [3], alternative software libraries for manipulating and running quantum circuits. Tasks such as obtaining the Hamiltonians and performing fermion-to-qubit mappings were done via OpenFermion [2] in the first case and Qiskit’s chemistry module in the second. In both, the underlying electronic structure information was extracted from PySCF [5]. Additionally, an independent noise-free UCCSD-VQE based on matrix algebra was implemented to verify correctness. Simulations were done to analyse the effect of sampling noise on UCCSD-VQE, as well as of CNOT gate errors, and of thermal relaxation and SPAM errors, using noise models created in Qiskit. Density matrix simulations were employed to analyse the purity of the final state, using noise models created from the specifications of IBMQ [42] quantum processors to mimic their behaviour. 3.2.1.1 Bond Dissociation Graph Much like in the previous section, we will start with the bond dissociation plots to get a general sense of performance before delving into the evolution of individual optimizations. In figure 15, we can see the VQE energies tracing out the bond dissociation curve of 𝐻2. Here, the simulations include no noise other than sampling. The effect of the number of shots clearly shows in the results, with the VQE energy getting closer and closer to the exact diagonalization curve from figure 15a to 15d. In the limit of infinite shots, the VQE energy essentially matches the FCI result (the observed error was of the order of 10−9). 61 CHAPTER 3. STATIC ANSÄTZE FOR VQE (a) QASM simulator, 256 shots (b) QASM simulator, 1024 shots (c) QASM simulator, 8192 shots (d) Statevector simulator Figure 15: Single-run UCCSD-VQE ground energy along the bond dissociation curve of the hydrogen molecule, obtained using IBMQ [42] simulators. The results plotted in figures 15a-15c, ordered by increasing shot count, were obtained using the QASM simulator. The results plotted in figure 15d were obtained using the state vector simulator (equivalent to infinite shots). The red and green curves represent the FCI energy and the Hartree-Fock energy, respectively. 62 3.2. PROBLEM TAILORED ANSÄTZE It is interesting to compare the VQE energy with that obtained by performing the Hartree-Fock selfconsistent field calculations. The state resulting from these calculations is the reference state in VQE, upon which the UCCSD ansatz performs a ‘correction’. It becomes evident that, in the presence of sufficiently weak noise, optimizing the UCCSD variational form is enough to bring the energy away from the HartreeFock approximation and into the FCI value. The variational freedom takes care of electronic correlation effects that go unaccounted for in the Hartree-Fock approximation. (a) QASM simulator, 256 shots (b) QASM simulator, 1024 shots (c) QASM simulator, 8192 shots (d) Statevector simulator Figure 16: Bond dissociation curves of figure 15, zoomed in on the region nearing the equilibrium bond length. The blue area marks the region of chemical accuracy (error of less than 1kcal/mol). To better analyse the results, figure 16 shows the same plots as figure 15, now magnified around the equilibrium bond length (the minimum of the curve). In these plots, it is visible that reaching chemical accuracy requires a considerable number of shots. With only 256 shots (figure 16a), UCCSD-VQE is often 63 CHAPTER 3. STATIC ANSÄTZE FOR VQE outperformed by the Hartree-Fock approach; when this happens, we know that the UCCSD correction didn’t at all improve upon the reference state, and the optimization was entirely in vain. 1024 shots (figure 16b) seem to be enough for UCCSD-VQE to consistently beat the Hartree-Fock method, but the single-run energy is still far from reaching chemical accuracy. In figure 16c, we can see that with 8192 shots, the energy is already often inside or nearing the region of chemical accuracy, even with a single run (taking the median over several runs would further refine the result). Of course, in the limit of infinite shots (figure 16d), the performance is nearly perfect. 3.2.1.2 Running the Algorithm on Cloud Quantum Computers So far, the simulations have included at most sampling noise - but they all assumed ideal behaviour from the quantum processor. In the previous section, we saw that the performance of the algorithm in real backends could prove significantly worse because of the other types of noise they inevitably suffer from. With UCCSD the difference becomes even more relevant, because the circuit is deeper: the 2-qubit ansatz from the previous section had a single CNOT and a total of only 7 gates. Even for 𝐻2, the simplest molecule possible, implementing the UCCSD operator requires significantly longer circuits. Figure 17 showcases the evolution of the VQE optimization for a single point (interatomic distance) in three different scenarios: noise-free (figure 17a), sampling noise only (figure 17b), real quantum processor (figure 17c). Here, the optimization of a single run was plotted instead of the final results of multiple runs (as in the bond dissociation graphs) because of the time required to run hybrid algorithms on cloud quantum computers with fair-share queuing. The difference between the noise-free and sampling noise only scenarios is already remarkable. While in the former the energy is fast to stabilize inside the region of chemical accuracy, in the latter the energy oscillates enough that it’s in and out of the region throughout the whole optimization. Even towards the last iterations, the energy does not meet chemical accuracy consistently: it is outside the region in the second-to-last iteration. However, if sampling noise is significant, it does not come close to the impact of other types of noise. Figure 17c shows that running the algorithm in IBMQ’s Belem backend, a real quantum processor, does not allow recovering any valid results. The VQE output energy oscillates about 1 Hartree over the true ground energy, without being remotely close to reaching chemical accuracy (of the order of the millesimals of Hartree). This is a sign that quantum information has been washed out by noise. The UCCSD ansatz for 𝐻2, obtained using Qiskit’s variational form and transpiled onto Belem, has 130 CNOTs; the circuit depth is approximately proportional to the CNOT count. The results show that the circuit in question is too deep for any meaningful data to be recovered by the end. An approximate density matrix representing the state at the end of the noisy UCCSD ansatz was obtained by simulating the algorithm on the QASM simulator, now with a noise model aiming to mimic the behaviour of the Belem processor. The noise model was obtained from Qiskit [3] and is built from data concerning the specific backend. This procedure allowed obtaining estimates to the purity of the final state 64 3.2. PROBLEM TAILORED ANSÄTZE (a) Statevector simulator (b) QASM Simulator (c) IBMQ Belem Figure 17: Evolution of the UCCSD-VQE optimization for 𝐻2at an interatomic distance of 0.74Å, using three different IBM Quantum [42] backends: the statevector simulator (with no noise of any kind), the QASM simulator (with sampling noise) and the Belem device (a 5-qubit Falcon processor). In the last two cases, 8192 shots were used. The energy is plotted in blue; the pale blue region marks the area of chemical accuracy. The same starting point, the Hartree-Fock reference state, was used for all backends. The chosen optimizer was COBYLA. 65 CHAPTER 3. STATIC ANSÄTZE FOR VQE and its fidelity with the ideal noiseless output. The fidelity of the final state with the ideal state was found to be only 0.081, in contrast with its fidelity with the fully mixed state, which was 0.40. The final purity of the state was found to be around 0.16. The STO-3G basis set was used; in a minimal basis set (as is the case), 𝐻2has four spin-orbitals, so that fermionic states of the molecule were mapped into states of four qubits. This implies a Hilbert space of dimension 16, so that the purity is lower bounded by 0.0625. The expectation value of the energy calculated in the fully mixed state was found to be -0.097 Hartree. In figure 17c, we can see that the energy of the state at the end of the UCCSD circuit in the Belem backend remains close to that value throughout the whole optimization. The optimizer spends 25 iterations trying to adjust the parameters that will lead to a minimum of the energy, but it is to no avail: the energy shows no more than slight deviations from the mean value, with no sign of approaching the ground energy. (a) QASM Simulator, 256 shots (b) QASM Simulator, 1024 shots (c) QASM Simulator, 8192 shots (d) IBMQ Manhattan, 256 shots (e) IBMQ Manhattan, 1024 shots (f) IBMQ Manhattan, 8192 shots Figure 18: Evolution of the VQE optimization under the same conditions of figure 17 (same molecule, interatomic distance and optimizer), plotted for three different shot counts and two different backends. The QASM simulator (18a-18c) portraits the behaviour of an ideal quantum processor, in contrast with IBMQ Manhattan (18d-18f), a 65-qubit Hummingbird processor. From figure 17, we can already conclude that sampling noise is far from being the limiting factor here. To further illustrate this, figure 18 compares VQE optimizations using the QASM simulator (with sampling noise only) and using the Manhattan backend (another real device) for 256, 1024, and 8192 shots. 66 3.2. PROBLEM TAILORED ANSÄTZE The improvement of the VQE performance with increased shot count is evident when the QASM simulator is used: the energy is visibly faster to stabilize. This is to be expected, since here sampling is the only source of noise. With 256 shots only (figure 18a) the energy shows changes of the order of the centesimals of Hartree up until the last iterations. With 8192 shots (figure 18c) the energy change is of the order of the millesimals of Hartree in every iteration from the seventh to the last one (30). When the Manhattan processor is used, the energy shows the same tendency for stabilization: from 256 shots (figure 18d) to 8192 shots (figure 15c), there is also a decrease in the order of magnitude of the oscillations in the energy towards the last iterations. However, here, this does not translate into any significant improvement in accuracy, because the energy is still very far from the ground energy, with an error of around 1 Hartree. Much like it happened with the Belem backend, the state at the end of the UCCSD circuit in the Manhattan backend is too corrupted for the optimizer to converge to the ground state: the energy evaluations are too far from the exact value. The purity is so low that changing the parameters is barely reflected on the expectation value. We can conclude that the circuits are too deep to be implemented in these quantum computers. 3.2.1.3 Analysing Noisy Instructions To better understand which sources of noise were most affecting the algorithm, the noise model emulating the behaviour of the Belem backend was decomposed into several noise models, each considering only the noise associated with a specific instruction. UCCSD-VQE was then simulated with each of the noise models; the results are plotted in figure 19. The curve that stands out is the one associated with the noise model that isolated CNOT noise. Even though this model ignored all other sources of noise, the performance is remarkably close to that of the complete Belem noise model (figure 17c). In here, the error in the final energy is once again close to 1 Hartree; further, the energy barely moves away from the starting energy throughout the whole optimization. One again, the impact of noise is felt too strongly for any meaningful results to be recovered. It is not surprising that this is the noisiest gate: it is the only two-qubit gate. In this particular backend (Belem), the average CNOT gate time is around 550 nanoseconds, while the single-qubit gate time is typically of the order of tens of nanoseconds. Among the available gates, the CNOT thus corresponds by far to the biggest cost in circuit depth. Since the noise associated with CNOT gates is the one that most impacts UCCSD-VQE, it will be interesting to analyse the relation between the magnitude of CNOT gate errors and the performance of the algorithm. A graphical representation of the impact of the average CNOT gate error on the final UCCSDVQE error can be found in figure 20. To obtain these results, a noise model was created to simulate CNOT gate errors. This model consisted of a two-qubit depolarizing error, succeeded by single qubit thermal relaxation errors on the involved qubits. This is a similar approach to the one employed in Qiskit’s backend noise models, mentioned previously and used in obtaining the results presented in figure 19. In these noise models, made available by the Qiskit 67 CHAPTER 4. PROBLEM TAILORED, DYNAMICALLY CREATED ANSÄTZE: ADAPT-VQE The only thing left to explain is when the algorithm terminates: in order to know when to stop adding operators and growing the ansatz, and finally output an attempted solution, a convergence criterion is necessary. There may be different possibilities, but the one employed in the original article (and reproduced in this project) was to terminate when the total gradient norm falls below a certain threshold 𝜖. At the beginning of each iteration, after calculating all of the derivatives, they are used to calculate this norm. If the norm is still over the threshold, the operator with the highest derivative is added to the ansatz; if not, the algorithm terminates. Evidently, the lower 𝜖is, the stricter the convergence requirement, the larger the final ansatz, and the higher the precision. We can then summarize the ADAPT-VQE algorithm in 7 steps (as labeled in figure 22). 1. Initialize the ansatz at identity, so that the initial state is simply the Hartree-Fock solution. The corresponding state preparation circuit consists of 𝑁parallel Pauli 𝑋gates, where 𝑁is the number of electrons. 2. Measure the gradient of each operator in the pool on the quantum computer. This can be done using formula 4.1. 3. Calculate the square norm of the gradient vector. If it is under a predetermined threshold 𝜖, terminate. If not, proceed to the following step. 4. Using the information obtained in step 2, select the operator with the largest gradient and append it to the ansatz. 5. Initialize the parameter vector with the values from the previous iteration. The newly added variational parameter should be initialized at zero. 6. Re-optimize all the parameters. This is a typical VQE optimization, with the current ansatz. 7. Move on to the next iteration, going back to step 2. It should be noted that before these steps take place, one must have done the necessary classical computations (common to the original VQE) in preparation, to obtain the Hamiltonian operator and the Hartree-Fock spin-orbitals. Further, one must have chosen the constitution of the pool, and the convergence threshold 𝜖. Steps 2 and 6 require preparing the ADAPT-VQE state multiple times, and measuring the expectation value of the Hamiltonian ˆ 𝐻as well as of the commutators ˆ 𝐻, ˆ 𝐴𝑖. A brief explanation on how this can be done follows. In iteration 𝑛, the ADAPT-VQE state is written 𝜓(𝑛)=𝑒𝜃𝑛ˆ 𝐴𝑛. . . 𝑒𝜃1ˆ 𝐴1|𝐻𝐹i. Under the JordanWigner transform 2.57 the excitation operators ˆ 𝐴𝑖are mapped to linear combinations of Pauli strings with imaginary coefficients. The operators 𝑒𝜃𝑖ˆ 𝐴𝑖can then be implemented as explained in subsection 2.1.3. Along with the preparation of the Hartree-Fock state presented in 2.3.3.2, this tells us how to create the state preparation circuit. 74 4.1. THE ADAPT-VQE ALGORITHM The measurements are also straightforward: the Hamiltonian ˆ 𝐻also consists of a linear combination of Pauli strings, so that its expectation value can be obtained using the procedure outlined in subsection 2.1.2 (possibly improved with the strategies for grouping commuting operators mentioned in 2.3.3.4). Since the product of two Pauli strings will also be a Pauli string, measuring ˆ 𝐻, ˆ 𝐴𝑖amounts to defining a new observable that can be measured in a similar way. In all cases, the different strings can be measured in parallel, if multiple quantum processors are available. It is important to note that the pool operators ˆ 𝐴𝑖are basic spin-complemented one and two-body excitations, that are Jordan-Wigner transformed into a linear combination of at most 16 Pauli strings (8 from each of the two anti-Hermitian two-body excitations with opposite spins), independently of the system size. Further, the Hamiltonian itself has been decomposed into a sum of polynomially-many Pauli strings (the number being O(𝑁4)on the number of spin-orbitals 𝑁). The observable corresponding to the commutator ˆ 𝐻, ˆ 𝐴𝑖can then be efficiently measured in a quantum computer. Thus, the extra measurements of the ADAPT-VQE algorithm do not in principle hinder efficiency, as long as the scaling of the size of the operator pool is reasonable. In the case considered in the original article, including only excitations up to order two, this scaling is evidently also quartic on the number of spin-orbitals of the chosen basis set. The measurements per iteration in ADAPT-VQE thus scale like O(𝑁8)(instead of O(𝑁4)like in VQE with static ansätze). It should be noted that, as mentioned previously, modified variants such as ADAPT-V and ADAPT-Vx manage to decrease this measurement overhead at the expense of an acceptable increase in the number of variational parameters and circuit depth [55]. While polynomial, the measurement overhead of ADAPT-VQE as compared to UCCSD-VQE or other predetermined ansätze is considerable. This is the price to pay for growing the ansatz using a procedure this customized: to choose the new operator, the gradients of all operators in the pool must be measured. ADAPT-VQE also implies an overhead in optimizations: instead of a single optimization with a fixed ansatz, we have multiple optimizations with a growing ansatz. The larger number of optimizations implies in itself a larger number of energy evaluations, as well as a bigger burden on the optimizer. All this sacrifice is done in favour of shallower circuits, and in some cases, greater accuracy. The ADAPT-VQE ansatz was designed to be a NISQ-friendly, high-accuracy ansatz for chemistry applications. 4.1.2 Qubit-ADAPT-VQE In 2020, Qubit-ADAPT-VQE was proposed in [85] with the purpose of improving upon the previous ADAPTVQE algorithm. The entirety of the modification lies in the operator pool. Spin-complemented fermionic single and double excitations can consist of a linear combination of up to 18 Pauli strings, and spin-adapted ones of up to 48 Pauli strings. Using the ladder-of-CNOTs implementation of Pauli exponentials [64], the circuit for one of such operators can require already thousands of CNOTs when as few as a dozen spin-orbitals are used to describe a molecule, even assuming all-to-all connectivity. 75 CHAPTER 4. PROBLEM TAILORED, DYNAMICALLY CREATED ANSÄTZE: ADAPT-VQE Qubit-ADAPT-VQE drastically reduces the circuit depth required for the implementation of each operator, by dispensing with the rigorous representation of fermionic excitations. Instead of using a linear combination of Pauli strings that respects fermionic symmetries (e.g. particle number), the operators are broken down into the individual Pauli strings. Because of this, the pool is called qubit pool, whereas the pool with the proper excitations is called fermionic pool. In addition to the one directly arising from the decomposition fermionic excitations, the original Qubit-ADAPT-VQE article introduced a few more pools, which will be covered in the next section. 4.2 Original Pool Choices This section provides a more thorough description of the different pools used in the original articles: as it was stated before, they can be of fermion or qubit inspiration. The former take inspiration from fermionic excitations, and the latter opt for more hardware-friendly operators, that are easier to implement in quantum computers. 4.2.1 Fermion Inspired Pools The first ADAPT-VQE article [35] proposed pools with operators of fermionic inspiration. The ansätze resulting from pools of this type are remarkably similar to the UCCSD operator, in that they also typically contain exponentiated single and double excitations. Of course, the ordering will usually be different, as it is dictated by the system in study; further, the number of operators in the ADAPT-VQE will vary (not only with respect to UCCSD-VQE, but also within the algorithm, depending on the molecule and on the convergence criterion). The operators in fermionic pools respect fermionic symmetries: namely, particle number, 𝑆𝑍and 𝑆2 preservation, and the anticommutation requirement. A simple 𝑁-order fermionic excitation consists of exciting 𝑁electrons from some 𝑁orbitals to other 𝑁 orbitals. Often, generalized excitations are used: this means that we include operators that excite fermions from occupied orbitals to occupied orbitals, for example, rather than limiting ourselves to occupied-to-virtual excitations. The fermionic pools used in this project were taken or adapted from the simulation code of the original ADAPT-VQE article [35], that can be found in the public GitHub repository [4]. 4.2.1.1 Spin Adapted Generalized Singles and Doubles One option is the SGSD pool, also called Spin-Adapted Generalized Singles and Doubles (SAGSD). This pool was compared against the qubit-inspired pool of qubit-ADAPT-VQE in the article that introduced the latter [85]. In this pool, the double excitation operators are split into their singlet and triplet components. 76 4.2. ORIGINAL POOL CHOICES In the first quantization formalism, spin-adapted single 𝜏1and double 𝜏2excitation operators can be written in the basis of eigenstates of the total 𝑆𝑧and 𝑆2operators as: ˆ𝜏1=|𝑠=1/2,𝑠𝑧=1/2i𝑝h𝑠=1/2,𝑠𝑧=1/2|𝑞 + |𝑠=1/2, 𝑠𝑧=−1/2i𝑝h𝑠=1/2, 𝑠𝑧=−1/2|𝑞 −ℎ.𝑐., ˆ𝜏2,𝐴 =−2 √12 (|𝑠=1, 𝑠𝑧=1i𝑟𝑠 h𝑠=1, 𝑠𝑧=1|𝑝𝑞 + |𝑠=1,𝑠𝑧=0i𝑟𝑠 h𝑠=1, 𝑠𝑧=0|𝑝𝑞 + |𝑠=1,𝑠𝑧=−1i𝑟𝑠 h𝑠=1, 𝑠𝑧=−1|𝑝𝑞) −ℎ.𝑐., ˆ𝜏2,𝐵 =−1 2|𝑠=0,𝑠𝑧=0i𝑟𝑠 h𝑠=0, 𝑠𝑧=0|𝑝𝑞 −ℎ.𝑐., (4.2) where 𝜏2,𝐴 is the triplet operator and 𝜏2,𝐵 is the singlet operator, and ℎ.𝑐. denotes the Hermitian conjugate of the expression before. Of course, not all sets of four spin-orbitals allow all these excitations. If 𝑝,𝑞or 𝑟,𝑠are spin-orbitals with the same spatial part, then there is only the singlet operator, since the triplet wave function is symmetric and thus can’t be the spin part of the wave function of two fermions with the same spatial part (as the antisymmetry principle would not be respected). Rewriting these operators in the basis of tensor products of single particle 𝑆𝑧and 𝑆2eigenstates, and using the convention that an overbar on the spatial orbital index implies a 𝛽(down) spin (while no overbar implies 𝛼, or up, spin), we get: ˆ𝜏1=|𝑝ih𝑞|+|¯ 𝑝ih¯𝑞|−ℎ.𝑐., ˆ𝜏2,𝐴 =−2 √12 |𝑟𝑠ih𝑝𝑞|+1 2(|𝑟¯𝑠i+|¯𝑟𝑠i) (h𝑝¯𝑞|+h¯ 𝑝𝑞|)+ |¯𝑟¯𝑠ih¯ 𝑝¯𝑞|−ℎ.𝑐., ˆ𝜏2,𝐵 =−1 2(|𝑟¯𝑠i−|¯𝑟𝑠i) (h𝑝¯𝑞|−h¯ 𝑝𝑞|)−ℎ.𝑐. (4.3) And finally, we can transform this into the second quantization formalism to obtain the operators in 4.4. ˆ𝜏1=𝑎† 𝑝𝑎𝑞+𝑎† ¯ 𝑝𝑎¯𝑞−ℎ.𝑐., ˆ𝜏2,𝐴 =1 √12 2𝑎† 𝑟𝑎𝑝𝑎† 𝑠𝑎𝑞+2𝑎† ¯𝑟𝑎¯ 𝑝𝑎† ¯𝑠𝑎¯𝑞+𝑎† 𝑟𝑎𝑝𝑎† ¯𝑠𝑎¯𝑞+𝑎† ¯𝑟𝑎¯ 𝑝𝑎† 𝑠𝑎𝑞+ 𝑎† 𝑟𝑎¯ 𝑝𝑎† ¯𝑠𝑎𝑞+𝑎† ¯𝑟𝑎𝑝𝑎† 𝑠𝑎¯𝑞−ℎ.𝑐., ˆ𝜏2,𝐵 =1 2𝑎† 𝑟𝑎𝑝𝑎† ¯𝑠𝑎¯𝑞+𝑎† ¯𝑟𝑎¯ 𝑝𝑎† 𝑠𝑎𝑞−𝑎† 𝑟𝑎¯ 𝑝𝑎† ¯𝑠𝑎𝑞−𝑎† ¯𝑟𝑎𝑝𝑎† 𝑠𝑎¯𝑞−ℎ.𝑐. (4.4) 77 CHAPTER 4. PROBLEM TAILORED, DYNAMICALLY CREATED ANSÄTZE: ADAPT-VQE For convenience, the fermionic anticommutation relations 2.53 were used to reorganize the order of the operators; for this purpose, they are more conveniently formulated as in 4.5. 𝑎† 𝑖𝑎† 𝑗=−𝑎† 𝑗𝑎† 𝑖, 𝑎𝑖𝑎𝑗=−𝑎𝑗𝑎𝑖, 𝑎𝑖𝑎† 𝑗=𝛿𝑖𝑗 −𝑎† 𝑗𝑎𝑖(4.5) The formulation of the pool operators using creation and annihilation operators 4.4 is particularly useful, as they are directly mapped to qubit operators via the Jordan-Wigner transformation 2.57. The OpenFermion [2] library allows creating instances of the FermionOperator class directly from expressions like this one. Under the Jordan-Wigner transformation, it would appear that each product of four fermionic ladder operators is mapped to a linear combination of 16 Pauli strings; however, the Pauli strings with real coefficients don’t change under conjugation and are thus canceled out by the Hermitian conjugate, which is why adding it forces unitarity in the first place. The Pauli strings with imaginary coefficients simply switch sign under conjugation, so that each anti-Hermitian sum of four body operators results in a linear combination of only 8 Pauli strings. As a consequence, the operators in the SGSD pool consist of a linear combination of at most 48 Pauli strings. Even though the order of the excitations is the same, because of the anticommutation string, these Pauli strings act on more qubits for larger systems: their length grows on average linearly with the size of the system. The number of operators in this pool scales as O(𝑁4), where 𝑁is the number of spin-orbitals. 4.2.1.2 Spin Complemented Generalized Singles and Doubles In the SCGSD pool, each excitation is grouped along with its spin-complement (i.e. the excitation acting on the opposite spin-orbitals of the same spatial orbitals). In the first quantization formalism, they can be written ˆ𝜏1=|𝑝ih𝑞|+|¯ 𝑝ih¯𝑞|−ℎ.𝑐., ˆ𝜏2,𝐴 =−(|𝑟𝑠ih𝑝𝑞|+|¯𝑟¯𝑠ih¯ 𝑝¯𝑞|)−ℎ.𝑐., ˆ𝜏2,𝐵 =−(|𝑟¯𝑠ih𝑝¯𝑞|+|¯𝑟𝑠ih¯ 𝑝𝑞|)−ℎ.𝑐., ˆ𝜏2,𝐶 =−(|𝑟¯𝑠ih¯ 𝑝𝑞|+|¯𝑟𝑠ih𝑝¯𝑞|)−ℎ.𝑐. (4.6) And if we transform them into the second quantization formalism, we obtain (using equations 4.5 to reorder the expressions once again) ˆ𝜏1=𝑎† 𝑝𝑎𝑞+𝑎† ¯ 𝑝𝑎¯𝑞−ℎ.𝑐., ˆ𝜏2,𝐴 =𝑎† 𝑟𝑎𝑝𝑎† 𝑠𝑎𝑞+𝑎† ¯𝑟𝑎¯ 𝑝𝑎† ¯𝑠𝑎¯𝑞−ℎ.𝑐., ˆ𝜏2,𝐵 =𝑎† 𝑟𝑎𝑝𝑎† ¯𝑠𝑎¯𝑞+𝑎† ¯𝑟𝑎¯ 𝑝𝑎† 𝑠𝑎𝑞−ℎ.𝑐., ˆ𝜏2,𝐶 =𝑎† 𝑟𝑎¯ 𝑝𝑎† ¯𝑠𝑎𝑞+𝑎† ¯𝑟𝑎𝑝𝑎† 𝑠𝑎¯𝑞−ℎ.𝑐. (4.7) 78 4.2. ORIGINAL POOL CHOICES The operators in this pool consist of a linear combination of at most 16 Pauli strings, whose length also scales linearly (on average) with 𝑁, the number of spin-orbitals; and once again, the number of operators in the pool scales like O(𝑁4). While the SCGSD pool has more operators consisting of less Pauli strings on average than the SGSD pool, the asymptotic scaling with the system size remains unchanged between the two. 4.2.2 Qubit Inspired Pools The first pools of qubit inspiration to be introduced broke the fermionic excitation operators down to their constituent Pauli strings. From that, other qubit pools were proposed with the goal of further decreasing the circuit depth per operator and the size of the pool. In general, operators in qubit pools don’t respect fermionic symmetries. They usually abdicate of particle number, 𝑆𝑍and 𝑆2preservation, and may also forgo anticommutation. The pools to follow were presented in the article introducing Qubit-ADAPT-VQE [85]. 4.2.2.1 Pauli String Pool The first pool to be considered in Qubit-ADAPT-VQE, then called simply Qubit Pool, was the one consisting of all individual Pauli strings appearing in the SGSD pool (or equivalently, the SCGSD one: they result in the same set of Pauli strings, even if doubles might differ). As a simple example, we can consider the following (fermionic) single excitation. 0.5𝑎† 0𝑎4+0.5𝑎† 1𝑎5−0.5𝑎† 4𝑎0−0.5𝑎† 5𝑎1 This excitation exists (when there are more than 6 spin-orbitals) both in the SCGSD and SGSD pool, since they don’t differ in single excitations. The two last terms are the Hermitian conjugates of the former two, and OpenFermion’s [2] orbital ordering was used, so that 0, 1 and 4, 5 are pairs of spin-complements corresponding to the same spatial orbital. Under the Jordan-Wigner transformation, this becomes 0.25𝑖·𝑋0⊗𝑍1⊗𝑍2⊗𝑍3⊗𝑌4+0.25𝑖·𝑋1⊗𝑍2⊗𝑍3⊗𝑍4⊗𝑌5 −0.25𝑖·𝑌0⊗𝑍1⊗𝑍2⊗𝑍3⊗𝑋4−0.25𝑖·𝑌1⊗𝑍2⊗𝑍3⊗𝑍4⊗𝑋5. This excitation then gives rise to four operators in the Qubit Pool: 𝑖·𝑋0⊗𝑍1⊗𝑍2⊗𝑍3⊗𝑌4, 𝑖 ·𝑋1⊗𝑍2⊗𝑍3⊗𝑍4⊗𝑌5, 𝑖·𝑌0⊗𝑍1⊗𝑍2⊗𝑍3⊗𝑋4, 𝑖 ·𝑌1⊗𝑍2⊗𝑍3⊗𝑍4⊗𝑋5. The norm of the coefficients no longer matters because now each Pauli string will have its own variational parameter. 79 CHAPTER 4. PROBLEM TAILORED, DYNAMICALLY CREATED ANSÄTZE: ADAPT-VQE The number of Pauli strings per fermionic excitation remains approximately the same on average as the size of the system grows. Thus, the size of the pool is still O(𝑁4),𝑁the number of spin-orbitals - albeit with a significantly larger prefactor that implies a larger number of measurements per iteration. The Qubit Pool operators also act on a number of qubits that grows like O(𝑁). However, the circuit depth associated with each operator is significantly changed (in this case, improved). While SGSD operators contain up to 48 Pauli strings, each operator in the Qubit Pool consists of a single one. This implies an up to almost 50-fold reduction in the circuit depth associated with a single operator, and is the motive behind the Qubit-ADAPT-VQE proposal. 4.2.2.2 Pauli String Pool Without the Jordan-Wigner String The Qubit-ADAPT-VQE article also introduced a modified version of the Qubit Pool, consisting of the same operators, now cleared of the Jordan-Wigner string. This is the Î𝑖−1 𝑘=1𝑍𝑘factor that appears in the JordanWigner mapping (formula 2.57) to account for the anticommutation of fermions. The four operators introduced by the same single fermionic excitation used as an example before would simply be 𝑖·𝑋0⊗𝑌4, 𝑖 ·𝑋1⊗𝑌5, 𝑖 ·𝑌0⊗𝑋4, 𝑖 ·𝑌1⊗𝑋5. The scaling of the total number of operators in thispool is once again O(𝑁4), with 𝑁being the number of spin-orbitals. However, the scaling of the length of the Pauli strings has changed from O(𝑁)to O(1). Without the Jordan-Wigner string, only qubits representing spin-orbitals directly involved in the original excitations will be acted on by the operators they give rise to. Since the excitations we are considering are at most doubles, the Pauli strings will have length four at most. 4.2.2.3 Minimal Pools Finally, the article that introduced Qubit-ADAPT-VQE proposed a way of creating minimal pools containing operators capable of transforming any real state into any real state. These minimal pools were proven to be complete and allow convergence when used on random real Hamiltonians. An example of one such pool is the pool G, presented in the original article and reproduced here in 4.8. 𝐺1=𝑖𝑍 ⊗𝑌⊗𝐼⊗(𝑁−2), 𝐺2=𝑖𝐼 ⊗1⊗𝑍⊗𝑌⊗𝐼⊗(𝑁−3), ..., 𝐺𝑁−2=𝑖𝐼⊗(𝑁−3)⊗𝑍⊗𝑌⊗𝐼⊗1, 𝐺𝑁−1=𝑖𝐼⊗(𝑁−2)⊗𝑍⊗𝑌, 𝐺𝑁=𝑖𝑌 ⊗𝐼⊗(𝑁−1), 𝐺𝑁+1=𝑖𝐼⊗1⊗𝑌⊗𝐼⊗(𝑁−2), ..., 𝐺2𝑁−3=𝑖𝐼⊗𝑁−3⊗𝑌⊗𝐼⊗2, 𝐺2𝑁−2=𝑖𝐼⊗𝑁−2⊗𝑌⊗𝐼⊗1 (4.8) The minimal pools proposed scale like O(𝑁)on the number of qubits, with a total of only 2𝑁−2 operators. This is a remarkable reduction from the quartic scaling from the previous pools: the number of gradients that must be evaluated per iteration is reduced from O(𝑁4)to O(𝑁). 80 4.2. ORIGINAL POOL CHOICES The minimal pools allow for a significant decrease in the circuit depth per operator, and can have a very local structure (as does 𝐺, that only ever acts on two followed qubits). Unfortunately, while they work for generic Hamiltonians (such as random real Hamiltonians), they are not adequate for the specific structure of molecular Hamiltonians that is imposed by fermionic symmetries. In particular, we can see that none of the operators in the G pool conserves the spin or particle number in the Hartree-Fock state, which has severe consequences. The Hartree-Fock state can be written using OpenFermion’s [2] orbital convention as |000...111i. Acting on this state with a parameterized 𝑌rotation acting on qubit 𝑘that is in a computational basis state |𝑥i, we get: |𝜓i=𝑒𝑗𝑌𝑘𝜃|000...i⊗|𝑥i𝑘⊗|...111i =cos 𝜃|000...i⊗|𝑥i𝑘⊗|...111i+ (−1)1−𝑥sin 𝜃|000...i⊗|1−𝑥i𝑘⊗|...111i =cos 𝜃|𝑛i+ (−1)1−𝑥sin 𝜃|𝑚i Clearly, |𝑛iand |𝑚irepresent computational basis states with a different number of electrons (different number of 1 states). Considering that the Hamiltonian matrix is real and Hermitian, so that 𝐻𝑚,𝑛 =𝐻𝑛,𝑚, the energy and gradient evaluated in this state are 𝐸|𝜓i=h𝜓|𝐻|𝜓i=cos2𝜃𝐻𝑛,𝑛 +sin2𝜃𝐻𝑚,𝑚 +2(−1)1−𝑥sin 𝜃cos 𝜃𝐻𝑚,𝑛, 𝑑𝐸|𝜓i 𝑑𝜃 𝜃=0 =2 cos 𝜃sin 𝜃(𝐻𝑚,𝑚 −𝐻𝑛,𝑛) +2(−1)1−𝑥cos (2𝜃)𝐻𝑚,𝑛𝜃=0 =2(−1)1−𝑥𝐻𝑚,𝑛 Since |𝑛i,|𝑚icorrespond to Slater determinants different numbers of electrons, 𝐻𝑚,𝑛 is zero and the gradient vanishes. Since |𝑛irepresents the Hartree-Fock state, the lowest-energy Slater determinant, its energy corresponds to the lowest value in the diagonal of the Hamiltonian. Thus, it is evident that a single parameterized 𝑌rotation (acting on the Hartree-Fock state) has null gradient and does not lower the energy beyond the Hartree-Fock energy. All operators in pool G are either 𝑌rotations or conditional 𝑌rotations, so that the the same applies. Since all gradients will be zero, the ADAPT-VQE algorithm will be prevented from starting in the first place, because the selection method will fail. This is a consequence of the minimal pool operators not respecting particle number and spin symmetries. For the ADAPT-VQE algorithm to be capable of starting (by adding the first operator to the ansatz), the pool must include operators that preserve symmetries in the Hartree-Fock state. However, this still does not guarantee proper convergence; for the ground state to be reached, the pool operators must obey extra restrictions. An explicit approach for creating symmetry-adapted minimal pools was introduced in 2021 in [75]. Because of the added restrictions, the symmetry-adapted version of the minimal pools can 81 CHAPTER 4. PROBLEM TAILORED, DYNAMICALLY CREATED ANSÄTZE: ADAPT-VQE actually have less operators than the minimal complete pools presented before. The size of the pool is still O(𝑁),𝑁the number of qubits/spin-orbitals. This means that the measurements per iteration scale as O(𝑁5)(against O(𝑁8)in fermionic-ADAPT-VQE or qubit-ADAPT-VQE with the original pools), a more modest increase against the O(𝑁4)scaling of static ansatz VQE. 4.3 Application to LiH In order to analyse the femionicand qubit-ADAPT-VQE performance, and compare ADAPT-VQE with UCCSD-VQE, the algorithm was applied to multiple molecules. 𝐿𝑖𝐻 is a molecule represented by 12 qubits in a minimal basis. Its Hartree-Fock ground energy is significantly distant from the FCI value, and the Hartree-Fock approximation can be significantly improved upon in this case. Unlike 𝐻2, a smaller molecule, more than a single operator is required to reach chemical accuracy in the case of 𝐿𝑖𝐻, resulting in a more interesting evolution of the ansatz. While the simulations are heavier computationally, ADAPT-VQE with a threshold 𝜖of 0.01 is still viable, and this suffices to reach an interesting accuracy. Because of this, 𝐿𝑖𝐻 was chosen to test the algorithm in a fully noise-free scenario. 4.3.1 Evolution of the Optimization Figure 23 shows the evolution of the ADAPT-VQE algorithm for 𝐿𝑖𝐻 at an interatomic distance of 1.45Å. It can be seen that the energy reaches chemical accuracy after four iterations; at this point, the ansatz has 4 spin-adapted operators and the corresponding 4 variational parameters. For reference, the UCCSD ansatz obtained from OpenFermion [2] has for the same molecule 44 variational parameters and 64 (simple) excitation operators. Evidently, this also corresponds to a greater accuracy than the 4-operator ansatz of ADAPT-VQE; but in the plot we see that as we let the algorithm evolve further, it manages to steadily improve the energy estimate. In the original article [35], through simulations including several molecules, it was shown that ADAPT-VQE is capable of outperforming UCCSD-VQE. 82 4.4. APPLICATION TO MOLECULAR HYDROGEN Figure 23: Evolution of the fermionic-ADAPT-VQE energy and gradient norm along the iterations, for the lithium hydride (𝐿𝑖𝐻) molecule at an interatomic distance of 1.45Å. The convergence criterion was set to a gradient treshold of 0.01. The SGSD pool and the COBYLA optimizer were used. The shaded blue area represents the region of chemical accuracy (error of less than 1kcal/mol). The simulation was done via matrix algebra, and does not include any type of noise. 4.4 Application to Molecular Hydrogen Like UCCSD-VQE, the Qubit-ADAPT-VQE (with the Pauli string pool cleared of the Z string) was applied to the hydrogen molecule. The small size of this molecule allows it to be represented by only 4 qubits, and results in shallower ansätze. Because of this, it was chosen to run the algorithm on real quantum computers and analyse the effect of noise. Since it is enough to reach chemical accuracy in this case, a single iteration of ADAPT-VQE was used to obtain the ground state. At this point, there is a single variational parameter in the ansatz, which can be implemented with a CNOT count and depth of 6 assuming all-to-all connectivity. In this section, ADAPT-VQE results will be presented and compared against UCCSD-VQE, with a focus on noise-resilience and (to a lesser extent) measurement costs. 83 CHAPTER 4. PROBLEM TAILORED, DYNAMICALLY CREATED ANSÄTZE: ADAPT-VQE Neither of the algorithms manages to reach chemical accuracy with this backend and shot count, at least in this particular iteration. Regardless, ADAPT-VQE secures the energy close to the Hartree-Fock energy, with a final error of the order of centesimals of Hartree. In contrast, the UCCSD-VQE outputs an energy with an error greater than unity. The appearance of the UCCSD-VQE optimization plot is a symptom of a noise-induced barren plateau : due to the excessive circuit depth, noise washes away the quantum information and erases the characteristics of the cost function, causing the optimizer to vainly search through a flat optimization landscape. When transpiled onto the Belem backend, the single-iteration ADAPT-VQE ansatz has 12 CNOTs, against 130 from UCCSD-VQE. This evidently impacts the purity of the state by the end of the circuit. As was done in the previous chapter, an approximate density matrix representing the state at the end of the noisy ADAPT-VQE ansatz was obtained by simulating the algorithm on the QASM simulator with a noise model aiming to mimic the behaviour of the Belem processor. The fidelity of the final state with the ideal state at the end of the Adapt ansatz was found to be 0.92 (against only 0.081 in UCCSD-VQE), while the fidelity with the fully mixed state was 0.087 (against 0.40). The final purity of the state was found to be around 0.84 (against 0.16). Figure 29 compares the noise resilience of ADAPT-VQE and UCCSD-VQE against different sources (thermal relaxation, SPAM errors, and sampling). The same noise sources were already explored individually in 21 and 27; however, now the scales were adjusted so that both curves start when ADAPT-VQE is still outside of chemical accuracy and finish when UCCSD-VQE is already inside. To reach chemical accuracy in more than 50% of the runs, UCCSD-VQE required roughly 850 times greater coherence times than ADAPT-VQE. The difference in performance is visible in the plot 29a, with ADAPT-VQE showing around one order of magnitude greater accuracy for any given value of T1 and T2. This is a consequence of the greater circuit depth required by the UCCSD ansatz. Figure 29c shows the impact of sampling noise, and figure 29b the impact of SPAM noise. In order to reach chemical accuracy in at least half of the runs in the presence of sampling noise exclusively, ADAPTVQE required only 5% of the number of shots required by UCCSD-VQE.ADAPT-VQE also tolerated roughly 150 times higher error probability in state preparation and measurements. It is interesting to see how both of these noise sources have greater impact in the performance of UCCSD-VQE, regardless of being independent of circuit depth. Once again it becomes clear that more variational parameters and more variational flexibility also imply an added difficulty in the optimization, especially in the presence of noise. The noise resilience of ADAPT-VQE against UCCSD-VQE has become evident, but it must be noted that this is not the only advantage of ADAPT-VQE against UCCSD-VQE: in the original article [35], it was shown that while the latter often fails to reach chemical accuracy for strongly correlated molecules, ADAPT-VQE manages to reach it as long as the convergence criterion is sufficiently ambitious. These simulations were not replicated due to the heavy computational demands involved. 90 4.4. APPLICATION TO MOLECULAR HYDROGEN (a) Thermal relaxation (b) SPAM errors (c) Sampling noise Figure 29: Comparison of the effect of different types of noise on ADAPT-VQE and UCCSD-VQE for 𝐻2at an interatomic distance of 0.74Å. The same noise models described in figures 21 and 27 were used. 91 CHAPTER 4. PROBLEM TAILORED, DYNAMICALLY CREATED ANSÄTZE: ADAPT-VQE While the results regarding accuracy and noise-resilience favour the latter, comparing the performance of UCCSD-VQE and ADAPT-VQE is a delicate matter that goes beyond the impact of noise or the final error in the energy. There are many different costs at play, and many links between them. ADAPT-VQE requires evaluating the expectation value of a significantly larger number of observables, because all gradients must be evaluated each iteration. Additionally, because each iteration requires re-optimizing all the parameters, there will be a significant accumulated number of optimizations by the end of the algorithm, unlike in UCCSD-VQE (which requires a single optimization). This makes it seem likely that ADAPT-VQE will imply a larger number of measurements in total. However, attention must be paid to the fact that since ADAPT-VQE creates the ansatz from scratch, most of the optimizations will have a significantly lower number of variational parameters. Because they’re lower dimensional, these optimizations are likely to be easier, and to require a lower number of function evaluations (i.e. calls to the quantum computer). Further, in the common case that the final ADAPT-VQE ansatz includes less variational parameters than the UCCSD ansatz, ADAPT-VQE will tolerate sampling noise better, and thus require a lower number of shots per term for the same accuracy. Additionally, the measurement overhead from the gradient measurements may be softened by the existence of common factors among the Pauli strings that occur in the [ˆ 𝐻, ˆ 𝐴𝑖]observables (and in the Hamiltonian itself). This allows reusing measurements. Figure 30 aims to carefully analyse all of the different factors that weigh in on the costs. The number of shots for each algorithm was chosen so that the precision was matched; this implied a significantly larger number of shots for UCCSD-VQE (32 times larger, in this case). As explained before, this is due to the fact that the optimization is higher-dimensional, thus more difficult, and more sensitive to noise. Figure 30b compares the number of energy evaluations throughout the optimizations. This number is larger for UCCSD-VQE; once more, we can attribute this to the larger number of variational parameters in the optimization. However, this plot does not consider the fact that the number of observables measured in ADAPT-VQE is larger, because the gradients of the pool operators will also be measured in each iteration. The additional measurements required for the evaluation of the gradients in ADAPT-VQE are factored in in figure 30c. Here, the number of observables includes all observables measured along a complete run of either algorithm. For ADAPT-VQE, this includes all energy and gradient measurements; for UCCSD-VQE, it matches the number of energy evaluations. In both cases, the utilized optimizer was gradient-free, so that all calls to the quantum computer throughout the optimization consisted of energy evaluations. As expected, the total number of measured observables was larger for ADAPT-VQE. In this case, the pool had 20 operators, implying an additional 20 observables to be measured each iteration (the [ˆ 𝐻, ˆ 𝐴𝑖]). Finally, figure 30d compares a very important metric: the total number of shots in a full run. This number includes all the shots used for the evaluation of all Pauli strings in all observables in all times they were measured. This plot thus reflects a multitude of factors: it weighs in the fact that UCCSD-VQE requires more shots per string and more energy evaluations per optimization, as well as the fact that ADAPT-VQE requires additional measurements per iteration to obtain the gradients. The total number of shots is the total number of circuits, repeated or not, that have to be executed on the quantum computer. 92 4.4. APPLICATION TO MOLECULAR HYDROGEN (a) Error in the energy (b) Number of energy evaluations (c) Number of measured observables (d) Total number of shots Figure 30: Plots comparing the costs associated with UCCSD-VQE and ADAPT-VQE for a similar final error, in a setting with no noise other than sampling noise. For each interatomic distance, 10 runs of the algorithms were performed. The median of the errors is plotted in figure 30a. Figures 30b,30c, and 30d present respectively the average number of energy evaluations, measured observables, and total shots per run. The number of shots per Pauli string was set to 213 and 218 for ADAPT-VQE and UCCSD-VQE respectively. These values were chosen so that the error was roughly matched between the two algorithms (the average error was 0.00038 a.u. and 0.00039 a.u. by the same order). 93 CHAPTER 4. PROBLEM TAILORED, DYNAMICALLY CREATED ANSÄTZE: ADAPT-VQE Interestingly, ADAPT-VQE required a smaller number of shots in total than UCCSD-VQE for the same final error. In addition to implying shallower circuits and a lower number of variational parameters, ADAPT-VQE does not in this case imply an overhead in measurements as compared to UCCSD-VQE - on the contrary. It required only a fraction of the number of shots in total (1.7%). In this example, the fact that sampling noise is better tolerated by ADAPT-VQE allows for a reduction of the number of shots per Pauli string that compensates for the fact that there will be a larger number of observables being measured. Additionally, it must be mentioned that as was conjectured before, there was a significant amount of repeated Pauli strings among these observables. For the molecule in question (𝐻2), the total number of Pauli strings in the [ˆ 𝐻, ˆ 𝐴𝑖]observables was 152. However, the set of Pauli strings that were unique among these and that didn’t occur in the Hamiltonian, i.e. the Pauli strings that should actually be measured, was found to consist of only 24 elements. Evidently, acknowledging this fact allows for a reduction of the total number of measured Pauli strings, which in turn reduces the total number of shots. This was already considered in obtaining the results plotted in figure 30d. Of course, these results should be interpreted with caution: we’re only considering a particular molecule, and a very small one. As the size of the system grows, so will the number of operators in the pool, increasing the measurement overhead of ADAPT-VQE. The total number of optimizations required for the same accuracy will grow as well. On the other hand, the difficulty that is brought on by the extra variational parameters in UCCSD-VQE will also increase with the size of the system. This algorithm will then become even more sensitive to sampling noise, and thus require an even larger precision on the energy evaluations (which will already require a larger number of shots for the same precision given that the system is bigger). Additionally, it has been established that ADAPT-VQE is better suited for real quantum processors, as the resulting circuits are shallower. The most evident benefit of this is that it improves the viability of the algorithm in near-term devices. But given that additional sources of noise also imply an additional difficulty in the optimizations, their presence is also likely to benefit the relative performance of ADAPT-VQE against UCCSD-VQE in what comes to measurement costs. With more qubits and deeper circuits, errors will accumulate faster, so this is another factor that is affected by the size of the system. Taking all of this into consideration, it is not certain how the two algorithms will compare in terms of the total number of required shots for larger molecules and in more realistic settings. Such simulations were not attempted on account of the prohibitive computational complexity. 4.5 Application to H4 The 𝐻4molecule is represented by 8 qubits in a minimal basis. While it is already too large a system to allow for meaningful results to be recovered in real quantum computers, it is still is significantly easier to simulate classically than 𝐿𝑖𝐻. Because it allows for more iterations to be simulated in useful time, 𝐻4 was used to analyse the effect of removing the Jordan-Wigner string from the Qubit Pool operators, and to compare fermionic-ADAPT-VQE against Qubit-ADAPT-VQE. 94 4.5. APPLICATION TO H4 4.5.1 Qubit Pool: Effect of the Jordan Wigner String Figure 31: Impact of the Jordan-Wigner string on the convergence of the qubit-ADAPT-VQE algorithm. The graphs compare the evolution of the error in the energy along the 20 first iterations of the algorithm, with and without the Z-strings in the pool operators. The molecule in study is the 𝐻4molecule at an interatomic distance of 1.45Å. In one of the Qubit Pools introduced before, the chain of Pauli Z operators responsible for faithful translation of the anticommutation of fermions into the qubit states was eliminated from these Pauli strings. The effect of its removal on the convergence of the algorithm for the 𝐻4molecule is plotted on figure 31. The plot shows that there is very little variation between qubit pools with and without the JordanWigner string. This is a remarkable fact. Firstly, it is remarkable because it means that representing the fermionic anticommutation accurately through the ansatz operators brings no benefit. But most of all, it is remarkable because it greatly aids in reducing circuit depth, as was explained before. These Z strings are non-local and act on a number of qubits that scales on average linearly with the size of the system. After relinquishing them, only the qubits corresponding to the one-body wave functions directly involved in the excitations will be acted on by the operators. As a consequence, if the employed excitations are at most 95 CHAPTER 4. PROBLEM TAILORED, DYNAMICALLY CREATED ANSÄTZE: ADAPT-VQE doubles, the pool operators will act on at most four qubits regardless of the molecule in study. We saw that the circuit depth required per operator will be O(1)instead of O(𝑁). 4.5.2 Fermionic vs Qubit Pools: Performance Comparison The purpose of the qubit pool is to make the ADAPT-VQE algorithm more NISQ-Friendly, in hopes of enabling an earlier practical implementation. As it was explained in the previous section, the operators in the qubit pool are in fact implemented by shallower circuits than those in the fermionic pool. However, because they don’t preserve fermionic symmetries and consist each of a single Pauli string, they are bound to make convergence slower as a function of the number of iterations. An 𝑁-operator ansatz from fermionic-ADAPT-VQE will typically imply a much greater accuracy than an 𝑁-operator ansatz from qubitADAPT-VQE. However, it will also imply a much longer circuit. In a fair comparison, circuit depth, the number of variational parameters, and the measurement overhead should all be taken into consideration. (a) Error as a function of the number of iterations of the algorithm. (b) Error as a function of the number of CNOT gates in the ansatz circuit. Figure 32: Evolution of the error in the ADAPT-VQE energy along the evolution of the algorithm for the SCGSD pool, the SGSD pool, and the qubit pool. The molecule in study is the 𝐻4molecule. Figure 32 shows the evolution of the error in the ADAPT-VQE energy for three types of pools, two fermionic ones (SCGSD and SGSD) and the qubit pool. The energy as a function of the iteration number, which corresponds to the number of operators in the ansatz, is plotted in figure 32a. As expected, we observe that qubit-ADAPT-VQE requires more iterations/operators to achieve the same accuracy as fermionic-ADAPT-VQE. In figure 32b, we can see that the tendency is reversed when we consider the error as a function of the number of CNOTs in the circuit, to which the circuit depth will be proportional. Qubit-ADAPT-VQE converges much faster than any of the fermionic pools in this case. Further, the molecule in study (𝐻4) 96 4.5. APPLICATION TO H4 is rather small: its representation in a minimal basis set requires only eight spin-orbitals (eight qubits). It was mentioned before that there is a difference in the scaling of the circuit depth necessary to implement qubit pool operators (constant) versus fermionic pool operators (linear). As such, the difference is bound to become even more notable with the increase of the system size (number of molecular orbitals). However, the error as a function of the number of CNOTs presented in 32b (and, in general, the circuit depth) does not represent the full picture, because measurement and variational parameter overhead are not considered. Because the qubit pool consists of fermionic operators broken into their constituent Pauli strings (and cleared of the anticommutation operator), it is significantly larger. For example, for 𝐻4, we already have 328 operators in the qubit pool, while SCGSD and SGSD (fermionic pools) have respectively 111 and 66 operators. Since in each iteration the gradients of all the pool operators must be measured, this represents a sizeable measurement overhead. Further, the scaling of the number of variational parameters is as in 32a, which means that qubit-ADAPT-VQE will, for a given accuracy, imply a higher-dimensional (harder) optimization than fermionic-ADAPT-VQE. This implies more effort of the classical optimizer, which in turn will demand more calls to the quantum computer, resulting in even more measurements being required. In essence, Qubit-ADAPT-VQE reduces the circuit depth necessary to reach a certain accuracy, at the expense of a greater burden on the classical optimizer and extra measurements on the quantum computer. 97 5 ADAPT-VQE: Exploring other Pool Options In this chapter, ADAPT-VQE will be tested with several pools. The pools will be introduced here and differ both from the ones used in fermionic-ADAPT-VQE [35] and from the ones used in qubit-ADAPT-VQE [85]. The purpose is to explore the importance of the similitude between pool operators and fermionic excitations, and investigate the origin of such importance, analysing whether it stems from preserving quantities such as particle number and spin, from respecting antisymmetry, or from other factors. The constitution of the ADAPT-VQE state in terms of Slater determinants, as well as the performance of the algorithm, will be analysed and compared between the different pools. 5.1 Motivation The operators in the pools from fermionic-ADAPT-VQE are fermionic excitations that, when transformed into qubit operators, can consist of a superposition of up to 48 Pauli strings, the average length of each scaling linearly with the size of the system. These operators preserve particle number and spin, respect fermionic anticommutation, and properly represent fermionic excitations. In qubit-ADAPT-VQE, the operators are decomposed into the individual Pauli strings, and the JordanWigner string is removed, causing their length to be independent of the system size. These operators do not preserve particle number or spin, do not respect fermionic anticommutation, and do not directly represent any proper fermionic excitation. In the previous chapter, we saw that while fermionic-ADAPT-VQE converges faster as a function of the iteration number, qubit-ADAPT-VQE converges faster as a function of circuit depth. The qubit pools have more operators than the fermionic pools, which implies a higher overhead from measuring the gradients in each iteration; they will also typically require a larger number of optimizations and variational parameters for the same accuracy, demanding more effort from the classical optimizer and further increasing the number of calls to the quantum computer. In compensation for this, qubit pools reduce the circuit complexity of each operator, and result in shallower circuits for a given accuracy. In the previous chapter, no options in-between the fermionic-inspired and qubit-inspired pools were tested. As well as spin-complemented and spin-adapted excitations, simple excitations can be used. These excitations can themselves be cleared of the Jordan-Wigner string while still preserving particle number and spin, and representing proper fermionic excitations up to the antisymmetry principle. And further, it is even possible to create operators that preserve (partially or in full) particle number and spin, but don’t correspond directly to a specific fermionic excitation. 98 5.2. ANALYSIS OF THE QUBIT POOL OPERATORS In order to assess the relative importance of each of these factors, several other pools are introduced and tested in this chapter. These pools find an intermediary position between those from fermionic-ADAPTVQE and those from qubit-ADAPT-VQE. Because they have more operators than the former and less than the latter, their overhead due to gradient measurements is in-between these two algorithms; and because the number of Pauli strings in each operator is lower than in the former and higher than in the latter, the depth of their circuit implementation finds also an intermediary position. We will begin by analysing the operators in the pool from qubit-ADAPT-VQE, as well as the impact that they have on the trial state. This is relevant for the definitions of the pools that ensue. After the pools have been introduced, the evolution of the state when each is used will be analysed. Finally, the performance associated with these pools will be compared. 5.2 Analysis of the Qubit Pool Operators From the Jordan-Wigner transform in definition 2.57 and the fact that excitations are subtracted of their Hermitian conjugate to assure unitarity of the exponential operators, it is evident that any indices q, p will be represented by at most two different Pauli strings in the Qubit Pool (without 𝑍strings): 𝑖·𝑌𝑞𝑋𝑝, 𝑖 ·𝑋𝑞𝑌𝑝.(5.1) And double excitations, 𝑞,𝑝,𝑠,𝑟will be represented by at most eight different Pauli strings: 𝑖·𝑌𝑞𝑋𝑝𝑋𝑠𝑋𝑟, 𝑖 ·𝑋𝑞𝑌𝑝𝑋𝑠𝑋𝑟, 𝑖 ·𝑋𝑞𝑋𝑝𝑌𝑠𝑋𝑟, 𝑖 ·𝑋𝑞𝑋𝑝𝑋𝑠𝑌𝑟, 𝑖·𝑋𝑞𝑌𝑝𝑌𝑠𝑌𝑟, 𝑖 ·𝑌𝑞𝑋𝑝𝑌𝑠𝑌𝑟, 𝑖 ·𝑌𝑞𝑌𝑝𝑋𝑠𝑌𝑟, 𝑖 ·𝑌𝑞𝑌𝑝𝑌𝑠𝑋𝑟. (5.2) Apart from the phase factor i , the coefficients do not matter because the final coefficients will be the variational parameters. These operators consist exclusively of Pauli X and Y operators, and the number of each is always odd. This is because the other terms are canceled out by the Hermitian conjugate, which is actually what leads to the unitarity of the exponentiated operators: the coefficient of these terms would be real. It was mentioned that single excitations will lead to at most two Pauli strings, and double excitations to at most eight Pauli strings in the pool. This is because not all sets of spin-orbitals lead to valid excitations. For example, single excitations will only exist between spin-orbitals of the same type (𝛼→𝛼or 𝛽→𝛽). A similar restriction exists in double excitations; if in a set there is an odd number of 𝛽orbitals and 𝛼 orbitals (e.g. 𝑝,𝑞,𝑟,𝑠are of types 𝛼,𝛼,𝛼,𝛽), no double excitations involving these four orbitals will exist. 99 CHAPTER 5. ADAPT-VQE: EXPLORING OTHER POOL OPTIONS 5.4.5.1 Considerations Because of the relations between Pauli matrices, the constituents of all of these pools are linear combinations of operators belonging to the original Qubit Pool. As compared to the latter, the number of operators in any of the One, Two, and Four Pools is decreased almost 8-fold, and almost 4-fold in the Eight Pool. The operators in these new pools consist of linear combinations of at most one, two, four and eight Pauli strings, in the same order; the increase in number implies an increase in affinity with proper fermionic excitations. It’s important to note that implementing the exponential of a linear combination of 𝑁Pauli strings does not necessarily correspond to an 𝑁-fold increase in circuit depth, as compared to the single Pauli strings in the Qubit Pool. Because they have a very specific form resembling controlled rotations, the exponentials of these newly introduced pools can be implemented more efficiently than those consisting of linear combinations of arbitrary Pauli strings. For example, reference [98] presented circuits for implementing the exponentiated Eight Pool operators ( qubit excitations ) with a CNOT count of 13 and depth of 11, against the 48 CNOT count and depth of the ladder-of-CNOTs [64] method. Figures 33 and 34 exemplify the two different circuit implementations for a specific operator (the operator in 5.8). The former illustrates the vanilla approach, while the latter shows the more efficient option. ˆ𝜏=𝑒−𝑖𝑡 8(+𝑋0𝑌1𝑋2𝑋3+𝑌0𝑋1𝑋2𝑋3+𝑌0𝑌1𝑌2𝑋3+𝑌0𝑌1𝑋2𝑌3−𝑋0𝑋1𝑌2𝑋3−𝑋0𝑋1𝑋2𝑌3−𝑌0𝑋1𝑌2𝑌3−𝑋0𝑌1𝑌2𝑌3)(5.8) 106 5.4. ALTERNATIVE POOLS Figure 33: Circuit for implementing the operator in 5.8 using the ladder-of-CNOTs method explained in subsection 2.1.3. The circuit was drawn in the IBM Quantum Composer [41]. Figure 34: Circuit for implementing the operator in 5.8 using the more efficient method introduced in [98]. The circuit was drawn in the IBM Quantum Composer [41]. 107 CHAPTER 5. ADAPT-VQE: EXPLORING OTHER POOL OPTIONS 5.5 Graphical Analysis of the State Evolution Figure 35: Evolution of the 20 first iterations of the ADAPT-VQE algorithm for three different pools: the XXYX One Pool, the XXYX Two Pool, and the original Qubit Pool (without the Jordan-Wigner string). The upper plot shows the evolution of the error in the energy. The lower plot shows the evolution of the total probability amplitude of computational basis states representing Slater determinants with an altered particle number or 𝑍spin projection. The molecule in consideration is the 𝐻4molecule at an interatomic distance of 1.45Å. Before comparing all of these pools, an illustration of the impact of ‘partial’ conservation of 𝑆𝑧and particle number on convergence is presented. Figure 35 shows a comparison of the performance of all of the pools that do not preserve in full any of these quantities. Other than the evolution of the error, the plot shows the significance of Slater determinants with altered 𝑍spin projection and particle number in the wave function, for all three pools. Even though the XXYX One Pool contains operators capable of bringing any Slater determinant into the wave function, it doesn’t seem to be capable of convergence, as the error stabilizes at around 10−3 108 5.5. GRAPHICAL ANALYSIS OF THE STATE EVOLUTION Hartree. The same does not happen for the original Qubit Pool: it becomes clear that the different Pauli string formats are important for indispensable interference effects. However, the efforts to partially preserve 𝑆𝑧and particle number in the XXYX Two Pool are clearly rewarded: this pool does not only match the performance of the Qubit Pool, it surpasses it. The upper and lower plots in figure 35 do indeed suggest a link between convergence and preservation of these quantities. It seems like wave functions that are more contaminated by Slater determinants with altered 𝑆𝑧and particle number imply a slower convergence. It is interesting to observe that the curves in the lower plot of figure 35 are non-monotonic. All curves have a minimum at iteration 13 or 14. It seems as though as while in the beginning it is possible to get closer to the ground state while decreasing the probability amplitude of Slater determinants with altered 𝑆𝑧or N, after a point the energy can’t be lowered further without increasing it. Regardless, this increase quickly comes to a halt, and does not seem to impede convergence. Figure 36 shows the evolution of the ADAPT-VQE algorithm for four different pools defined in the previous section: the XXYX One Pool, the XXYX Two Pool, the XXYX Four Pool, and the Eight Pool. The evolution of the algorithm with the original Qubit Pool (without the Jordan-Wigner string) is also plotted for reference. Other than the evolution of error, the plots show the evolution of the number of Slater determinants (computational basis states) in the superposition state, split through three groups: those with altered 𝑆𝑍, those with altered particle number, and those with unaltered 𝑆𝑍and particle number. The purpose is to contrast the convergence speed against the composition of the state in terms of Slater determinants. An interesting aspect all plots have in common is that by the end of these iterations, there are twenty fully correct Slater determinants in the state. In the case of the Eight Pool (figure 36d), this is enough to reach a very high accuracy (almost 10−7Hartree, well inside chemical accuracy, of the order of 10−3 Hartree). Further, this number stabilizes after a certain number of iterations that ranges from 11 (in the Four Pool, figure 36c) to 15 (in the Qubit Pool, figure 36e); this stabilization does not seem to imply a stabilization of the error in general. This seems to confirm that the relevant limitation after this point is not the lack of computational basis states in the superposition that correspond to Slater determinants with relevance for the ground state, rather the lack of variational freedom with respect to their coefficients. The error only does seem to stabilize after the 20 correct Slater determinants are reached when the XXYX One Pool is used (figure 36a). The reason was explained before: while this pool allows for enough degrees of freedom to bring all necessary Slater determinants into the superposition, it does not allow for enough degrees of freedom that play with interference effects. Restricting the format of Pauli strings to XXYX forbids the use of different operators in the same four spin-orbitals that would lead to different local phases (in this case, the only options are ±1, because the state is real). This should be contrasted with the results obtained using the original Qubit Pool (figure 36e). The composition of state in terms of correct and incorrect Slater determinants is remarkably similar between these two pools. However, the Qubit Pool continues to decrease the error after the number of Slater determinants has stabilized at 20. This is because the existence of different formats of Pauli strings for the same spin-orbitals (the eight in 109 CHAPTER 5. ADAPT-VQE: EXPLORING OTHER POOL OPTIONS (a) XXYX One Pool (b) XXYX Two Pool (c) XXYX Four Pool (d) Eight Pool (e) No 𝑍Qubit Pool Figure 36: Evolution of the 20 first iterations of the ADAPT-VQE algorithm for different pools. The upper plot shows the evolution of the error in the energy. The lower plot shows the number of Slater determinants in the superposition state with altered particle number, 𝑆𝑧, or neither. In 36c and 36d, altered 𝑆𝑍and particle number correspond to the same colour because the curves coincide. When the green curve is not visible, it lies under the orange one. The molecule in consideration is 𝐻4at an interatomic distance of 1.45Å. 110 5.6. GENERAL COMPARISON OF POOL PERFORMANCE 5.2) empowers ADAPT-VQE with the ability to better control interference effects, and use them to converge: namely, it can leverage them to decrease the probability amplitude of wrong Slater determinants. When we compare the XXYX One Pool against the corresponding Two Pool (figure 36b), we can find a reason for the better performance of the latter: while it does not fully preserve particle number and 𝑆𝑍, it decreases the amount of wrong Slater determinants in the state. While the One Pool finishes with ten Slater determinants with wrong 𝑆𝑧and particle number in the state, the Two Pool finishes with only two of each. The error decreases steadily up until the last iteration. Unlike in the Qubit Pool, the improved convergence is not due to different operators in the pool implying different local phases; the Two Pool operators associated with double excitations consist of a linear combination of Qubit Pool operators, but there is still only one operator for each set of four spin-orbitals (the size of the Two Pool is almost an eight of the size of the Qubit Pool). The effect of multiplying the individual Pauli strings from the One Pool by (1+𝑍𝑞𝑍𝑝𝑍𝑠𝑍𝑟) is to preserve particle number and 𝑆𝑧in most computational basis states. While this does not provide ADAPT-VQE with more control over interference effects, it reduces the presence of incorrect Slater determinants in the state. Presumably, this decreases the need for the extra options of variational freedom, leading to improved convergence. In fact, the Two Pool even outperforms the Qubit Pool. Finally, we can contrast the composition of the state associated with the Four Pool (figure 36c) against that of the Eight Pool (figure 36d). Their ease of convergence will be compared in the next section. Both the Four Pool and the Eight Pool preserve 𝑆𝑍and particle number; the difference is that the operators in the latter are in direct correspondence with proper fermionic excitations. The former is actually faster to reach the 20 correct Slater determinants in the state, even though it has less and smaller operators in the pool (with an almost 2-fold reduction in both the size of pool and the number of Pauli strings per operator). The reason why this happens is clear from the previous section: the operators in the Four Pool and in the Eight Pool have non-trivial action on four and two Slater determinants, respectively. This is also the reason why the number of Slater determinants in the Eight Pool increases at most by two upon the addition of a single operator. 5.6 General Comparison of Pool Performance Up to now, we have been analysing the relation between convergence and composition of the state in terms of Slater determinants. This section will abstract from such details and directly compare the performance of different pools, in light of the acquired knowledge. Firstly, figure 37 compares multiple One Pools and multiple Two Pools. The purpose is to verify that the evolution of the error is not overly dependent on the chosen Pauli string format, and so the results presented previously for the One, Two, and Four pools constructed from XXYX strings hold in general. In fact, while there is evidently some variation between pools (because the operators do imply different local phases), the tendency seems identical between all plotted One Pools (YXXX, XYXX, XXYX, XXXY) and between all corresponding Two Pools. 111 CHAPTER 5. ADAPT-VQE: EXPLORING OTHER POOL OPTIONS (a) One Pools (b) Two Pools Figure 37: Comparison of the convergence of the ADAPT-VQE algorithm using several different One Pools (figure 37a) and Two Pools (figure 37b) against the full Qubit Pool (without Jordan-Wigner strings). In figure 37b, a curve corresponding to a One Pool was also plotted for reference. The error in the energy is plotted as a function of the iteration number. The molecule in consideration is 𝐻4at an interatomic distance of 1.45Å. In figure 37a, we can see that all these four One Pools perform worse than the original, full Qubit Pool. In fact, all of them seem to nearly stop decreasing the error after around iteration 10, while the Qubit Pool is steadily improving the approximation to the ground state. This seems to confirm that the different formats of Pauli string are necessary for interference effects to allow convergence, when a single string is used per operator. In contrast with the performance of the One Pools, all the corresponding Two Pools (figure 37b) converge faster than the original Qubit Pool. This does indeed suggest that the preservation of 𝑆𝑍and particle number in most computational basis states accelerates convergence in general, regardless of the string format used to construct the pool. In figure 38, we can compare the performance of the ADAPT-VQE algorithm with the GSD Pool and with the Eight Pool. It should be remembered that the difference between the respective operators is that in the latter, all Jordan-Wigner strings responsible for the anticommutation of fermions are removed. As such, this figure is comparable with figure 31, that compared the performance of ADAPT-VQE using the full Qubit Pool with the Jordan-Wigner strings and without them. Once again, removing such strings barely seems to have an effect on the algorithm. As happened with the individual Pauli strings from the Qubit Pool, removing the Jordan-Wigner strings here represents a relevant saving in the circuit depth. Both using the ladder-of-CNOTs method [64] and using the more efficient method of [98], the scaling of the depth necessary to implement the Eight Pool and the GSD Pool operators is respectively O(1)and O(𝑁),𝑁the number of spin-orbitals. 112 5.6. GENERAL COMPARISON OF POOL PERFORMANCE Figure 38: Comparison of the convergence of the ADAPT-VQE algorithm using the GSD Pool and the Eight Pool. The error in the energy is plotted as a function of the iteration number. The molecule in consideration is 𝐻4at an interatomic distance of 1.45Å. Finally, figure 39 compares the performance of all four newly introduced pools. The twenty first iterations of ADAPT-VQE for 𝐻4are plotted in 39a and the ten first for 𝐿𝑖𝐻 are plotted in 39b. As expected, the observed tendency is an improvement of convergence with increased similarity with fermionic excitations. However, it is interesting to see how the XXYX Four Pool performs similarly to the Eight Pool, even though the operators consist of half the Pauli strings and the pool has nearly half the operators. In the case of 𝐿𝑖𝐻 (figure 39b), the curve corresponding to the XXYX Four Pool is not visible because it lies under the one corresponding to the Eight Pool. In the case of 𝐻4(figure 39a), the curves only overlap in some iterations. In most iterations in which this doesn’t happen, the XXYX Four Pool has a visibly better performance (iterations 7-11 and 17-19). Only in iteration 20 does the Eight Pool ansatz prepare a better approximation to the ground state. As there is no clear trend, tests on more molecules would be required to draw conclusions. It seems possible that as the energy plunges well into chemical accuracy, the Eight Pool allows fine-tuning of the wave function by allowing selection of specific double excitations, unlike the XXYX Four Pool, which only guarantees preservation of 𝑆𝑍and particle number. However, since it allows action on more computational basis states at once, there is also reason to believe that the XXYX Four Pool may hasten convergence until points of extreme accuracy. For 𝐻4, this pool outperforms the larger Eight Pool until the error is lowered 113 CHAPTER 5. ADAPT-VQE: EXPLORING OTHER POOL OPTIONS (a) 𝐻4molecule at an interatomic distance of 1.45Å. (b) 𝐿𝑖𝐻 molecule at an interatomic distance of 1.5Å. Figure 39: Impact of the choice of pool (between One, Two, Four, and Eight Pools) in the ADAPT-VQE algorithm, for two different molecules. The error in the energy is plotted against the iteration number. In figure 39b, the XXYX Four Pool and the Eight Pool correspond to the same color because the curves coincide. These curves also overlap in some iterations in figure 39a. to under 10−6, which is roughly 0.07% of what is actually necessary to be inside chemical accuracy. 114 6 ADAPT-VQE: Further Ansatz Manipulation This chapter aims to explore the possibility of introducing extra flexibility and deliberation into the evolution of the ADAPT-VQE ansatz, instead of having it defined by a single criterion of selection. Two options will be outlined, having in common the intent of keeping the preparation circuit shallower than ADAPT-VQE (while attempting to avoid a prohibitive number of additional optimizations). The first one focuses on removing operators from the ansatz along the evolution of the algorithm, with the purpose of removing the ones that performed worse than expected. The second one attempts to be more conservative in the addition of new operators, in order to avoid getting bad performing operators in the ansatz to begin with. The approaches will be compared against each other, as well as against the original ADAPT-VQE algorithm, for multiple molecules and multiple pools. This comparison will be done by using the different approaches to obtain ansätze with the same number of operators, and analysing the error resulting from each in light of the number of optimizations accumulated along the execution (which is strongly related to the measurement cost). 6.1 Removing Operators 6.1.1 Motivation The ADAPT-VQE algorithm selects the operators based on a single quantity: the magnitude of the derivative of the energy with respect to the operator coefficients 𝜃𝑖, at point zero. This is a convenient selection method, since the gradient can be measured in a quantum computer. The circuits used for these measurements at a given iteration are barely deeper than the ansatz itself at that same iteration, thus guaranteeing that the selection method itself isn’t the bottleneck of the algorithm. What is more, there being no further information, how much an operator is impacting the energy locally, at 𝜃𝑖=0, is the best possible indicator of how much it will have lowered it once it has been added to the ansatz and had its coefficient optimized. However, this is not a faultless selection method. While it seems likely that the best performing operator will be among those with the largest gradients, it is certainly possible that its gradient is not the largest one of all. Conversely, the operator with the largest gradient is not guaranteed to have a great impact on the energy. It might even be that the derivative of the selected operator was large at 𝜃𝑖=0, but upon optimizing the coefficient and exploring the optimization landscape further, the optimizer will find that the change in energy quickly comes to a halt, as a local minimum is reached, and return a state with little energy difference as compared to the previous one. Of course, there is no way to know in advance, and this criterion allows for a much more informed choice than random selection. 115 CHAPTER 6. ADAPT-VQE: FURTHER ANSATZ MANIPULATION crucial to convergence after a certain point. In dealing with terms that have been removed, there are thus two main requirements: they can’t be allowed to be added back too soon, and they can’t be forbidden from being added back forever. An approach that seems reasonable is blocking the operators that are removed until an operator is added that performs similarly. Since operators are removed when they have little impact on the energy as compared to operators added after them, it seems natural that they would be allowed back once a newcomer has an impact comparable to the one they had. There are, however, a few problems with this method, that are better illustrated by an example. The following is a depiction of the evolution of an actual simulation of ADAPT-VQE with term removal (again for the LiH molecule at an interatomic distance of 1.45Å, using the SGSD pool). Aside from showing the drawbacks of blocking operators, this will serve to demonstrate how the term removal procedure flows. The criterion for removing operators is as defined before, and a removed operator 𝑖is allowed to be selected again once an operator 𝑘added to the ansatz later produces an energy change that is no more than the double the one operator 𝑖had produced (Δ𝐸𝑖<0.5Δ𝐸𝑘). Figure 42: Relevant data structures at iteration 3 of the ADAPT-VQE algorithm, with an additional term removal feature based on blocking and unblocking operators. The molecule in consideration is 𝐿𝑖𝐻, and the pool the SGSD pool. The units of energy, omitted to avoid overcrowding of the figures, are Hartree throughout the whole chapter. Because the variational parameters are dimensionless, the gradient shares the energy units. Figure 42 represents the relevant data structures at iteration three of the algorithm. The most important one is the ansatz, that specifies the current state preparation circuit. As the name suggests, position represents the position of each operator in the ansatz: the column labeled with 0 concerns the first operator in the ansatz, the one labeled with 1 concerns the second one, and so on. The first operator is both the one that was first added and the one that comes first in the circuit. The index identifies which operator in the pool was selected; what fermionic operator this corresponds to is irrelevant, as the purpose is just attributing an identifier to each pool operator. For a matter of simplicity, this identifier is a number rather than a formula. The first three iterations are uneventful: the criterion for attempting to remove a term is never met. As can be seen in figure 42, the energy change always decreases from an iteration to the following one. As 122 6.1. REMOVING OPERATORS such, no operator is suitable to be removed, and the algorithm proceeds just as the original ADAPT-VQE. The set of blocked operators is empty. Figure 43: Relevant data structures at iteration 4. Finally, at iteration 4 (figure 43), the condition necessary to attempt to remove an operator is met. The condition reads Δ𝐸𝑖>𝑟Δ𝐸𝑗⇔ −9>0.5× (−62), which does indeed hold: the operator 191 ( j ), added in this iteration, produces a change in energy significantly higher than that of operator 25 ( i ), added to the ansatz in the previous iteration. As such, operator 25 is removed from the ansatz, and the coefficients of the remaining operators are re-optimized. From this, it is possible to calculate the increase in energy caused by removing this operator, which turns out to be 9×10−5a.u.. The condition for effectively removing the operator reads −Δ𝐸0 𝑖>𝑡Δ𝐸𝑖⇔ −9>1.5× (−9), which once again holds. It can be concluded that the impact of this operator in the energy isn’t much changed by the posterior addition of operator 191 into the ansatz. This will not necessarily always be the case: sometimes, one might find that removing an operator from the middle of the ansatz will cause a drastic increase in energy, even if adding that operator in the first place had barely lowered it. This is a consequence of the presence of the operator in the ansatz being crucial for the impact on the energy of the next ones, and possibly decisive for their selection. 123 CHAPTER 6. ADAPT-VQE: FURTHER ANSATZ MANIPULATION By the end of iteration 4 operator 25 is blocked, so as to seek operators that might have a greater impact on the energy, even if their gradient is lower. Figure 44: Relevant data structures at iteration 7. The algorithm proceeds without noteworthy developments until iteration 7. At iterations 5 and 6, operators 188 and 184 are added. They lower the energy roughly 7 and 4 times more than 25 (the blocked operator) respectively. This is a sign that the removal was worthwhile. Finally, at iteration 7, an operator is added that causes a significantly smaller change in energy. The condition for unblocking operator 25 reads Δ𝐸𝑖<0.5Δ𝐸𝑘⇔ −9<0.5× (−4), so that operator 25 is released and allowed to be selected again in future iterations. The set of blocked operators is once again found empty. Unsurprisingly, 25 is the operator selected in the next iteration (iteration 8). It is interesting to note how the gradient and energy change of this operator are now, at iteration 8 (figure 45), very similar to those back at iteration 3 (figure 42), when the operator was first added. Despite the position in the ansatz being changed (from 2 before, to 6 now), the effect of the operator on the ansatz seems to be about the same. Since the previous operator (64) had produced an energy change of −4×10−5a.u., roughly half of that seen in iteration 8, that operator is removed from the ansatz. After re-optimizing the coefficients, the energy has only increased by 4×10−5a.u., so the algorithm goes forward with the removal. At the end of iteration 8, operator 64 is blocked. By this iteration, one can already notice undesirable behaviour. Operator 25 had been removed and blocked, but hadn’t it been for that, it would precede operator 64. The removal criterion is based on subsequently added operators, but if they are only added later because they had been removed from an earlier position in the ansatz, this seems rather inadequate. 124 6.1. REMOVING OPERATORS Figure 45: Relevant data structures at iteration 8. An easy solution for this would be not allowing operators to be removed upon the addition of a previously blocked operator. However, this is not the only problem with the approach, as will become clear soon. Figure 46: Relevant data structures at iteration 9. By iteration 9 (figure 46), operator 59 is added. This operator produces an energy change similar to the one operator 64, added in iteration 7 (figure 44) and blocked in the ensuing iteration (figure 45). As such, the operator is unblocked. Nothing remarkable happens at iteration 10: operator 29 is added, and no operator is suitable to be removed upon this addition. That changes by iteration 11 (figure 47), when the previously removed operator 64 is added back to the ansatz. At this point, operator 29 is removed and blocked. Once again, it should be pointed out that this is really inefficient. Operator 64 had been removed from being outperformed by an operator that belonged to an earlier point of the ansatz, and this proved 125 CHAPTER 6. ADAPT-VQE: FURTHER ANSATZ MANIPULATION Figure 47: Relevant data structures at iteration 11. to be a bad choice: it was added shortly after, with no benefits and a cost of two extra optimization steps. Now, this very operator causes the same problem: it removes operator 29 from the ansatz, even though it would have come before it if we hadn’t removed it. Figure 48: Relevant data structures at iteration 12. In fact, this operator (29) is unblocked right in the following iteration (iteration 12, figure 48), once another operator is added (32) that causes about the same impact in the energy. And it will be added again immediately, by iteration 13. Once again, removing the operator costed us two extra optimization steps, and had no effect other than moving it a couple of positions forward in the ansatz (with no benefit in the energy). At iteration 14 (figure 49), something very interesting happens: the selected operator (35) causes an energy change of −34×10−5a.u., a change so large that it is only surpassed by that of iteration 6. Despite having a lower gradient than any other operator so far, it largely outperforms any of the five last operators, 126 6.1. REMOVING OPERATORS Figure 49: Relevant data structures at iteration 14. with an energy change of up to almost 20 times that of its predecessors. Thus, unprecedentedly, there are five operators that are liable to be removed. They are removed one by one, the coefficients being re-optimized each time; and each time, it is verified that the operator can be effectively removed. This is remarkable: removing these five operators caused an increase in energy roughly equal to the sum of the decreases that they had individually caused. This means that operator 35, added only in this iteration, was by far a better candidate to hold position 5 of the ansatz that any of those in positions 5-9 up to now. However, it was not picked until as late as iteration 14, because its gradient was not large enough. Looking at figure 49, and focusing especially on the data regarding the ansatz at the begging of the iteration, one can notice that operators in positions 6-10 all had very similar gradients upon being selected. However, there are relevant differences in the energy change that they produced: and in fact, it was the operator with the lowest gradient among these (35) that caused the greatest change in energy. Moreover, selected with a gradient of more than double and holding position 5 in the ansatz, operator 25 produced only a fraction of this energy change. This further highlights the downsides of using the gradient as the selection method, and proceeding in disregard of the energy change that the selected operator actually produced. At this point, we have an ansatz with only 6 operators that prepares a state with lower energy than any before. It is a significantly closer approximation to the ground state than the ansatz with 10 operators we had back at iteration 13, despite having a 40% decrease in the number of variational parameters and (approximately) circuit depth. But in spite of the the ansatz looking favorable at this point, the way it was reached was certainly not the best. By iteration 14, three operators (25, 29 and 64) have been removed from the ansatz twice. 127 CHAPTER 6. ADAPT-VQE: FURTHER ANSATZ MANIPULATION Removing an operator and adding it back represents an overhead in optimization costs (coefficients must be re-optimized upon removal and upon re-addition), so it is certainly beyond optimal to have them removed twice. There could be no greater sign that removed operators are being allowed back into the ansatz too soon than seeing them removed a second time. In conclusion, it seems that the strategy of blocking them until a new operator produces a change in energy similar to the one they had produced is not the most appropriate. After examining this example run, this became exceedingly evident. This is not quite so unpredictable. After we remove an operator, we know that it didn’t perform as well as could be expected given its gradient. In a first analysis, blocking it until there is a similar energy change seems adequate. However, unblocking operators based on the energy change produced by a single new operator is dangerous: there is the possibility that the new operator also impacted the energy significantly less than another one with the same gradient could. What is more, once unblocked, the operator will go back to competing with the remaining operators in the pool on the same standing. The selection method is still the gradient, so it is likely that it will be picked very soon after it’s unblocked (as did happen multiple times in the example). The information that an operator had little impact in energy is being underused: we are only taking it into account to remove and block the operator. Once it is unblocked, we allow it to be selected again by its gradient, even though we now know that it is an inadequate indicator for that operator especially. In the previous exposition, we can see that operator 25, the first being removed (iteration 4, figure 43), is allowed back into the ansatz by the addition of operator 64 (iteration 7, figure 44). However, operator 64 is later removed, signaling that it is not a good reference for unblocking others. What is more, this last operator is itself unblocked by the addition of 59 (iteration 9, figure 46), later removed (iteration 14, figure 49). By iteration 14 all the aforementioned operators are removed once again, because an operator with a smaller gradient turns out to be a far better choice than any of them. However, this operator did not have a chance to be added to the ansatz before all the removed operators were added again, because they did still have higher gradients. This is particularly bad considering how close the gradient of this operator (35) was to that of the bad performing operators. For example, at iteration 11 (figure 47), previously removed operator 64 is added back to the ansatz. In this iteration, the gradient of operator 64 is only 107% that of operator 35, the above average performing one. And even though we have no way of knowing that this operator will perform so well, we do know that operator 64 did not perform as well as it could be expected from its gradient. As such, it would definitely be interesting to give an opportunity to others with close, but slightly smaller gradients. The goal of removing terms was precisely to explore operators with lower gradients, and this goal is not being met: blocking removed operators was proved to be an inadequate option. As soon as a bad performing operator is added to the ansatz, all of the previously blocked ones are unblocked and likely to be selected. After this discussion, the solution that seems the most adequate is to somehow change the selection 128 6.1. REMOVING OPERATORS criterion so as to incorporate information regarding previously removed operators: ideally, an operator that has been removed would be less likely to be selected than others. Something as binary as blocking and unblocking operators doesn’t seem to take advantage of the information that we have to the full extent. 6.1.2.3 Performance Penalty The best way to incorporate extra factors into the selection criterion without drastically changing it would be to impose a penalty on the gradient of the removed operators. The selected operator would be the one having the highest ‘effective’ gradient, with this quantity reflecting the below average performance of those operators that had been previously removed. Instead of just having the operator blocked and later unblocked, we’d have its effective gradient lower than the actual gradient, in accordance with how poorly it performed. Evidently, it would be interesting to have some quantification of just how bad the energy change was, so as to apply a proportional penalty. The purpose would be signaling to the algorithm that, even if a removed operator has a high gradient, we know that it doesn’t change the energy as much as one could expect. As such, we want to attempt to add other operators to the ansatz, even if they have lower gradients, in hope of finding some that have a greater impact on the energy. The goal is using the knowledge that we now have (regarding the removed operators) to make a more informed selection: this penalty would be a way of correcting, in part, the flaws of the gradient as an indicator of the potential energy change. In those cases in which the gradient was a particularly poor indicator, to the point we removed the operator, the effective gradient would be lowered accordingly. The aforementioned performance ratio seems interesting for this purpose. After we select a given operator and re-optimize the ansatz, we can calculate the associated performance ratio from its gradient and the resulting energy change (as compared to the previous iteration) from formula 6.1. The energy change is what we actually want from an operator: that it brings us closer to the ground energy. As the selection method, the gradient quantifies in a way how much we expect it to change the energy: we select the operator with the highest gradient hoping that it is also the one with the greatest impact on the energy. So intuitively, this ratio can be interpreted as something like performance ratio =actual performance expected performance. While the performance ratio seems useful and has the nice property of being dimensionless, it doesn’t allow us in itself to establish a penalty: we don’t know what is a good or bad performance ratio. Once again, the easiest solution seems to be assessing performance relative to that of the remaining operators. If we take the average over all known performance ratios (meaning, over the ratios of all of the operators ever added to the ansatz), we get a very interesting quantity: the average (absolute) energy change per unit gradient. At last, we know what energy change it is reasonable to expect from an operator with a given gradient. 129 CHAPTER 6. ADAPT-VQE: FURTHER ANSATZ MANIPULATION One can then define a dimensionless penalty for removed operators as penalty =performance ratio standard performance ratio,(6.2) where the standard performance ratio is our reference - it can just be the average taken over all known performance ratios. This metric tells us just how below average the performance of the operator was. If the energy change per unit gradient of the removed operator was half of the average, the effective gradient will be half of the actual gradient. This ‘sanction’ will prevent it from being selected back into the ansatz as soon. As long as there are other operators in the pool that, considering their gradient, will have a greater impact on the energy assuming that their energy change per unit gradient is average , the removed operator will not be selected again. Only when, from the gradients of the other operators in the pool and under the same assumption, we don’t expect that any will impact the energy more, will the removed one be selected. As it was exposed before in figure 41, the performance ratio has a tendency to decrease along the iterations; that is why it wasn’t a suitable metric to decide which operators to remove. Taking the average will already soften this effect: the value will be constantly updated as more performance ratios are known, and as such will follow their tendency to decrease. To further avoid over-penalizing the removed operators, it will be interesting to use the moving average, rather than the average, as the standard performance ratio. This will prevent data from too old iterations to affect the calculation of the penalty. The window over which we calculate the average shouldn’t fall in the other extreme and be too small. If we only consider the ratio of two or three of the last added operators, and one of them performs particularly bad, it will have a significant impact on the average - thus causing the penalty to be too soft, and the removed operator to be added back too soon. It is better to be cautious, because the resulting optimization overhead can significantly worsen how fast the algorithm is to converge, by wasting iterations adding and removing the same operator multiple times. In the following examples, a window of 10 was used. The final version of the procedure of operator removal can now be stated. The basic algorithm is the same as the original ADAPT-VQE, except at every iteration 𝑗, when a new operator is added, the following must be done: • Calculate the performance ratio of the new operator (equation 6.1). Store this value (associated with the operator) and update the standard performance ratio. • Go through the previous iterations and compare the energy changes produced in them with that produced in this one. If an operator added at iteration 𝑖caused an energy change Δ𝐸𝑖such that Δ𝐸𝑖>𝑟Δ𝐸𝑗, attempt to remove this operator and re-optimize the coefficients of the ansatz without it. Calculate the energy change produced by this action, Δ𝐸0 𝑖. • Compare the increase in energy caused by removing the operator with the decrease that had been caused by adding it. If −Δ𝐸0 𝑖>𝑡Δ𝐸𝑖, remove the operator. In the following iterations, multiply 130 6.1. REMOVING OPERATORS its gradient by a penalty calculated as in equation 6.2. The penalty should be recalculated each iteration, to use the most recent value of the standard performance ratio. As it was mentioned before, one should choose 0<𝑟<1and 𝑡>1. The examples to follow used 𝑟=0.5and 𝑡=1.5. The 10-iteration moving average of the performance ratio was used as the standard. In order to compare this approach with the previous one (blocking / unblocking operators), the execution of the algorithm will be exemplified for the same test case as before, albeit more succinctly. Figure 50: Relevant data structures at iteration 3 of the ADAPT-VQE algorithm, with an additional term removal feature based on a performance penalty. The molecule in consideration is 𝐿𝑖𝐻, and the pool the SGSD pool. At iteration 3 (figure 50), the ansatz is exactly the same as previously (figure 42). Until this point, as was pointed out before, no operators have met the removal criterion. The only difference is that now we’re keeping track of the performance ratio of the operators, easily calculated from their gradient upon being selected and the energy change that they caused (equation 6.1). This ratio has decreased from iteration to iteration up to now. It was already known that the produced energy change and the gradient of the selected operator were strictly decreasing until iteration 3, but it is also interesting to note that the absolute value of the change in energy per unit gradient is also decreasing. At this point, the standard performance ratio sits at 41 ×10−3. Evidently, it has also been on a monotonous decrease. As before, the monotony is broken at iteration 4 by an operator that impacts the energy significantly more than its predecessor. We attempt to remove it and do so successfully, as the impact in the energy is still low after the addition of the last operator. But now we don’t ‘block’ operator 25: we just mark its bad performance, and keep a note that it should be penalized in the following iterations. The developments of iterations 5 to 13, with a focus on the penalized operator, are represented in figure 52. 131