Fast computation of power series solutions of systems of - - PowerPoint PPT Presentation

Fast computation of power series solutions

• f systems of differential equations

Alin Bostan ALGO, INRIA Rocquencourt joint work with

• F. Chyzak, F. Ollivier,
• B. Salvy, ´
• E. Schost, A. Sedoglavic

Context

find fast algorithms for solving (in power series) ordinary differential equations applications: combinatorics, numerical analysis (control), algebraic complexity ◮ fast means using few operations (±, ×, ÷) in the base field K.

More precisely

Given a linear differential equation with power series coefficients ar(t)y(r)(t) + · · · + a1(t)y′(t) + a0(t)y(t) = 0 compute the first N terms of a (basis of) power(s) series solution(s) y(t). ◮ naive algorithm (undetermined coefficients), O(rN 2). ◮ Best that can be hoped: complexity linear in N and polynomial in r. ◮ Same problem for linear systems Y ′ = AY and for non-linear systems.

Fast polynomial (matrix) multiplication

M(N) = complexity of polynomial multiplication in degree < N = O(N 2) by the na¨ ıve algorithm = O

• N 1.58

by Karatsuba’s algorithm = O

• N logα (2α−1)

by Toom-Cook algorithm = O

• N log(N) loglog(N)
• by Sch¨
• nhage–Strassen FFT

MM(r, N) = complexity of polynomial matrix mult., deg < N, size = r = O(rω M(N)) by the Cantor–Kaltofen algorithm = O(rωN + r2 M(N)) by the B.–Schost algorithm

Previous results

r = 1 Brent & Kung (1978): reduction to exponentials of power series +

Newton iteration, O(M(N)).

r = 2 Brent & Kung (1978): reduction to Riccati non-linear equations +

linearization, O(M(N)).

r > 1 Brent & Kung (1978), van der Hoeven (2002): order reduction +

linearization, O(rr M(N)).

r > 1 van der Hoeven (2002) 1 solution in O(r M(N) log N), model of

relaxed multiplication

New results

Problem constant polynomial power series

• utput

(input, output) coefficients coefficients coefficients size (eqn, basis) O(M(r)N) O(dr2N) O(MM(r, N)) O(rN) (eqn, 1 sol) O(M(r)N/r) O(drN) O(r M(N) log N) O(N) (sys, basis) O(rM(r)N) O(drωN) O(MM(r, N)) O(r2N) (sys, 1 sol) O(M(r)N) O(dr2N) O(r2 M(N) log N) O(rN)

Experimental results

"MatMul.dat" 2 3 4 5 6 7 8 9 10

• rder

1000 1500 2000 2500 3000 3500 4000 4500 5000 precision 1 2 3 4 5 6 7 8 time (in seconds) "Newton.dat" 2 3 4 5 6 7 8 9 10

• rder

1000 1500 2000 2500 3000 3500 4000 4500 5000 precision 5 10 15 20 25 30 time (in seconds)

Experimental results

N ... r 2 4 8 16 256 0.02 vs. 2.09 0.08 vs. 6 0.44 vs. 28 3 vs. 169 512 0.04 vs. 8 0.17 vs. 25 1 vs. 113 6.41 vs. 688 1024 0.08 vs. 32 0.39 vs. 104 2.30 vs. 484 15 vs. 2795 2048 0.18 vs. 128 0.94 vs. 424 5.54 vs. 2025 36 vs. > 3h ⋆ 4096 0.42 vs. 503 2.26 vs. 1686 13.69 vs. 8348 92 vs. > 1/2 day⋆ Basis, system with r equations, precision N: new vs. naive

Experimental results

2 4 6 8 10 12 14 16 500 1000 1500 2000 2500 3000 3500 4000 time (in seconds) precision ’DAC.dat’ 10 20 30 40 50 60 70 200 400 600 800 1000 1200 1400 1600 time (in seconds) precision ’dac.dat’ ’naive.dat’

Left: DAC computation of one solution for LDE of orders 2, 4, and 8 Right: DAC vs. naive, one solution of a LDE of order 2

The divide-and-conquer algorithm

Problem: solve Ly = 0, where L =

i ai(t)δi, with δ = t d dt.

Idea: the lowest terms of y(t) only depend on the lowest terms of ai. Proof: if y = y0 + tmy1, then L(δ)y = L(δ)y0 + tmL(δ + m)y1 DAC algorithm to solve L(δ)y = 0 mod t2m: determine y0, by recursively solving L(δ)y0 = 0 mod tm; compute the “error” R, such that L(δ)y0 = tmR mod t2m; compute y1, by recursively solving L(δ + m)y1 = −R mod tm. C(N) = 2 C(N/2) + O(rM(N)) ➥ C(N) = O(r M(N) log N) ◮ Newton method: computing a whole basis of solutions in 1. will allow to determine y1 in 3. from y0 and R alone, without a second recursive call: ˜ C(N) = ˜ C(N/2) + O(M(r, N)) ➥ ˜ C(N) = O(M(r, N))

The divide-and-conquer algorithm

Problem: solve Ly = 0, where L =

i ai(t)δi, with δ = t d dt.

Idea: the lowest terms of y(t) only depend on the lowest terms of ai. Proof: if y = y0 + tmy1, then L(δ)y = L(δ)y0 + tmL(δ + m)y1 DAC algorithm to solve L(δ)y = 0 mod t2m:

• 1. determine y0, by recursively solving L(δ)y0 = 0 mod tm;
• 2. compute the “error” R, such that L(δ)y0 = tmR mod t2m;
• 3. compute y1, by recursively solving L(δ + m)y1 = −R mod tm.

C(N) = 2 C(N/2) + O(rM(N)) ➥ C(N) = O(r M(N) log N) ◮ Newton method: computing a whole basis of solutions in 1. will allow to determine y1 in 3. from y0 and R alone, without a second recursive call: ˜ C(N) = ˜ C(N/2) + O(M(r, N)) ➥ ˜ C(N) = O(M(r, N))

The divide-and-conquer algorithm

Problem: solve Ly = 0, where L =

i ai(t)δi, with δ = t d dt.

Idea: the lowest terms of y(t) only depend on the lowest terms of ai. Proof: if y = y0 + tmy1, then L(δ)y = L(δ)y0 + tmL(δ + m)y1 DAC algorithm to solve L(δ)y = 0 mod t2m:

• 1. determine y0, by recursively solving L(δ)y0 = 0 mod tm;
• 2. compute the “error” R, such that L(δ)y0 = tmR mod t2m;
• 3. compute y1, by recursively solving L(δ + m)y1 = −R mod tm.

C(N) = 2 C(N/2) + O(rM(N)) ➥ C(N) = O(r M(N) log N) ◮ Newton method: computing a whole basis of solutions in 1. will allow to determine y1 in 3. from y0 and R alone, without a second recursive call: ˜ C(N) = ˜ C(N/2) + O(M(r, N)) ➥ ˜ C(N) = O(M(r, N))

Newton iteration: real case

x N(x)

xκ+1 = xκ − (x2

κ − 2)/(2xκ),

x0 = 1.5 x1 = 1.4166666666666666666666666666667 x2 = 1.4142156862745098039215686274510 x3 = 1.4142135623746899106262955788901 x4 = 1.4142135623730950488016896235025

Newton iteration: power series case

Let ϕ : K[[t]] → K[[t]]. To solve ϕ(g) = 0 in K[[t]], one can apply Newton’s tangent method: gκ+1 = gκ − ϕ(gκ) ϕ′(gκ) mod t2κ+1. ◮ The number of correct coefficients doubles after each iteration. ◮ Total cost = 2 ×

• the cost of the last iteration
• .

Theorem [Cook (1966), Sieveking (1972) & Kung (1974), Brent 1975] Division, logarithm and exponential of a power series in K[[t]] can be computed at precision N using O(M(N)) operations in K.

Division and logarithm of power series

◮ To compute the reciprocal of f ∈ K[[t]], choose ϕ(g) = 1/g − f: g0 = 1 f0 and gκ+1 = 2 gκ − fg2

κ

mod t2κ+1 for κ ≥ 0. ◮ division of power series at precision N in O(M(N)). ◮ log(f) = −

i≥1 1 i (1 − f)i of f ∈ 1 + tK[[t]] in O(M(N)) by:

• 1. computing the Taylor expansion of h = f ′/f modulo tN−1;
• 2. taking the antiderivative of h.

Exponentials of power series

◮ To compute the first N terms of the exponential exp(f) =

i≥0 1 i!f i,

choose ϕ(g) = log(g) − f. Iteration: g0 = 1 and gκ+1 = gκ−gκ (log(gκ) − f) mod t2κ+1 for κ ≥ 0. ◮ First order differential equations: compute the first N terms of f ∈ K[[t]] such that af ′ + bf = c.

• if c = 0 then the solution is

f0 = exp(−

• b/a)
• else, variation of constants: f = f0g, where g′ = c/af0.

Intermezzo: constant coefficients case

Problem: solve y′(t) = A y(t), y(0) = v, where A is a scalar matrix. Equivalent question: compute y(t) = exp(

• Av).

Warning: this equivalence is no longer true in the non-scalar case! Idea: the Laplace transform zN = N−1

i=0 Aivti of y(t) is the truncation

at order N of z =

i≥0 Aivti = (I − tA)−1v.

Algorithm (sketch):

• 1. compute z as a rational function of degre ≤ r (indep. on N);
• 2. deduce its Taylor expansion modulo tN: O(N/r M(r)) to expand

each coordinate of z. a b = c + d b = c +

• e + f

b

Brent & Kung’s algorithm for non-linear equations

• 1. reduce y′′(t) + a(t)y′(t) + b(t)y(t) = 0 to a 1st-order equation

◮ factor D2 + a(t)D + b(t) as (D + S(t))(D + T (t)): S + T = a, T ′ + ST = b, thus T ′ + aT − T 2 − b = 0

• 2. reduce the Riccati equation to a linear equation (by linearization)
• 3. find S, T and solve two linear 1st order equations to get y(t)

◮ generalizes to arbitrary orders: Lin(r, N) = NonLin(r − 1, N) + O(r M(N)) + Lin(r − 1, N) NonLin(r − 1, N) = O(Lin(r − 1, N)) ➥ Lin(r, N) = O(rr M(N)).

Newton iteration hits again

Suppose we have to solve a “functional” equation φ(Y ) = 0, where φ : Mr×r(K[[t]]) → Mr×r(K[[t]]) is differentiable. Define the sequence Yκ+1 = Yκ − Uκ+1, where

• Uκ+1 is a solution of valuation ≥ 2κ+1 of the linearized equation

Dφ|Yκ · U = φ(Yκ),

• Dφ|Yκ is the differential of φ at Yκ.

Then, the sequence Yκ converges quadratically to the solution Y .

First application: matrix inversion

To compute the inverse Z of a matrix Y of power series:

• choose the map φ : Z → I − Y Z with differential Z → −Y Z
• the equation for U becomes

−Y U = I − Y Zκ mod t2κ+1

• solution U = −Y −1 (I − Y Zκ) = −Zκ (I − Y Zκ) mod t2κ+1

This yields the Newton–Schulz iteration for Y −1 [Schulz, 1933] Zκ+1 = Zκ + Zκ(Ir − Y Zκ) mod t2κ+1. Cinv(N) = Cinv(N/2) + O(M(r, N)) ➥ Cinv(N) = O(M(r, N))

Second application: solving differential equations

To compute the solution Y of the system Y ′ = AY

• choose the map φ : Y → Y ′ − AY , with differential φ.
• the equation for U is

U ′ − AU = Y ′

κ − AYκ

mod t2κ+1

• using Lagrange’s method of variation of parameters, solution

U = YκVκ mod t2κ+1, Y ′

κ − AYκ = YκV ′ κ mod t2κ+1.

This yields the BCOSSS iteration for Y : Yκ+1 = Yκ − Yκ

• Y −1

κ

• Y ′

κ − AYκ

• mod t2κ+1.

Csolve(N) = Csolve(N/2) + O(M(r, N)) ➥ Csolve(N) = O(M(r, N))

Further questions

constant factor improvements: middle products of polynomial matrices small characteristic case: Pad´ e approximants? p-adic lifting? faster Newton for the case of a single equation: exploit companion form bit complexity analysis