Quantum algorithm for linear differential equations with exponentially improved dependence on precision

0.21 MB·

by  Dominic W. Berry

Quantum algorithm for linear differential equations with
exponentially improved dependence on precision
DominicW. Berry ∗ Andrew M.Childs†,‡ Aaron Ostrander‡,§ GuomingWang ‡
7
1
0
2

b Abstract
e
F Wepresentaquantumalgorithmforsystemsof(possiblyinhomogeneous)linearordinarydifferentialequations
withconstantcoefficients. Thealgorithmproducesaquantumstatethatisproportional tothesolutionatadesired
7
final time. The complexity of the algorithm is polynomial in the logarithm of the inverse error, an exponential
1
improvementoverpreviousquantumalgorithmsforthisproblem.Ourresultbuildsuponrecentadvancesinquantum

linearsystemsalgorithmsbyencodingthesimulationintoasparse,well-conditionedlinearsystemthatapproximates
]
h evolutionaccordingtothepropagatorusingaTaylorseries.Unlikewithfinitedifferencemethods,ourapproachdoes
p notrequireadditionalhypothesestoensurenumericalstability.
-
t
n
a 1 Introduction
u
q
OneoftheoriginalmotivationsfordevelopingaquantumcomputerwastoefficientlysimulateHamiltoniandynamics,
[
i.e.,differentialequationsoftheform d~x =A~xwhereAisanti-Hermitian.GivenasuitabledescriptionofA,acopyof
dt
2 theinitialquantumstate x(0) ,andanevolutiontimeT,thegoalistoproduceaquantumstatethatisε-closetothe
v | i
finalstate x(T) . Thefirstalgorithmsforthisproblemhadcomplexitypolynomialin1/ε[13,2,10,4]. Subsequent
4 | i
8 workgaveanalgorithmwithcomplexitypoly(log(1/ε))—anexponentialimprovement—whichisoptimalinablack-
6 boxmodel[5].Morerecentworkhasstreamlinedthesealgorithmsandimprovedtheirdependenceonotherparameters
3 [6,7,15,8,14,16].
0 WhileHamiltoniansimulationhasbeenafocusofquantumalgorithmsresearch,themoregeneralproblemofsim-
1. ulatinglineardifferentialequationsoftheform d~x =A~x+~bforarbitrary(sparse)Aislesswellstudied.Reference[3]
dt
0 solvesthisproblemusingaquantumlinearsystemsalgorithm(QLSA)toimplementlinearmultistepmethods,which
7
representthedifferentialequationswithasystemoflinearequationsbydiscretizingtime. Thecomplexityofthisap-
1
proachispoly(1/ε). ConsideringtherecentimprovementstothecomplexityofHamiltoniansimulation,itisnatural
:
v toaskwhetherlineardifferentialequationscanbesolvedmoreefficientlyasafunctionofε.
i
X HamiltoniansimulationisacentralcomponentoftheQLSA,andthetechniquesunderlyingpoly(log(1/ε))Hamil-
tonian simulation have been adapted to give a QLSA with complexity poly(log(1/ε)) [11]. However, even if this
r
a improvedQLSA is used to implementthe algorithmof Ref. [3], the overallcomplexityis still poly(1/ε), since the
multistepmethoditselfisasignificantsourceoferror.
In a similar vein, the QLSA of Ref. [12] has poly(1/ε) complexity even when using a Hamiltonian simulation
algorithmwithpoly(log(1/ε))complexity,simplybecausephaseestimationhascomplexitypoly(1/ε).Reference[11]
provides a QLSA with poly(log(1/ε)) complexity by avoiding phase estimation and instead directly inverting the
linearsystemusingalinearcombinationofunitaries(LCU).Thus,onemightconsiderrealizingthesolutionof d~x =
dt
A~x+~basalinearcombinationofunitaries. Unfortunately,inthegeneralcasewhereAisnotanti-Hermitian,thebest
implementationofthisapproachthatweareawareofhasanexponentiallysmallsuccessprobability.
∗DepartmentofPhysicsandAstronomy,MacquarieUniversity
†DepartmentofComputerScienceandInstituteforAdvancedComputerStudies,UniversityofMaryland
‡JointCenterforQuantumInformationandComputerScience,UniversityofMaryland
§DepartmentofPhysics,UniversityofMaryland
1
In this paper, we circumventthese limitations and present a quantum algorithm for linear differentialequations
withcomplexitypoly(log(1/ε)),anexponentialimprovementoverRef.[3]. AsinRef.[3],ourapproachappliesthe
QLSA.However,insteadofusingalinearmultistepmethod,weencodeatruncationoftheTaylorseriesofexp(At),
thepropagatorforthedifferentialequation,intoalinearsystem. Sinceiteffectivelyimplementsalinearcombination
ofoperations,ourapproachisconceptuallysimilartoquantumsimulationvialinearcombinationsofunitaries,butwe
achievesignificantlybetterperformancebyconstructingthislinearcombinationstepwisethroughasystemoflinear
equations.ThisalternativetodirectapplicationofLCUmethodsmightbeadvantageousforotherquantumalgorithms.
In addition to scaling well with the simulation error, our algorithm has favorable performance as a function of
other parameters. The complexity is nearly linear in the evolution time, which is a quadratic improvement over
Ref. [3] and is nearly optimal [4]. The complexity is also nearly linear in the sparsity of A and in a parameter
characterizing the decay of the solution vector. The latter dependence is necessary since producing a normalized
versionofasubnormalizedsolutionvectorisequivalenttopostselection,whichiscomputationallyintractable[1],as
discussed furtherin Section8. Alongsimilar lines, we assume thattheeigenvaluesof A havenon-positiverealpart
sinceitisintractabletosimulateexponentiallygrowingsolutions.(ThisimprovesuponRef.[3],wheretheeigenvalues
λ of A must satisfy arg( λ) α for some constant α dependingon the stability of the multistep method.) For a
| − |≤
precisestatementofthemainresult,seeTheorem9.
This paper is organizedas follows. Section 2 describes how we encode the solution of a system of differential
equationsintoasystemoflinearequations. Thefollowingthreesectionsanalyzepropertiesofthissystem: Section3
bounds its condition number, Section 4 analyzes how well it approximates the differential equation, and Section 5
shows that a measurement of its solution vector provides a solution of the differential equation with appreciable
probability. In Section 6 we explain how to prepare the state that is input to the QLSA using black boxes for the
initialconditionandinhomogeneoustermofthedifferentialequation. Weformallystateandproveourmainresultin
Section7. Finally,weconcludeinSection8withadiscussionoftheresultandsomeopenproblems.
2 Constructing the Linear System
AsinRef.[3]weconsideradifferentialequationoftheform
d~x
=A~x+~b (1)
dt
whereAand~baretime-independent.Thishastheexactsolution
~x(t)=exp(At)~x(0)+(exp(At) I)A 1~b. (2)
−
−
Define
k zj
T (z):= ∑ exp(z) (3)
k
j! ≈
j=0
and
k zj 1
S (z):= ∑ − (exp(z) 1)z 1 (4)
k −
j! ≈ −
j=1
wheretheapproximationsholdforlargek. Thenforshortevolutiontimeh(namelyh 1/ A ,where denotesthe
≤ k k k·k
spectralnorm)andlargek,wecanapproximatethesolutionby
~x(h) T (Ah)~x(0)+S (Ah)h~b. (5)
k k
≈
Thisapproximatesolutioncanbeusedinturnasaninitialconditionforanotherstepofevolution,andwecanrepeat
thisprocedureasdesiredforatotalnumberofstepsm.
Weencodethisprocedureinalinearsystemusingthefollowingfamilyofmatrices.
2
Definition1. LetAbeanN Nmatrix,andletm,k,p Z+. Define
× ∈
d m 1 k
C (A):= ∑ j j I ∑− ∑ i(k+1)+j i(k+1)+j 1 A/j
m,k,p
| ih |⊗ − | ih − |⊗
j=0 i=0 j=1
m 1 k d
∑− ∑ (i+1)(k+1) i(k+1)+j I ∑ j j 1 I, (6)
− | ih |⊗ − | ih − |⊗
i=0 j=0 j=d p+1
−
whered:=m(k+1)+p,andIistheN Nidentitymatrix.
×
Nowconsiderthelinearsystem
m 1
C (Ah)x = 0 x +h ∑− i(k+1)+1 b , (7)
m,k,p in
| i | i| i | i| i
i=0
where x , b CN andh R+. ThefirstregisterlabelsanaturalblockstructureforC . Forexample,thesystem
in m,k,p
C2,3,2(A|h)i|xi| =i∈|0i|xini+h∑∈1i=0|4i+1i|biisasfollows:
I x
in
| i
Ah I hb
−   | i
Ah/2 I 0
−
 Ah/3 I   0 
 −   
 I I I I I   0 
 − − − −   
C2,3,2(Ah)x = Ah I x =hb . (8)
| i  − | i  | i
 Ah/2 I   0 
 −   
 Ah/3 I   0 
 −   
 I I I I I   0 
 − − − −   
 I I   0 
 −   
 I I  0 
 −   
   
AfterperformingmstepsoftheevolutionapproximatedwithaTaylorseriesoforderk,thesolutioniskeptconstant
for psteps. Thisensuresasignificantprobabilityofobtainingthesolutionatthefinaltime, similarlyasinRef.[3].
NotethatC (Ah)isnonsingular,sinceitisalower-triangularmatrixwithnonzerodiagonalentries. IfN=1,then
m,k,p
C (Ah)isa(d+1) (d+1)matrix.
m,k,p
×
ThesolutionofEq.(7)is
m 1
x =C (Ah) 1 0 x +h ∑− i(k+1)+1 b , (9)
m,k,p − in
| i "| i| i i=0| i| i#
whichcanbewrittenas
m 1 k p
x = ∑− ∑ i(k+1)+j x +∑ m(k+1)+j x (10)
i,j m,j
| i | i | i
i=0 j=0 j=0
(cid:12) (cid:11) (cid:12) (cid:11)
forsome x CN. BythedefinitionofC (Ah),the(cid:12)se x ssatisfy (cid:12)
i,j m,k,p i,j
∈
(cid:12)(cid:12) (cid:11) x0,0 = xin , (cid:12)(cid:12) (cid:11) (11)
| i | i
k
x = ∑ x , 1 i m, (12)
i,0 i 1,j
| i − ≤ ≤
j=0
(cid:12) (cid:11)
x =Ah(cid:12)x +hb , 0 i<m, (13)
i,1 i,0
| i | i | i ≤
x =(Ah/j) x , 0 i<m, 2 j k, (14)
i,j i,j 1
− ≤ ≤ ≤
(cid:12)xm,j(cid:11)= xm,j 1(cid:12), (cid:11) 1 j p. (15)
(cid:12) − (cid:12) ≤ ≤
(cid:12) (cid:11) (cid:12) (cid:11)
(cid:12) (cid:12)
3
Fromtheseequations,weobtain
x = x , (16)
0,0 in
| i | i
x =((Ah)j/j!)x +((Ah)j 1/j!)hb , 1 j k, (17)
0,j 0,0 −
| i | i ≤ ≤
x =T (Ah)x +S (Ah)hb
(cid:12) 1,0(cid:11) k 0,0 k
(cid:12)| i | i | i
exp(Ah)x +(exp(Ah) I)A 1 b , (18)
in −
≈ | i − | i
x1,j =((Ah)j/j!)x1,0 +((Ah)j−1/j!)hb , 1 j k, (19)
| i | i ≤ ≤
x =T (Ah)x +S (Ah)hb
(cid:12) 2,0(cid:11) k 1,0 k
(cid:12)| i | i | i
exp(2Ah)xin +(exp(2Ah) I)A−1 b , (20)
≈ | i − | i
.
.
.
x =T (Ah)x +S (Ah)hb
m 1,0 k m 2,0 k
| − i | − i | i
exp(Ah(m 1))x +(exp(Ah(m 1)) I)A 1 b , (21)
in −
≈ − | i − − | i
x =((Ah)j/j!)x +((Ah)j 1/j!)hb , 1 j k, (22)
m 1,j m 1,0 −
− | − i | i ≤ ≤
x =T (Ah)x +S (Ah)hb
(cid:12) m,0(cid:11) k m 1,0 k
(cid:12) | i | − i | i
exp(Ahm)x +(exp(Ahm) I)A 1 b , (23)
in −
≈ | i − | i
x = x
m,j m,0
| i
(cid:12) (cid:11) exp(Ahm)xin +(exp(Ahm) I)A−1 b , 1 j p. (24)
(cid:12) ≈ | i − | i ≤ ≤
Intheseapproximations,weassumekissufficientlylargethatwecanneglectthetruncationerrors T (Ah) exp(Ah)
k
and S (Ah)h (exp(Ah) I)A 1 (we make this more precise in Section 4). Note that x (dkefined b−y Eq. (10)k)
k −
− − | i
includesapiecethatcanbeinterpretedasthehistorystateoftheevolution
(cid:13) (cid:13)
(cid:13) (cid:13)
d~x
=A~x+~b (25)
dt
withtheinitialcondition~x(0)=~x . Moreprecisely, x isagoodapproximationofthesystem’sstateattimeih,for
in i,0
| i
anyi 0,1,...,m . Furthermore, x = x = = x isagoodapproximationof
m,0 m,1 m,p
∈{ } | i | i ··· | i
~x(t)=exp(At)~x +(exp(At) I)A 1~b (26)
in −
−
fort=mh.
Ifwemeasuredthefirstregisterofthestate x / x inthestandardbasis,wewouldobtainthestate x / x
i,j i,j
| i k| ik
for random i,j. By choosing a large p, we can ensure there is a high probability of obtaining x / x ,
m,0 m,0
| (cid:12) i(cid:11)k(cid:13)|(cid:12) i(cid:11)k(cid:13)
xm,1 / xm,1 , ..., or xm,p / xm,p (as we show in Section 5). Then this probability can be raise(cid:12)d to Ω(cid:13)((cid:12)1) b(cid:13)y
| i k| ik | i k| ik
using amplitude amplification (or by classical repetition). This is how we prepare a state close to~x(t)/ ~x(t) for
k k
t=mh.
NotethatthestatewegenerateisnotofthesameformasthehistorystateinRef.[3],whichonlyencodes~x(t)at
intermediatetimes. Thesolutionofourlinearsystemnotonlyencodes~x(t)atintermediatetimes(via x )butalso
i,0
| i
encodes((Ah)j/j!)~x(t)+((Ah)j 1/j!)h~batintermediatetimes(via x ).
− i,j
| i
To analyze the performance of this approach to solving differential equations, we establish three properties of
thissystem oflinearequations. First, since the complexityofthe bestknownQLSAsgrowslinearlywithcondition
number(andsublinearcomplexityisimpossibleunlessBQP=PSPACE[12]), weanalyzetheconditionnumberof
C (Section3). Second,weshowthatthesolutionofthelinearsystemincludesapiecethatisclosetothesolution
m,k,p
oftheassociateddifferentialequation(Section4).Third,weshowthatthispiececanbeobtainedfromameasurement
thatsucceedswithappreciableprobability(Section5).
4
3 Condition Number
In this section, we upperboundthe conditionnumberof the matrixC (A) undermild assumptionsaboutA. We
m,k,p
beginwithatechnicallemmathatupperboundsthenormsofthecolumnsoftheinverseofthismatrix.
Lemma2. Letλ Csuchthat λ 1andRe(λ) 0. Letm,k,p Z+ suchthatk 5and(k+1)! 2m,andlet
∈ | |≤ ≤ ∈ ≥ ≥
d=m(k+1)+p.Thenforanyn,l 0,1,...,d ,
∈{ }
C (λ) 1 l 1.04eI (2)(m+p) (27)
m,k,p − 0
| i ≤
withI0(2)<2.28amodifiedBesselfun(cid:13)(cid:13)ctionofthefirs(cid:13)(cid:13)tkinpd,and
nCm,k,p(λ)−1 l √1.04e. (28)
h | | i ≤
(cid:12) (cid:12)
(cid:12) (cid:12)
Proof. RecallthedefinitionsofT (z)andS (z)inEqs.(3)and(4),respectively.Wealsodefine
k k
k b!zj b
T (z):= ∑ − (29)
b,k j!
j=b
forb k.
≤
Fixanyl 0,1,...,d . Supposethesolutionofthelinearsystem
∈{ }
C (λ)x = l (30)
m,k,p
| i | i
is
m 1 k p
x = ∑− ∑x i(k+1)+j +∑x m(k+1)+j (31)
i,j m,j
| i | i | i
i=0 j=0 j=0
forsomex C. BythedefinitionofC (λ),thex sshouldsatisfy
i,j m,k,p i,j
∈
k
x ∑x =δ , 1 i m, (32)
i,0 i 1,j i(k+1),l
− − ≤ ≤
j=0
x (λ/j)x =δ , 0 i<m, 1 j k, (33)
i,j i,j 1 i(k+1)+j,l
− − ≤ ≤ ≤
x x =δ , 1 j p, (34)
m,j m,j 1 m(k+1)+j,l
− − ≤ ≤
whereδ =1ifi= j,and0otherwise.
i,j
Weconsiderthecases0 l<m(k+1)andm(k+1) l dseparately.
≤ ≤ ≤
Case1: 0 l<m(k+1). Supposel=a(k+1)+bforsome0 a<mand0 b k. Inthiscase,Eq.(34)
• ≤ ≤ ≤ ≤
implies
x =0, 0 i<a, 0 j k, (35)
i,j
≤ ≤ ≤
x =0, 0 j<b, (36)
a,j
≤
x =b!λj b/j!, b j k, (37)
a,j −
≤ ≤
x =T (λ), (38)
a+1,0 b,k
x =(λj/j!)x , 1 j k, (39)
a+1,j a+1,0
≤ ≤
x =T (λ)x =T (λ)T (λ), (40)
a+2,0 k a+1,0 k b,k
.
.
.
x =T (λ)x =(T (λ))m a 1T (λ), (41)
m,0 k m 1,0 k − − b,k
−
5
x =x =(T (λ))m a 1T (λ), 1 j p. (42)
m,j m,0 k − − b,k
≤ ≤
Since λ 1,foranyb j k,wehave
| |≤ ≤ ≤
xa,j =b!λ j−b/j! b!/j! 1. (43)
| | ≤ ≤
(cid:12) (cid:12)
Furthermore,since λ 1andRe(λ) (cid:12)0,b(cid:12)yLemma10andLemma11inAppendixA,wehave
| |≤ ≤
1 1
T (λ) 1+ 1+ , (44)
k
| |≤ (k+1)! ≤ 2m
T (λ) √1.04. (45)
b,k
≤
Consequently,wehave (cid:12) (cid:12)
(cid:12) (cid:12)
x = T (λ)i a 1T (λ) (1+1/2m)m√1.04 √1.04e, a+1 i m, (46)
i,0 k − − b,k
| | ≤ ≤ ≤ ≤
(cid:12) (cid:12)
and (cid:12) (cid:12)
x = (λj/j!)x x √1.04e, a+1 i m, 1 j k. (47)
i,j i,0 i,0
≤| |≤ ≤ ≤ ≤ ≤
Itfollowsthat (cid:12) (cid:12) (cid:12) (cid:12)
(cid:12) (cid:12) (cid:12) (cid:12)
x = x √1.04e, 0 j p. (48)
m,j m,0
| |≤ ≤ ≤
(cid:12) (cid:12)
Usingthesefacts,weobtain (cid:12) (cid:12)
m 1 k p
x 2= ∑− ∑ x 2+∑ x 2
i,j m,j
k| ik
i=0 j=0 j=0
(cid:12) (cid:12) (cid:12) (cid:12)
= ∑k b!λ(cid:12)j b/(cid:12)j! 2+(cid:12)m∑−1(cid:12)∑k (λj/j!)x 2+(p+1)x 2
− i,0 m,0
| |
j=b i=a+1j=0
(cid:12) (cid:12) (cid:12) (cid:12)
∑k ((cid:12)(cid:12)b!/j!)2+(cid:12)(cid:12)m∑−1 ∑k (1/j(cid:12)!)2 x 2+((cid:12)p+1)x 2
i,0 m,0
≤ | | | |
j=b i=a+1j=0
k m 1 k
∑(b!/j!)2+1.04e ∑− ∑(1/j!)2+1.04e(p+1)
≤
j=b i=a+1j=0
I (2)+1.04eI (2)(m a 1)+1.04e(p+1)
0 0
≤ − −
1.04eI (2)(m+p), (49)
0
≤
whereinthefifthstepweusethefacts
k ∞
∑(1/j!)2 ∑(1/j!)2=I (2), (50)
0
≤
j=0 j=0
k k b 1 1 4
∑(b!/j!)2 ∑− 1/(b+1)2s =1+ <I (2), 1 b k. (51)
≤ ≤ 1 (b+1) 2 b(b+2) ≤ 3 0 ≤ ≤
j=b s=0 − −
Case2:m(k+1) l d. Supposel=m(k+1)+bforsome0 b p. Inthiscase,Eq.(34)implies
• ≤ ≤ ≤ ≤
x =0, 0 i<m, 0 j k, (52)
i,j
≤ ≤ ≤
x =0, 0 j<b, (53)
m,j
≤
x =1, b j p. (54)
m,j
≤ ≤
6
Itfollowsthat
m 1 k p
x 2= ∑− ∑ x 2+∑ x 2
i,j m,j
k| ik
i=0 j=0 j=0
(cid:12) (cid:12) (cid:12) (cid:12)
=p b+(cid:12)1 (cid:12) (cid:12) (cid:12)
−
p+1. (55)
≤
Inbothoftheabovecases,wehave x 1.04eI (2)(m+p)and nx √1.04eforanyn 0,1,...,d ,
0
k| ik≤ |h | i|≤ ∈{ }
asclaimed.
p
Nowwearereadytoupperboundthenormoftheinverseofthematrix.
Lemma 3. Let A =VDV 1 be a diagonalizable matrix, where D = diag(λ ,λ ,...,λ ) satisfies λ 1 and
− 0 1 N 1 i
Re(λ) 0fori 0,1,...,N 1 .Letm,k,p Z+ suchthatk 5and(k+1)! 2m. Th−en | |≤
i
≤ ∈{ − } ∈ ≥ ≥
C (A) 1 3κ √k(m+p), (56)
m,k,p − V
≤
whereκV = V V−1 istheconditionn(cid:13)(cid:13)umberofV.(cid:13)(cid:13)
k k·
Proof. Forconven(cid:13)ience(cid:13),wewilldropthesubscriptsm,k,p,anduseC()todenoteC (). WediagonalizeC(A)as
(cid:13) (cid:13) · m,k,p ·
C(A)=V˜C(D)V˜ 1, (57)
−
where V˜ := ∑dj=0|jihj|⊗V has condition number κV˜ = κV. Then we may upper bound C(A)−1 in terms of
C(D) 1 as
− (cid:13) (cid:13)
(cid:13) (cid:13)
(cid:13)(cid:13) (cid:13)(cid:13) C(A)−1 = V˜C(D)−1V˜−1
(cid:13) (cid:13) (cid:13)V˜ C(D)−1(cid:13) V˜−1
(cid:13) (cid:13)≤(cid:13) · (cid:13)·
=(cid:13)(cid:13)C((cid:13)(cid:13)D)(cid:13)(cid:13)−1 ·κV˜.(cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13) (58)
(cid:13) (cid:13)
Wehave (cid:13) (cid:13)
C(D) 1 ψ
C(D) 1 =max − | i (59)
−
(cid:13) (cid:13) |ψi (cid:13)(cid:13) k|ψik (cid:13)(cid:13)
wherewemaximizeoverallstates|ψi∈C((cid:13)d+1)NTh(cid:13)estate|ψicanbewrittenas|ψi=∑dl=0|li|ψliforsome|ψli∈CN.
Thenwehave
2
d
C(D) 1 ψ 2= ∑C(D) 1 l ψ
− − l
| i (cid:13)l=0 | i| i(cid:13)
(cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13)(cid:13)(cid:13)(d+1)∑d C(D)−1(cid:13)(cid:13)(cid:13)(cid:13)l ψl 2. (60)
≤ | i| i
l=0
(cid:13) (cid:13)
(cid:13) (cid:13)
Nowlet|ψli=∑Nj=−01ψj,l|jiforsomeψj,l∈C.WithD=∑Nj=−01λj|jihj|,wehaveC(D)=∑Nj=−01C(λj)⊗|jihj|.Hence,
byLemma2,weobtain
2
N 1
C(D) 1 l ψ 2= ∑− ψ C(λ ) 1 l j
− l j,l j −
| i| i (cid:13)j=0 | i| i(cid:13)
(cid:13)(cid:13) (cid:13)(cid:13) =(cid:13)(cid:13)(cid:13)(cid:13)N∑−1 ψj,l 2 C(λj)−1 l (cid:13)(cid:13)(cid:13)(cid:13)2
| i
j=0
(cid:12) (cid:12) (cid:13) (cid:13)
(cid:12) (cid:12) (cid:13) (cid:13)
7
N 1
1.04eI (2)(m+p) ∑− ψ 2
0 j,l
≤
j=0
(cid:12) (cid:12)
=1.04eI (2)(m+p) ψ(cid:12) 2.(cid:12) (61)
0 l
k| ik
UsingthisexpressioninEq.(60)gives
d
C(D) 1 ψ 2 1.04eI (2)(d+1)(m+p)∑ ψ 2
− 0 l
| i ≤ k| ik
l=0
(cid:13) (cid:13)
(cid:13) (cid:13) =1.04eI (2)(m(k+1)+p+1)(m+p) ψ 2
0
k| ik
6
1.04eI (2)k(m+p)2 ψ 2. (62)
0
≤ 5× k| ik
UsingthisresultinEq.(59)yields
C(D) 1 3√k(m+p). (63)
−
≤
CombiningthisexpressionwithEq.(58)then(cid:13)givesEq(cid:13).(56),asclaimed.
(cid:13) (cid:13)
Itremainstoupperboundthenormofthematrix.
Lemma4. LetAbeanN Nmatrixsuchthat A 1. Letm,k,p Z+,andk 5. Then
× k k≤ ∈ ≥
C (A) 2√k. (64)
m,k,p
≤
(cid:13) (cid:13)
(cid:13) (cid:13)
Proof. ObservethatC:=C (A)canbewrittenasthesumofthreematrices:
m,k,p
C=C +C +C (65)
1 2 3
where
d
C := ∑ j j I, (66)
1
| ih |⊗
j=0
m 1 k
C := ∑− ∑ (i+1)(k+1) i(k+1)+j I, (67)
2
− | ih |⊗
i=0 j=0
m 1 k d
C := ∑− ∑ i(k+1)+j i(k+1)+j 1 A/j ∑ j j 1 I, (68)
3
− | ih − |⊗ − | ih − |⊗
i=0 j=1 j=d p+1
−
whered=m(k+1)+p. Onecaneasilycheckthat C =1, C =√k+1,and C =max A ,1 =1(thisis
1 2 3
k k k k k k {k k }
trivialforC ,andfollowsdirectlyfromacalculationofC C†andC C† fortheothercases). Consequently,
1 2 2 3 3
C C + C + C
1 2 3
k k≤k k k k k k
√k+1+2
≤
2√k (69)
≤
asclaimed.
CombiningLemma3andLemma4,weobtainthefollowingupperboundontheconditionnumberofC (A):
m,k,p
Theorem5. LetA=VDV 1beadiagonalizablematrixsuchthat A 1,D=diag(λ ,λ ,...,λ )andRe(λ) 0,
− 0 1 N 1 i
for i 0,1,...,N 1 . Let m,k,p Z+ such that k 5 andk(kk+≤1)! 2m. Let C :=C −(A), and let κ≤=
m,k,p C
C ∈C{ 1 bethe−con}ditionnumber∈ofC. Then ≥ ≥
−
k k·
(cid:13) (cid:13) κ 6κ k(m+p), (70)
(cid:13) (cid:13) C≤ V
whereκ = V V 1 istheconditionnumberofV.
V −
k k·
(cid:13) (cid:13)
(cid:13) (cid:13)
8
4 Solution Error
Inthissection,weprovethatthesolutionofthelinearsystemdefinedbyEq.(7)encodesagoodapproximationofthe
solutionofthedifferentialequationdefinedbyEq.(25)withtheinitialcondition~x(0)=~x .
in
Theorem6. LetA=VDV 1 beadiagonalizablematrix, where D=diag(λ ,λ ,...,λ )satisfiesRe(λ) 0for
− 0 1 N 1 i
i 0,1,...,N 1 . Leth R+ suchthat Ah 1. Let x , b CN, andlet x(t) −be definedby Eq. (2≤6). Let
in
m∈,k{,p Z+suc−hth}atk 5a∈nd(k+1)! 2mk.Lekt≤x bed|efinied|biy∈Eqs.(9)and(1|0).Tihenforany j 0,1,...,m ,
i,j
∈ ≥ ≥ ∈{ }
x(jh) x (cid:12)2.8κ(cid:11) j( x +mh b )/(k+1)!, (71)
j,0 (cid:12) V in
| i− ≤ k| ik k| ik
whereκV = V V−1 isth(cid:13)(cid:13)econdition(cid:12)(cid:12)num(cid:11)b(cid:13)(cid:13)erofV.
k k·
Proof. Notethat x(cid:13)(jh)(cid:13)(thesolutionofthedifferentialequation)satisfiestherecurrencerelation
|(cid:13) i(cid:13)
x((j+1)h) =exp(Ah)x(jh) +(exp(Ah) I)A−1 b , (72)
| i | i − | i
while x (inthesolutionoftheassociatedlinearsystem)satisfiestherecurrencerelation
j,0
(cid:12) (cid:11) x =T (Ah) x +S (Ah)hb . (73)
(cid:12) j+1,0 k j,0 k
| i
RecallthatTk(λ)=∑kj=0λj!j andSk(λ)=(cid:12)(cid:12)∑kj=1(cid:11)λjj−!1. Inadd(cid:12)(cid:12)tion(cid:11),wehave|x(0)i=|x0,0i=|xini.
Define y(t) :=V 1 x(t) and y =V 1 x . We will give an upperboundon δ := y(jh) y and
− i,j − i,j j j,0
| i | i | i−
convertit into an upper boundon ε := x(jh) x . Since A=VDV 1, we have exp(Ah)=Vexp(Dh)V 1,
(cid:12)j (cid:11) | (cid:12)i−(cid:11) j,0 − (cid:13) (cid:12) (cid:11)(cid:13) −
T (Ah)=VT (Dh)V 1,andS (Ah)(cid:12)=VS (Dh)V(cid:12) 1. ThenEq.(72)implies (cid:13) (cid:12) (cid:13)
k k − k (cid:13)k − (cid:12) (cid:11)(cid:13)
(cid:13) (cid:12) (cid:13)
y((j+1)h) =exp(Dh)y(jh) +(exp(Dh) I)D 1 c , (74)
−
| i | i − | i
where c =V 1 b . Meanwhile,Eq.(73)implies
−
| i | i
y =T (Dh) y +S (Dh)hc . (75)
j+1,0 k j,0 k
| i
Inaddition,wehave y(0) = y0,0 = y(cid:12)(cid:12)in :=V(cid:11)−1 xin . (cid:12)(cid:12) (cid:11)
| i | i | i | i
Now since Re(λ) 0 for any i 0,1,...,N 1 , we have exp(Dh) 1. Moreover, since Ah 1, we
i
≤ ∈{ − } k k≤ k k≤
have λh 1 for any i 0,1,...,N 1 . Then since Re(λh) 0 for any i 0,1,...,N 1 , by Lemma 10 in
i i
| |≤ ∈{ − } ≤ ∈{ − }
AppendixA,weget
exp(Dh) T(Dh) 1/(k+1)!, (76)
k
k − k≤
andbyLemma12inAppendixA,weget
S (Dh) (exp(Dh) I)D 1h 1 1/(k+1)!. (77)
k − −
− − ≤
Theerrorδ = y(jh) y (cid:13) canbeboundedasfollows. N(cid:13)otethatδ =0,andusingthetriangleinequality,
j j,0(cid:13) (cid:13) 0
| i−
wehave
(cid:13) (cid:12) (cid:11)(cid:13)
(cid:13) (cid:12) (cid:13)
δj+1= exp(Dh)y(jh) +(exp(Dh) I)D−1 c Tk(Dh) yj,0 Sk(Dh)hc
| i − | i− − | i
(cid:13)exp(Dh)y(jh) Tk(Dh) yj,0 + (exp(Dh) I(cid:12))D−(cid:11)1 c Sk(Dh)h(cid:13)c
≤(cid:13) | i− − (cid:12) | i− (cid:13)| i
(cid:13)exp(Dh)y(jh) exp(Dh(cid:12)) yj,(cid:11)0(cid:13) +(cid:13) exp(Dh)yj,0 Tk(Dh)yj,0 (cid:13)
≤(cid:13) | i− (cid:12) (cid:13) (cid:13) | i− | i (cid:13)
+(cid:13) (exp(Dh) I)D−1 c S(cid:12)k(Dh(cid:11))(cid:13)hc(cid:13) (cid:13)
(cid:13) − | i− (cid:12) (cid:13)| i(cid:13) (cid:13)
e(cid:13)xp(Dh) y(jh) yj,0 + exp(D(cid:13)h) Tk(Dh) yj,0
≤k (cid:13) k | i−| i k (cid:13) − k | i
+ (exp(Dh(cid:13)) I)D−1 Sk((cid:13)Dh)h c (cid:13) (cid:13)
(cid:13) − − (cid:13) k| ik (cid:13) (cid:13)
1
δ (cid:13)+ y +h c .(cid:13) (78)
≤ j(cid:13) (k+1)! j,0 k| ik (cid:13)
(cid:2)(cid:13)(cid:12) (cid:11)(cid:13) (cid:3)
(cid:13)(cid:12) (cid:13)
9
Thisimplies
j
δ max y +h c . (79)
j i,0
≤ (k+1)! 0 i jk| ik k| ik
(cid:20) ≤≤ (cid:21)
Next,wegiveanupperboundonmax y . Fixanyi 0,1,...,m . Thenwehave
0 i m i,0
≤≤ k| ik ∈{ }
yi,0 = i(k+1)Cm,k,p(D)−1 z , (80)
| i h | | i
wherethebra i(k+1) actsonthefirstregister,and
h |
m 1
z = 0 y +h ∑− j(k+1)+1 c . (81)
in
| i | i| i | i| i
j=0
Therefore,bythetriangleinequality,
m 1
yi,0 i(k+1)Cm,k,p(D)−1 0 yin +h ∑− i(k+1)Cm,k,p(D)−1 j(k+1)+1 c . (82)
k| ik≤ h | | i| i h | | i| i
j=0
(cid:13) (cid:13) (cid:13) (cid:13)
(cid:13) (cid:13) (cid:13) (cid:13)
To bound both terms, we consider the general expression i(k+1)C (D) 1 s β for integer s. Let β =
m,k,p −
h | | i| i | i
∑Nl=−01βl|liforsomeβl ∈C. ThensinceCm,k,p(D)=∑Nl=−01Cm(cid:13),k,p(λl)⊗|lihl|,weget (cid:13)
(cid:13) (cid:13)
N 1
i(k+1)Cm,k,p(D)−1 s β = ∑− βl i(k+1)Cm,k,p(λl)−1 s l . (83)
h | | i| i h | | i| i
l=0
Then,byLemma2,wehave
N 1
i(k+1)C (D) 1 s β 2= ∑− β 2 i(k+1)C (λ) 1 s 2
m,k,p − l m,k,p l −
h | | i| i | | h | | i
l=0
(cid:13) (cid:13) (cid:12) (cid:12)
(cid:13) (cid:13) N (cid:12)1 (cid:12)
1.04e ∑− β 2
l
≤ | |
l=0
=1.04e β 2. (84)
k| ik
Henceweobtain
i(k+1)C (D) 1 0 y √1.04e y , (85)
m,k,p − in in
h | | i| i ≤ k| ik
i(k+1)(cid:13)Cm,k,p(D)−1 j(k+1)+1 c (cid:13) √1.04e c . (86)
h (cid:13)| | i| i(cid:13)≤ k| ik
(cid:13) (cid:13)
Usingthesetwofactsandthe(cid:13)triangleinequality,Eq.(82)implies (cid:13)
y √1.04e( y +mh c ). (87)
i,0 in
k| ik≤ k| ik k| ik
Sincethisholdsforanyi 0,1,...,m ,weget
∈{ }
max y √1.04e( y +mh c ). (88)
i,0 in
0 i mk| ik≤ k| ik k| ik
≤≤
NowusingEqs.(79)and(88),weget
δ (√1.04e+1)j( y +mh c )/(k+1)!
j in
≤ k| ik k| ik
=(√1.04e+1)j V 1 x +mh V 1 b /(k+1)!
− in −
| i | i
(√1.04e+1)j(cid:0)(cid:13)V−1 ( xi(cid:13)n +m(cid:13)h b /(cid:13)((cid:1)k+1)!). (89)
≤ (cid:13) k| (cid:13)ik (cid:13) k| ik (cid:13)
(cid:13) (cid:13)
(cid:13) (cid:13)
10