Solving Differential Equations On Quantum Computers
MULTISCALE STRUCTURAL SIMULATIONS LABORATORYAEROSPACE ENGINEERING UNIVERSITY OF MICHIGAN Solving differential equations on quantum computersProf. Veera SundararaghavanDepartment of Aerospace Engineering, University of MichiganSid Srivastava (PhD candidate)Keynote Talk: Modeling and Computation session16th Pan-American Congress of Applied MechanicsMay 23, 2019Acknowledgments: USRA Quantum information Sciences Program
What is Quantum computing?Quantum computersemploy superposition ofstatesQubit stateDur and Heusler, Arxiv2013“It’s not just a question of moving more quickly. It’s a question of moving in different ways.”“It’s as though you’re Houdini trying to pick a lock and escape from an underwater cabinet. If you werefree to move your hands wherever you’d like, you could do so much more efficiently than if you werehandcuffed.”Aephraim SteinbergProfessor of physics at the University of Toronto andCentre for Quantum Information and Quantum Control
The Quantum annealerTunable field on thequbitTunable interactionbetween qubits3
Optimization using quantum annealersClassicalHikingEnergyQuantum TunnelingState4
No. of PublicationsAlgorithms were developed before the machines arrived1960-80:Quantum Information theory1980-90:Feynman’s talk – Cannot simulate Quantum system efficiently on a classical computerBenioff proposes a theoretical framework for QCDeutsch describes the first universal quantum computerEkert invents entangled based secure communication1990-98:Shor’s algorithm for factorizationGrover’s search algorithm1998:2-qubit NMR quantum computer(Present)Quantum Information TheoryShor’sAlgorithmData from Web of science
Recent developmentsPrimary focus: Identify and solve problemswhich are very hard to be solved onClassical supercomputers (e.g. NP hardcombinatorial problems).Neural Network[1]Quantum key distribution [3]Secondary focus: Solve well-establishedclassical problems on a Quantum systemmore efficiently.National Quantum Initiative Act (Jan 2019)at a funding level of 1.2 billion impactsgrowth of quantum computationalsciencesMolecular energy estimation[2]Playing Battleship [4][1] Rebentrost, Patrick, et al. "Quantum Hopfield neural network." Physical Review A 98.4 (2018): 042308.[2] Kandala, Abhinav, et al. "Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets." Nature 549.7671 (2017): 242.[3] Islam, Nurul T., et al. "Securing quantum key distribution systems using fewer states." Physical Review A 97.4 (2018): 042347.[4] https://www.research.ibm.com/ibm-q/ (Last accessed on 20 Feb 2019)
Some issues with near-term quantum computersLimited QubitsIn classical computers, with similar binary (0/1) encoding:Data type size : 32 bits (float) to 80 bits (long double)1 GB memory 12 million high precision variablesIn contrast, currently available quantum annealers have a limited number of physical qubits.Limited connectivity (Gate operations):Schematic of IBM Q5 processorUnit graph structure for D-Wave 2000Q7
Quantum annealing: Ising spin systemsHow does a Quantum annealer work?D-Wave ProcessorAluminum loopJosephson junctionNTT Basic Research Laboratories (2005)Energy in spin systems is dependent on 3 things:1. Topology of graph.2. Parameters H (Field strength) and J (Coupling strength)3. Labeling of vertices
Ising spin systemsSpin systems A spin model defines a Hamiltonian (Energy) on a simple undirected Graph for a given set of labelling. An Undirected Graph G(V,E) is a set of vertices (V) and edges (E) with no orientation. It is ‘simple’ if it does not containany multi-edge or self loop. A Vertex labeling is a function of V to a set of labels ({ 1,-1} in our case)Label 1Label -1The energy (E) for a given labeling (S) :𝐸 𝑆 𝐻𝑖 𝑆𝑖 𝑖𝐽𝑖𝑗 𝑆𝑖 𝑆𝑗⟨𝑖,𝑗⟩ 𝐻𝑖 𝐻𝑗 𝐽𝑖𝑗𝐻𝑖𝐽𝑖𝑗𝐻𝑗 𝐻𝑖 𝐻𝑗 𝐽𝑖𝑗 𝐻𝑖 𝐻𝑗 𝐽𝑖𝑗 𝐻𝑖 𝐻𝑗 𝐽𝑖𝑗
Annealing procedure The annealing procedure is conducted by varying the fieldas𝐸 𝑡 𝐴 𝑡𝑥𝑖 𝑆𝑖 𝐵(𝑡)(𝑧𝑖 𝐻𝑖 𝑆𝑖 𝑧 𝑧 𝑖,𝑗 𝐽𝑖𝑗 𝑆𝑖 𝑆𝑗) The fridge temperature used in D-Wave Vesuviusprocessor is 12mK. The total annealing time is in range of 5 μs - 2000 μsClassicalHikingEnergyWhy is this better than classicalcomputing?QuantumTunnelingState*Mishra, A., Albash, T. and Lidar, D.A., 2018. Finite temperature quantum annealing solving exponentially small gap problem with non-monotonic success probability. Naturecommunications, 9(1), p.2917.
Key topic of this talk:Mapping physics to Ising models𝐽11qi11𝑎𝑖𝑎𝑗𝐽212Node 𝑗Node 32𝐽231qj23qj3𝐽33Gate based quantum computers orQuantum annealers11
Scope of this talkPrimary objective: Formulate and test Quantum annealing based algorithms for differential equations.Adiabatic Quantum Computer(Black box)For most part, we will treat Quantum annealer as a black box which solves graph labeling in one step.Also in scope: Quantum approximate optimization on gate based quantum computersNot in Scope: Quantum linear solver-based procedureA great amount of work has been based on QLSA solver developed by Seth Lloyd [1,2]. QLSA is very promising but isnot as robust to noise which becomes important in near term quantum computers.[1] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd, Quantum algorithm for linear systems of equations, Physical Review Letters (2009), 103, no. 15, 150502, arXiv:0811.3171. (QLSA)[2] Childs, Andrew M., and Jin-Peng Liu. "Quantum spectral methods for differential equations." arXiv preprint arXiv:1901.00961 (2019).
A simple exampleCase Study I: 1-D truss problemyAxf(x)u 1B𝑑𝑑𝑢𝐸𝐴(𝑥) 𝑓 𝑥 0𝑑𝑥𝑑𝑥Dirichlet boundary conditions:L0 𝑥 𝐿𝑢 0 0𝑢 𝐿 1Goals Introducing Box algorithm for solving differential equations* Familiarizing with the D-wave quantum annealing architecture*Srivastava, Siddhartha, and Veera Sundararaghavan. "Box algorithm for the solution of differential equations on a quantum annealer." accepted forpublication in Physical Review A.
Solving differential equationCase Study I: 1-D truss problem𝑑𝑑𝑢𝐸𝐴(𝑥) 𝑓 𝑥 0𝑑𝑥𝑑𝑥yAxu 1Bf(x)Dirichlet boundary conditions:LEnergy methods:Solution obtained by minimizing the potential energy given as:𝐿min 𝜋 𝑢 01𝑑𝑢𝐸𝐴2𝑑𝑥2 𝑓 𝑢 𝑑𝑥0 𝑥 𝐿𝑢 0 0𝑢 𝐿 1
Discretization and compact basisFinite element approximation:Pick some appropriate finite dimensional space, 𝑉ℎ with basis 𝜙1 , 𝜙2 , , 𝜙𝑛 :𝑢 𝑛𝑖 1 𝑎𝑖 𝜙𝑖with 𝑎𝑖 ℝSolve min 𝜋 𝑎 to get the best approximation of 𝑢𝐿𝜋 𝑎 01𝐸𝐴22𝑛𝑎𝑖 𝜙′𝑖𝑖 1Choice of 𝜙:‘Hat functions’ (Compact support)This gives a sparse structure to the minimization problem𝑛 𝑓𝑎𝑖 𝜙𝑖 𝑑𝑥𝑖 1
Nodal graphWe want to ensure that there is only ‘ONE’ 1 and ‘TWO’ -1 on each node𝒒𝒊𝟏 𝒒𝒊𝟐 𝒒𝒊𝟑 Energy H 1, J 11113H 3J6-1 -1 -1 -3H 3J0-111H-J01-11H-J011-1H-J0-1 -11-H-J-2-11-1-H-J-21-1 -1-H-J-2𝐻𝑖 𝑎 on all nodes𝐽𝑖𝑗 𝑏 on all edgesHigher energy states3 Minimum energydegenerate states
Representation of Solution spaceThree energy minimizers (symmetric) for eachnode Three values of 𝑎𝑖 (coefficient of linearexpansion) for 𝑖 𝑡ℎ node
Element graphElement graph encodes the physics of the problem𝑎𝑖𝑎𝑗1Node 𝑗Node 𝑖Element2𝐽212 Each node can take 3 values Each element can have one of 9states (𝑎𝑖 , 𝑎𝑖 1 )1𝜋𝑒 𝐸𝐴 𝑎𝑖 1 𝑎𝑖2𝐽11qi1qi2𝐽12Estimate 𝐽 (edge strength) such that:𝐽313qi3𝐽22(𝑎𝑖 1 𝑎𝑖 ) 𝑓2qj1𝐽132𝐽32𝐽23qj13qj32𝐽33Element graph has 9 edges and 9 valid coloringsE 𝑖𝑗 𝐽𝑖𝑗 𝑆𝑖 𝑆𝑗 𝜋𝑒 (𝑎𝑖 , 𝑎𝑖 1 )System of 9 linear equations in 9 variables
Element graph : Example19
Element graph : Example (contd)20
Graph Embedding We want to map all nodes from the graph model onto the physical graph. This mapping should preserve the minimum energy states.Required ConnectivityPhysical connectivityGood News: D-wave provides you with API’s to search embedding in a heuristic fashion
Energy minimizationBoundary Conditions: 𝑎1 𝑢1 and 𝑎𝑛 𝑢3Choose 𝐻𝑖 (Field term) corresponding to 𝑞11 and 𝑞𝑛3 as large negative values.Solve for 𝑬𝑨(𝒙) 𝟏 and 𝒇(𝒙) 𝟎Low Energy SolutionHigh ProbabilityHigh Energy SolutionLow Probability100 labels per node are required to get a precision of 0.01 for a bounded displacement between [0,1].
Introduce slack variables23
Box algorithm: Iterative Procedure Define 𝑢𝑐𝑖 as the displacement corresponding to center label of 𝑖 𝑡ℎ node.𝒖𝒄 {𝑢1𝑐 , 𝑢𝑐2 , . . , 𝑢𝑐𝑛 } And a parameter, ‘𝑟’ (called slack variable) so that the displacements of 𝑖 𝑡ℎ node forcorresponding to labels {1,2,3} are {𝒖𝒊𝒄 𝒓, 𝒖𝒊𝒄 , 𝒖𝒊𝒄 𝒓}, respectively.𝑢𝑐5 𝑟𝑢𝑐2 𝑟𝑢𝑐2𝑢1𝑐 𝑟𝑢1𝑐𝑢1𝑐 𝑟𝑢𝑐2 𝑟𝑢𝑐3 𝑟𝑢𝑐3𝑢𝑐3 𝑟𝑢𝑐4 𝑟𝑢𝑐4𝑢𝑐4 𝑟𝑢𝑐5𝑢𝑐5 𝑟Definition: ‘Box’ is the highdimensional representation of allpossible outcomes. In this case5D space. Box length 2rBox center: 𝑢𝑐1 , 𝑢𝑐2 , . . , 𝑢𝑐5Approach:1. If box corner is chosen,recenter the box to the corner.2. If box center is chosen, shrinkthe box.24
Iterative solutionThe link weights are modified based on thecurrent choices of center and slack variable foreach node.25
Numerical solutionyEA 2EA 1xL 1 The yellow region representsthe space between 𝒖𝒄 𝒓 and𝒖𝒄 𝒓. The result converges to theexact solutionDisplacementAu 1Bx
Quantum computer’s solution𝐸𝐴(𝑥) 2 𝑥𝑓 𝑥 4𝑥 6
Convergence𝑎2𝑛𝑣 𝑎𝑖 𝜙𝑖𝑖 1𝑒3with 𝜙𝑖 𝑉ℎ and 𝑎𝑖 ℝBest solution in 𝑉ℎ𝑒2𝑒1Now we sample 𝑎𝑖 from {𝑢𝑐𝑖 𝑟, 𝑢𝑐𝑖 , 𝑢𝑐𝑖 𝑟}2𝑟𝑎1
Convergence𝑎2Consider the iteration, when the minimum energy pointcorresponds to 𝑢𝑐 .Observe: All other points in the sample lie outside the energy contourcorresponding to 𝐹(𝑢𝑐 ) The length of the major axis is bounded for a given matrix M2𝑟2𝑟𝑢𝑐𝑒 𝑢𝑐𝑇 𝑀𝑢𝑐𝑎1Only considering the horizontal and vertical node, you canbound the major axis of the ellipse as:𝑑𝑚𝑎𝑥 2 1 𝜆𝑚𝑎𝑥 /𝜆𝑚𝑖𝑛 𝑟Following the same logic the bound can be extended to ℝ2 :𝑑𝑚𝑎𝑥 2 1 (𝑛 1)i.e. for a finite discretization, lim 𝑑𝑚𝑎𝑥 0𝑟 0This mean as 𝑟 0 , 𝑢𝑐 approaches the best approximation for 𝑢 in 𝑉ℎ𝜆𝑚𝑎𝑥𝜆𝑚𝑖𝑛𝑟𝑛
Another exampleCase Study II: Advection-Diffusion problemHomogeneous Advection-diffusion equation with Dirichlet boundary conditions on both endsGoverning Equation: 𝑢′′ 𝑣𝑢′ 00 𝑥 10𝑢 0 0Finite Element Method (Galerkin approach)Weak-form:𝑊 𝑢, 𝑢 : ,𝑢 10 1(𝑢′ 𝑢′ 𝑣𝑢′ 𝑢) 𝑑𝑥 0Ω Observe that the weak form is non-symmetric This results in unstable solution for high values of ′𝑣′ i.e. inhighly advective flows.Goals Application of energy minimization in a restrictive way Some tweaks in the Box algorithm and speed of convergence Error correction measuresFEM (Linear elements) calculationfor 𝑣 4
Box algorithm for A-D equationFirst, we need a functional minimization formPotential flow assumption:1𝐹𝑣 [𝑢] 2𝐿𝑣 𝜙𝑒 𝑣.𝑥 𝑢′ 2 𝑑𝑥0This energy can be written for n-dimensionswith non-homogenous terms with Dirichletand flux boundary conditions*Restriction (Necessary condition for 𝑛 2) 𝑣 0Application of Box algorithm is same as truss problem*Auchmuty, Giles. "Variational principles for advection–diffusion problems." Computers & Mathematics with Applications 75.6 (2018): 1882-1886.
Slow convergence Box algorithm performs well for smaller velocities For higher velocities we use step size selection𝑟𝑛𝑒𝑤 𝛼𝑟𝑜𝑙𝑑0 𝛼 1𝑣 0.5𝑣 1𝛼 0.85𝑣 2𝑣 4
Δ𝐸Δ𝐸Summary of Case Study II𝐼𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛𝑠𝛼 0.5𝐼𝑡𝑒𝑟𝑎𝑡𝑖𝑜𝑛𝑠𝛼 0.85 We formulated and implemented Box algorithm for Advection-Diffusion (w/ Potential flow) We showed that step size selection can be used to approach the global minima
Error correctionSo far, we have treated the quantum computer as a black box which outputs the correct minima.When 𝛿 𝐸 𝐸 i.e. the relative energy of all states are same then the solver can output sub-optimal solutions.𝑬 1𝑬 0.98Incorrect minimaIdea for error correction: Scale energy to maximize thegap between the different states: (𝐻, 𝐽) (𝐻′ , 𝐽′ )𝐸′ 1𝑎(𝐸𝑖 𝑏𝑖 )𝑖 𝑒𝑙𝑒𝑚𝑒𝑛𝑡Choose a, 𝑏𝑖appropriately
BifurcationCase Study III: Beam-buckling problemyP4th Order differential equation with critical behaviorxGoverning Equation:𝐸𝐼𝑤 ′′′′ 𝑃𝑤 ′′ 00 𝑥 𝐿LBoundary conditions at x 𝑥𝑏 :𝑤 𝑥𝑏 0 (Displacement) or𝑤′ 𝑥𝑏 0 (Slope)or𝑉 𝑥𝑏 𝑤′′′ 𝑥𝑏 0 (Shear)𝑀 𝑥𝑏 𝑤′′ 𝑥𝑏 0 (Moment)Energy form1𝐹𝑤 2𝐿𝐸𝐼 𝑤 ′′0Goals Introducing higher order derivatives Non convex energy form when 𝑃 𝑃𝑐𝑟2 𝑃 𝑤′2 𝑑𝑥
BifurcationFEM discretization:We enforce continuity of slope on element boundary using Hermite cubic interpolation(𝑤𝑖 , 𝑤𝑖′ )(𝑤𝑗 , 𝑤𝑗′ )At each node there are 2 DOF’s : (𝑤𝑖 , 𝑤𝑖′ )Node 𝑗Node 𝑖ElementNote: Having 2 DOF’s per node means we need to construct new nodal and element graphsNondimensionalized form:2𝜋𝐿3 𝐹 1 𝐸𝐼21𝑤𝑒0′′𝑙𝑒2 𝑃𝑐𝐿2𝑤′ 2𝑑𝑧withPL2𝑃𝑐 EI
BifurcationNodal and Element graphs:𝐻 2𝐽 1𝐽 19 Energyminimizing states 2(bilateral symmetry)Each minimizing state of the nodal graph is mapped to one of the 9 solutions of the node{𝑤𝑐𝑖 𝑟, 𝑤𝑐𝑖 , 𝑤𝑐𝑖 𝑟} {𝑤′𝑖𝑐 𝑟, 𝑤′𝑖𝑐 , 𝑤′𝑖𝑐 𝑟}Element graph is constructed as a complete bipartite graph between consecutive nodes.Total connections (element graph) 81Exactly solvable weights!!!Total possible states for the element 81
Why is non-convexity a problem?2-element problem examplePSymmetry of the problem:𝑤1′ 𝑤3′𝑤2′ 0Essentially two degrees of freedom:𝑤1′ 𝑤𝑝 , 𝑤2 𝑤Observe single slack variable box algorithm cannot resolve thedownward hill of the saddle and gives a false stable point. Weaugment another slack variable for the slope.Red region constitutes the downward hillof the saddle
Multiple slack variablesA heuristic remedy:Naturally identifiescritical loads via energyminimizationConsider different box sizes for different variables𝑖{𝑤𝑐𝑖 𝑟1 , 𝑤𝑐𝑖 , 𝑤𝑐𝑖 𝑟1 } {𝑤 ′ 𝑐 𝑟2 , 𝑤′𝑖𝑐 , 𝑤′𝑖𝑐 𝑟2 }In the post buckling solution (for 2 element case) thesolution tends to choose the up/down solution with similarlikelihoods.More work is needed in this direction to extend the algorithm for non-convex problems.
Summary Formulated and implemented an energy-based algorithm for solving differential equations on quantumannealer. Showed applications in following different types of equations:– Truss mechanics Convergence of the method for conv
Solving differential equations on quantum computers Prof. Veera Sundararaghavan Department of Aerospace Engineering, University of Michigan Sid Srivastava (PhD candidate) Keynote Talk: Modeling and Computation session 16th Pan-American Congress of Applied Mechanics May 23, 2019 Acknowledgments: USRA Quantum information Sciences Program
9.1 Properties of Radicals 9.2 Solving Quadratic Equations by Graphing 9.3 Solving Quadratic Equations Using Square Roots 9.4 Solving Quadratic Equations by Completing the Square 9.5 Solving Quadratic Equations Using the Quadratic Formula 9.6 Solving Nonlinear Systems of Equations 9 Solving Quadratic Equations
(iii) introductory differential equations. Familiarity with the following topics is especially desirable: From basic differential equations: separable differential equations and separa-tion of variables; and solving linear, constant-coefficient differential equations using characteristic equations.
Lesson 2a. Solving Quadratic Equations by Extracting Square Roots Lesson 2b. Solving Quadratic Equations by Factoring Lesson 2c. Solving Quadratic Equations by Completing the Square Lesson 2d. Solving Quadratic Equations by Using the Quadratic Formula What I Know This part will assess your prior knowledge of solving quadratic equations
CONCEPT IN SOLVING TRIG EQUATIONS. To solve a trig equation, transform it into one or many basic trig equations. Solving trig equations finally results in solving 4 types of basic trig equations, or similar. SOLVING BASIC TRIG EQUATIONS. There are 4 types of common basic trig equations: sin x a cos x a (a is a given number) tan x a cot x a
1 Solving Linear Equations 1.1 Solving Simple Equations 1.2 Solving Multi-Step Equations 1.3 Solving Equations with Variables on Both Sides 1.4 Rewriting Equations and Formulas Mathematical Thinking: Mathematically proficient students can apply the mathematics they know to solve problems arising in ev
Introduction to Advanced Numerical Differential Equation Solving in Mathematica Overview The Mathematica function NDSolve is a general numerical differential equation solver. It can handle a wide range of ordinary differential equations (ODEs) as well as some partial differential equations (PDEs). In a system of ordinary differential equations there can be any number of
Introduction to Differential Equation Solving with DSolve The Mathematica function DSolve finds symbolic solutions to differential equations. (The Mathe- matica function NDSolve, on the other hand, is a general numerical differential equation solver.) DSolve can handle the following types of equations: † Ordinary Differential Equations
Andhra Pradesh State Council of Higher Education w.e.f. 2015-16 (Revised in April, 2016) B.A./B.Sc. FIRST YEAR MATHEMATICS SYLLABUS SEMESTER –I, PAPER - 1 DIFFERENTIAL EQUATIONS 60 Hrs UNIT – I (12 Hours), Differential Equations of first order and first degree : Linear Differential Equations; Differential Equations Reducible to Linear Form; Exact Differential Equations; Integrating Factors .
