Tensor numerical approach to space-time approximation of - - PDF document

SLIDE 1

Tensor numerical approach to space-time approximation

• f multi-dimensional parabolic equations

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016

Max-Planck-Institute for Mathematics in the Sciences, Leipzig, Germany

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 1 / 38

Tensor methods vs. the “curse of dimensionality”: Nd → dN → d log N ?

PDE based applications in higher dimensions: ◮ Elliptic (parameter-dependent) BV/spectral problems: u ∈ H1

0 (Ω) :

Hu := − div (a grad u + uv) + Vu = F (= λu) in Ω ∈ Rd. ◮ Parabolic-type dynamics: Find u : Rd × (0, T) → R, s.t. u(x, 0) ∈ H2(Rd) : σ ∂u ∂t + Hu = F, σ ∈ {1, i}. Computational challenges: Multi-dimensionality, multi-parametric problems, huge grids in Rd and in t ∈ [0, T], multi-dimensional integral transforms (convolution), high oscillations in coefficients, redundancy in data representation. Big Data ≇ Big Knowledge ! (curse of redundancy) ⇒ Knowledge = o(Data) ? Parallelization issues

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 2 / 38

SLIDE 2

Separability concept in comput. quantum chemistry: spectra/dynamics 1929, P.A.M. Dirac: The fundamental laws necessary for the mathematical treatment of large part of physics and the whole of chemistry are thus completely known, and the difficulty lies only in the fact that application of these laws leads to equations that are too complex to be solved. 1961, R. Bellman: In view of all that we have said in the foregoing sections, the many obstacles we appear to have surmounted, what casts the pall over our victory celebration? It is the curse of dimensionality, a malediction that has plagued the scientist from earliest days. ◮ Low-rank nonlinear tensor approximation vs. information redundancy. 1927-1966, Hitchcock, Tucker: Canonical (CP) and Tucker tensor formats in chemometrics. 1992, S. White: DMRG and Matrix Product States for quantum spin systems. 2003, Dirac-Frenkel molecular dynamics on rank-1 tensor manifold. Since 2006: Tensor numerical methods for PDEs in Rd.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 3 / 38

Outline of the talk

1

Big-data compression via separation of variables: basic tensor formats

2

Range-separated tensor format: reduced model for many-particle interactions

3

Quantized-TT tensor approximation (QTT) of logarithmic complexity

4

O(log n) numerical quadratures by simple QTT algebraic operations

5

Elliptic PDEs with highly oscillating coefficients in log-complexity: tensor approximation vs. asymptotic homogenization

6

PDEs on complex geometry: Low-rank approximation in iso-geometric analysis

7

d-dimensional dynamics: time stepping vs. equations in 2 + 1, 3 + 1,...,d + 1

QTT integration for retarded potential IE (wave equation) Rank estimates Heat eqn., Fokker-Planck, chemical master equations

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 4 / 38

SLIDE 3

Separable representation in tensor-product Hilbert space

◮ Euclidean vector space Vn = Rn1×...×nd = d

ℓ=1 Rnℓ, n = (n1, ..., nd),

V = [vi] ∈ Vn : W, V =

• iwivi,

i = (i1, ..., id) : iℓ ∈ {1, ..., nℓ}. Separable representations in Vn via rank-1 tensors: V = [vi1...id ] = v(1) ⊗ . . . ⊗ v(d) ∈ Vn : vi1...id = d

ℓ=1v (ℓ) iℓ .

◮ Storage: Stor(V) = d

ℓ=1 nℓ ≪ dim Vn = d ℓ=1 nℓ

(nd ≪ nd). ◮ The scalar product reduces to univariate operations W, V = w(1) ⊗ . . . ⊗ w(d), v(1) ⊗ . . . ⊗ v(d) = d

ℓ=1w(ℓ), v(ℓ)Vℓ.

Fast multiliear algebra: d-dimensional Hadamard product, contraction, convolution product, addition etc. all reduce to 1D-operations!

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 5 / 38

Low-parametric separable representations: canonical, Tucker, MPS (TT) formats ◮ Most commonly used tensor formats extend rank-R matrices: d = 2 : V = R

k=1ukvT k ≡ G1G T 2 ≡ UDV T ∈ Rm×n.

– Canonical (CP) tensors (multidimensional rank-R representation) [Hitchcock ’27]. – The orthogonal Tucker decomposition (multidimensional truncated SVD) [Tucker ’66] – Matrix product states (MPS) factorization [S. White ’92 et al.] ⇒ Variants of MPS: Tensor train (TT) [Oseledets, Tyrtyshnikov ’09]; Hierarchical Tucker (HT) [Hackbusch, Kühn ’09].

• Def. Canonical R-term representation in Vn: V ∈ CR(Vn),

V = R

k=1v(1) k

⊗ . . . ⊗ v(d)

k ,

v(ℓ)

k

∈ Rnℓ. ◮ Advantages: Storage = dRN, simple multilinear algebra, but hard for approximation.

+ b

A

1

b V V V V V V V V V + = ... +

1 1 2 2 2

r r r

(1) (1) (1) (2) (2) (2) 2 1 (3) (3) (3)

r

b

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 6 / 38

SLIDE 4

Orthogonal Tucker decomposition

• Def. Rank r = (r1, . . . , rd) Tucker tensors: V ∈ d

ℓ=1 Tℓ ⊂ Vn , V (ℓ) ∈ Rnℓ×rℓ

V =

r

• k1,...,kd =1

bk1...kd v(1)

k1 ⊗. . .⊗v(d) kd = B ×1 V (1) ×2 · · · ×d V (d),

Tℓ = span{v(ℓ)

kℓ }rℓ kℓ=1 ⊂ Rnℓ

◮ Storage: drN + r d, r = max rℓ ≪ N = max nℓ. For functional tensors r = O(log N).

B r3 I3 I2 r2 (3) (2) I1 I I I 2 3 1 r2 r3 r1 r1 (1)

V T T T

5 10 5 10 0.2 0.4 0.6 0.8 1 2 4 6 8 10 12 14 10

−10

10

−5

10

Tucker rank error Slater function, AR=10, n = 64

EFN EFE EC

Basics of tensor numerical approximation [BNK, CMAM ’06], [BNK, Khoromskaia ’07]: Low-rank decomposition of functional tensors, e.g. for v(x) = e−x,

1 x, e−x x :

V − Vr ≤ e−cr.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 7 / 38

Matrix Product States (Tensor Train) factorization

• Def. MPS (TT) format: Given r = (r1, ..., rd), rd = 1. V = [vi] ∈ TT[r] ⊂ Vn is

parametrized by contracted product of tri-tensors in Rrℓ−1×nℓ×rℓ, vi1...id =

• αG (1)

α1 [i1]G (2) α1α2[i2] · · · G (d) αd−1[id] ≡ G (1)[i1]G (2)[i2]...G (d)[id],

G (ℓ)[iℓ] is an rℓ−1 × rℓ matrix, 1 ≤ iℓ ≤ nℓ. Storage: dr 2N ≪ Nd.

.

r r r r r r r 1 1 2 3 3 4 4 i i i i 1 2 4 5 n n n n n 1 2 3 4 5 2 r 3 i

a

i2

=

1 i i3i4 i5

d = 5: A = [ai1i2i3i4i5] ∈ Rn1×n2×n3×n4×n5. For r = 1 all formats coincide ! 1D (univariate) function related vectors f = {f (xi)} ∈ RN are non-compressible ? ◮ Find a hidden low-rank tensor structure in large functional vectors and matrices !

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 8 / 38

SLIDE 5

Range-separated (RS) tensor format for approximating funct. with multiple cusps

Figure: Villin protein (left). Short- and long-range parts of the reference Newton kernel 1/x.

[Benner, Khoromskaia, BNK ’16]

◮ Long range interaction potential in the large N0-particle system, xk ∈ R3, k = 1, ..., N0, and the radial function p(x) (say, p(x) = 1/x) P(x) = N0

ν=1zk p(x − xk),

zk ∈ R, xk, x ∈ Ω = [−b, b]3, ◮ The interaction energy of the system of N0 charged particles EN = EN(x1, . . . , xN) = 1 2

N

• j=1

zj

N

• k=1,k=j

zk xj − xk. (1)

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 9 / 38

RS canonical/Tucker/TT formats: beneficial properties

[Benner, Khoromskaia, BNK, arXiv: 1606.09218, 2016]

Definition

(RS-canonical (Tucker, TT) tensors). The RS-canonical tensor format defines the class

• f d-tensors A ∈ Rn1×···×nd , represented as a sum of a rank-R CP tensor U and a

cumulated CP tensor generated by localized U0, s.t. rank(Uν) = rank(U0) ≤ R0, Uν = Replica(U0), A = R

k=1ξku(1) k

⊗ · · · ⊗ u(d)

k

+ N0

ν=1cνUν,

with diam(suppUν) ≤ n0 = O(1).

Theorem

The storage size for RS-canonical tensor is estimated by Stor(A) ≤ dRn + (d + 1)N0 + dR0n0. Each entry of an RS-CP tensor can be calculated at O(dR + dR0n0) cost. For ε-rank of the Tucker approximation to the long-range CP tensor U, r0 = rankTuck(U): |r0| ≤ C b log3/2(| log(ε/N0)|), rankCan(U) ≤ |r0|2.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 10 / 38

SLIDE 6

Tensor approximation of the quantized vectors in log-complexity !

The ”curse of dimensionality“ ⇒ "blessing dimensions" !

F

N=23 L=log N=3 + ...

◮ [BNK ’2009] Quantized Tensor Approximation (QTT) of functional N-vectors (N = 2L). Isometry QL : [xi]N

i=1 =: x → A := [aj] ∈ QL :=

L

ℓ=1R2,

aj := xi, i → j ∈ {1, 2}⊗L : i − 1 = L

ℓ=1(jℓ − 1)2ℓ−1

(binary coding). TT approximation of quantized image in QL ⇒ QTT method ◮ Storage in QTT tensor format scales logarithmically in N = 2L, 2r 2L ≪ 2L.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 11 / 38

QTT-approxim. theory: Why the super-compressed QTT-approximation does a job?

• Thm. [BNK ’09-’10]. QTT-approximation of functional N-vectors in O(log N)-op. N = 2L.

◮ For quantized exponential N-vector: rankQTT(x) = rankCan(QL(x)) = 1, x := {zn−1}N

n=1 ∈ CN → ⊗L p=1

• 1

z2p−1

• ∈

L

• p=1

C2, z ∈ C. ◮ For quantized sin and cos N-vector: rankQCan,C(x) = rankQTT,R(x) = 2, x := {sin(ωh(n − 1))}N

n=1 ∈ CN,

h = 1/(N − 1), ∀ω ∈ C. ◮ QTT of a polynomial of degree m, rankQTT(Pm) ≤ m + 1. ⇒ P.w. polynomials. ◮ Smoothly modulated multiply replicated functions. ◮ Gaussian e−λx2 has small ε-rank in QTT and QCan formats: r = O(| log ε|). ◮ Approx. I(ω, f ) =

• Ω f (x)eiωg(x)dx as a function of frequency ω [BNK, Sauter, Veit 12’-’15]

◮ Contribution to QTT approx. of functions [BNK, Oseledets, Dolgov, Grasedyck, Kazeev, Schwab] ◮ Numerical observation: 2L × 2L 1D-Laplacian reshapes to a low TT-rank [Oseledets ’09]

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 12 / 38

SLIDE 7

Tensor representation of operators (large multidimensional matrices)

Example: d-dimensional discrete Laplacian as the Kronecker product sum ∆d = ∆1 ⊗ I ⊗ ... ⊗ I + I ⊗ ∆1 ⊗ I... ⊗ I + ... + I ⊗ I... ⊗ ∆1 ∈ RN⊗d ×N⊗d , ∆1 = tridiag{−1, 2, −1} ∈ RN×N, I is the N × N identity. ◮ Canonical/Tucker representation: rankCP(∆d) = d, rankTuck(∆d) = 2. ◮ TT/QTT representation: rankTT(∆d) = 2, rankQTT(∆d) ≤ 4, ∀d. ∆d =

• ∆1

I

• ⋊

⋉

• I

∆1 I ⋊

⋉(d−2)

⋊ ⋉ I ∆1

• .

◮ [Kazeev, BNK ’10] Explicit QTT representation: rankQTT(∆1) = 3, rankQTT(∆−1

1 ) ≤ 5.

∆1 =

• I

J′ J

• ⋊

⋉   I J′ J J J′  

⋊ ⋉(d−2)

⋊ ⋉   2I − J − J′ −J −J′   .

I =

• 1

1

• ,

J =

• 1
• .

“⋊ ⋉” is a regular matrix product of block core matrices, blocks being multiplied by means of tensor product. Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 13 / 38

Beyond Galerkin: Nonlinear approximation in low-rank tensor formats

Approximation problem: Given A ∈ Vn (in general, A ∈ S0 ⊂ Vn), find Tr(A) := argmin

X∈S

X − A, where S ⊂ {Tr, CR, MPS/TT[r]]}. Quasi-optimal (nonlinear) tensor approx. is based on: SVD, QR, Cholesky + ALS iter. ◮ SVD-based HOSVD for Tucker tensors [De Lathauwer et al. 2000] ◮ SVD-based HOSVD for MPS/TT tensors [Vidal ’06] Tucker and MPS/TT ranks: Tr := {A ∈ Vn : rank A(p) ≤ rp} : A(p) = [a(j1 j2 . . . jp−1

• row index

; jp

• column index

; jp+1 . . . jd

• row index

)]. TT[r] := {A ∈ Vn : rank A[p] ≤ rp} : A[p] = [a(j1 j2 . . . jp

• column index

; jp+1 . . . jd

• row index

)]. Canonical rank can not be presented as a matrix rank ⇒ unstable approximation :( ◮ Rank reduction in the canonical format. Reduced HOSVD: Canonical → Tucker → canonical ALS for the Tucker core [BNK, Khoromskaia, SISC ’08].

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 14 / 38

SLIDE 8

• I. QTT based quadratures of O(log n) = O(| log ε|) complexity

◮ Quantized integrand f (x) and weight function w(x) with moderate QTT-ranks. The rectangular n-point quadrature, n = 2L, |I − In| = O(n−α), Cost = O(r 2 log n). 1

−1

w(x)f (x)dx ≈ In(f ) := h

n

• i=1

w(xi)f (xi) = W, FQTT, W, F ∈ ⊗L

ℓ=1R2.

1 2 3 4 5 6 7 8 9 10 −2 2 4 6 8 10 12 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −2 −1.5 −1 −0.5 0.5 1 1.5 2

◮ Retarded potentials: [BNK, Sauter, Veit, CMAM 2011]. ◮ Twofold-QTT: fast quadratures for highly oscillating integrals + QTT in ω. Compute funct. of ω ∈ [ωmin, ωmax] for various f (x), g(x), [BNK, Veit, CMAM 2015]. I(ω, f ) :=

• Ω

f (x)ei ωg(x)dx.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 15 / 38

• II. Asymptotic homogenization for elliptic PDEs with oscillating coefficients

[BNK, Repin: RJNAMM ’15]

Traditional FEM-Galerkin p.w.l. approximation in Ω = (0, 1)d: Find uǫ ∈ H1

0(Ω) :

− div (aǫ∇uǫ) = f in Ω, f ∈ L2 (Ω) , (2) with oscillating or lattice-structured/replicated coefficients in Ω = ∪iΠǫ

i ,

aǫ(x) := A(x − xi ǫ ) for x ∈ Πǫ

i – scaled unit cell.

◮ Homogenization methods suggest asymptotically uǫ − uhomo ≤ C√ǫ as ǫ → 0. Challanges: "Defected" periodic coeff.; Moderate ǫ; Huge grids: nd ≈ (n0/ǫ)−d. Assumptions: (A) aǫ(x) and “homogenized coeff.” a0(x) = 1

2(a+(x) + a−(x)), have low QTT rank.

a+(x) and a−(x) are the upper and lower majorants of aǫ(x). (B) Solution of the “homogenized” equation A0u0 = f is cheap. a0(x) → A0.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 16 / 38

SLIDE 9

QTT approximation of PDEs with “defected” periodic coefficients: modulated functions ◮ computational scheme: Solve the FEM-Galerkin approximation of (2) Aǫuǫ = f, Aǫ ∈ RNd ×Nd , uǫ, f ∈ RNd , N = 2L (3) in quantized tensor space QL by the rank-truncated preconditioned iteration: u0 = A−1

0 f :

(E + βA−1

0 Aǫ − E)u ≡ (E + B)u = βu0,

β > 0, assuming that B = q < 1 for B = βA−1

0 Aǫ − E, and rankQTT(B) is small.

Given tol. δ > 0, u0 ∈ Q(r)

L

: uk+1 = TQTT,δ (βu0 − Buk) → uδ, k = 0, 1, 2, ... uk+1 − uk ≤ qku0 up to the truncation level δ > 0.

500 1000 1500 2000 0.5 1 1.5 2 500 1000 1500 2000 0.5 1 1.5 2

1 2 3 4 5 6 x 10

4

0.2 0.4 0.6 0.8 1 1.2 1.4

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 17 / 38

QTT rank-truncated iteration of log-cost O(log N) vs. homogenization method

grid size 2L 213, (it) 214, (it) 215, (it) 216, (it) 217, (it) r(aǫ) r(uδ) C + sin(ωx) 0.97, (5) 1.2, (5) 1.3, (5) 2.0, (6) 2.1, (6) 2.67 3.7 4-steps coef. 3.4, (9) 4.3, (9) 4.5, (9) 6.7, (9) 14.3, (14) 2.9 4.96 C + sin(ωx3) 5.3, (5) 10.0, (6) 9.95, (6) 11.98, (6) 16.2, (5) 7.53 8.24 CPU sec., iter. history, QTT-ranks, ω = 1/ǫ = 2π64, δQTT = 10−7: uǫ − uδ0 ⋍ 10−7, uǫ − uδ1 ⋍ 10−6.

1 1.5 2 2.5 3 3.5 4 4.5 5 10

−8

10

−7

10

−6

10

−5

10

−4

10

−3

1 2 3 4 5 6 7 8 9 10

−8

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 10

−8

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

• Precond. Steepest Descent in QTT-truncated iterations for the periodic, step-, and cubically-oscillating coef.

500 1000 1500 2000 0.5 1 1.5 500 1000 1500 2000 0.5 1 1.5 500 1000 1500 2000 0.5 1 1.5

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 18 / 38

SLIDE 10

From asymptotic to stochastic homogenization and beyond

2D diffusion coefficients aǫ(x), rhs, and the respective solutions: 400x400-grid in 5 sec. f1(x, y) = sin(2x) sin(2y); f2(x, y) = sin(2x).

1 0.8 0.6 0.4 0.2 0.2 0.4 0.6 0.8 1 0.5

• 0.5
• 1

1 1 0.8 0.6 0.4 0.2 0.2 0.4 0.6 0.8 4 2

• 2
• 4

1 1 0.5 1 0.8 0.6 0.4 0.2 0.5

• 0.5
• 1

1 1 0.8 0.6 0.4 0.2 1 0.8 0.6 0.4 0.2

• 10

10 5

• 5

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 19 / 38

• III. Tensor approach to Iso-Geometric Analysis (IGA): toward dynamics

[Mantzaflaris, Jüttler, BNK, Langer: ’15 - ’16] – Low-rank IGA in 2D and 3D: elliptic PDEs on complicated geometry.

Figure: The yeti footprint domain, parameterized by 21 quadratic B-spline patches. Patch segmentation (left),

control grids (middle) and the (absolute) Jacobian determinant are shown. On the left picture rank(det ∇G) and max(rank(qk,ℓ), k, ℓ ∈ {1, 2}) are given. ǫ = 10−10. Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 20 / 38

SLIDE 11

IGA: Matrix generation via patch-wise low-rank representation

FEM-Galerkin p.w.l. approximation in Ω ∈ Rd: Find u ∈ H1

0(Ω) :

− div (a(x)∇u) = f in Ω, f ∈ L2 (Ω) .

Figure: The Jacobian determinant values of the selected patch (Fig. 2) on the domain (left). The graph and

control net in parameter domain (right). Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 21 / 38

• IV. Temporal problems: Time stepping vs. global x − t schemes in d + 1

Problem 1. Complex-time molecular Schrödinger eq. in QMD, d = 3N. i ∂ψ ∂t = Hψ = (−1 2∆d + V )ψ, ψ(x, 0) = ψ0(x), x ∈ Rd, V : Rd → R is (given) approximation to the potential energy surface (PES). Problem 2. Real-time evolution. The Fokker-Planck equation ψ(0) = ψ0, dψ dt = −Aψ; Aψ = −ε∆ψ + div(ψv), ψ : Rd → R, v : Rd → Rd is a given velocity field. ψ(t) → ψ∗ : Aψ∗ = 0 (small eigs. of A)

Example 1’. Heat equation in Rd × [0, T]:

dψ dt

= −∆ψ + f , ψ(0) = ψ0.

Problem 3. Chemical master equation for joint probability density P(x, t), P(x, 0) = P0, dP(x, t) dt = AP(x, t), x ∈ Rn1×...×nd . Problem 4. Wave equation in unbounded exterior domain (Dirichlet problem) ∂2

t ψ − ∆ψ = 0

in Ω × [0, T], ψ(·, 0) = ∂tψ(·, 0) = 0 in Ω ∈ R3.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 22 / 38

SLIDE 12

Main approaches (all methods approximate e−tHψ0!)

◮ Time integrators Sparse grids in (x, t): [Schwab et al.; Griebel et al.] Dirac-Frenkel projection onto Tucker/TT/QTT-manifold S, dy dt − Hy, δy = 0, δy ∈ TyS.

[Meyer et al. ’03; Lubich et al. ’07-’15]

Wave equation in exterior domain: Tensor quadrature for retarded potentials

[BNK, Sauter, Veit ’11]

Time stepping by implicit scheme in TT/QTT + ALS/DMRG local solver,

[Dolgov, BNK, Oseledets ’11; Schwab et al. ’12]

◮ Global time-space schemes in d + 1 Time-space separation by QTT-Cayley transform, rank bounds [Gavrilyuk, BNK ’11] QTT-Tucker + ALS-type solver on global (x, t) tensor manifold

[Dolgov, BNK ’12-’14] Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 23 / 38

Wave equation: Retarded potential integral equation via QTT integration

[BNK, Sauter, Veit: CMAM ’11] + work in progress ...

Ω ⊂ R3 – a Lipschitz domain, Γ = ∂Ω. Homogeneous wave equation ∂2

t u − ∆u = 0

in Ω × [0, T] ; u = g

• n Γ × [0, T] .

In applications, Ω is often the unbounded exterior of a bounded domain. We employ an ansatz as a single layer potential for the solution u u(x, t) := Sϕ(x, t) :=

• Γ

ϕ(y, t − x − y) 4πx − y dΓy, (x, t) ∈ Ω × [0, T] . S is referred to as retarded single layer potential (the retarded time argument t − x − y which connects time and space variables). The boundary integral equation for unknown density ϕ,

• Γ

ϕ(y, t − x − y) 4πx − y dΓy = g(x, t) ∀(x, t) ∈ Γ × [0, T] . (4)

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 24 / 38

SLIDE 13

QTT-Cayley transform method

◮ Dynamics and spectrum of high-dim. Hamiltonians [Gavrilyuk, BNK ’11] ψ(t) =

∞

• p=0

L(0)

p (t)up ≡ i(H + iI)−1 ∞

• p=0

L(0)

p (t)G pψ0,

t ∈ [0, T], G = H(H + iI)−1, Lp - Laguere polynomials. Recursion for up, u0 = i(H + iI)−1ψ0, up+1 = Gup, p = 0, 1, ... ◮ The m-term truncated series representation ψm(t) =

m

• p=0

L(0)

p (t)up.

(5) Def. f = akϕk(x) is called H-analytic if there is C = C(f ) > 0, s.t. Hnf =

• ∞
• k=0

a2

kλ2n k ≤ C nn!

for n = 1, 2, 3, ... where {ϕk, λk} is the eigenpear of H = H∗.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 25 / 38

Rank bounds

• Rem. For H-analytic vector f the power series

∞

• n=0

sn n! Hnf =: f s,H possesses a

positive convergence radius r > 0, i.e., f s,H < ∞ if 0 ≤ s < r.

• Lem. Let ψ0 be H-analytical, then for every fixed s > 0, and fixed T > 0, s.t.

ψ0s,H < ∞, the ψm(t) converges exponentially in m, ψ(t) − ψm(t) ≤ cm−1/12e−c1 3

√mψ0s,H,

t ∈ [0, T] where c, c1 > 0 are independent of m. Lem. QTT-rank of a tensor U = [ψm(t1), ..., ψm(tNt )]Nt

k=1, tk = k∆t, is bounded by

rankQTT(U) ≤

m

• p=0

(p + 1)rankQTTG p(ψ0). For the harmonic oscillator: rankQTT(U) ≤ Cm2rankQTT(ψ0).

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 26 / 38

SLIDE 14

Spectral calculus, complexity, computational difficulties

◮ The energy spectrum is recovered by FFT of autocorrelation function a(t) = ψm(t), ψ(0), S(E) = ∞ a(t)eiEtdt. ◮ Simultaneous time-space QTT decomposition of complexity O(log Nt log N m4 (rankQTT(G mψ0))2), m = O(logp ε). ◮ Difficulties on examples of the molecular Schrödinger, Fokker-Planck and master equations: Even the low-rank Hamiltonian does not guarantee the low rank of the evolving solution ! The Hamiltonian operator might not be well-separable (convection-diffusion). Boundary conditions (layer) may enforce rank instabilities (approximate b.c.). Large mode size (quantization).

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 27 / 38

High-dim. spectral calculations for the Henon-Heiles Hamiltonians (QMD)

[BNK, Oseledets ’11] Perturbed harmonic oscillator

H = −1 2∆ + 1 2

d

• k=1

q2

k + λ d−1

• k=1
• q2

kqk+1 − 1

3q3

k

• ,

rankTT(H) ≤ 3, rankQTT(H) ≤ 4. ◮ The energy spectrum is recovered by FFT of autocorrelation function a(t) = ψm(t), ψ(0), S(E) = ∞ a(t)eiEtdt.

Spectrum of the Henon-Heiles Hamiltonian (d = 4, 10), n = 27, QTT-rank ≤ 10. Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 28 / 38

SLIDE 15

Discretization on global QTT-tensor manifold in space-time

Given A, fk in the TT/QTT format, apply Crank-Nicolson, (I + τ 2 A)yk+1 = (I − τ 2 A)yk + τ 2 (fk + fk+1) =: Fk+1. (A) Time stepping by DMRG-TT iter. for solving (I + τ 2 A)yk+1 = Fk+1. (B) Global O(log Nx log Nt) block-AMEn solver in QTT-Tucker format: yk+1 − yk + τ 2 Ayk+1 + τ 2 Ayk = τ 2 (fk + fk+1). Solve the huge global Nd

x × Nt system in QTT format by AMEn iteration

     I + τ

2 A

−I + τ

2 A

I + τ

2 A

... ... −I + τ

2 A

I + τ

2 A

          y1 y2 . . . yN      =     

• I − τ

2 A

• y 0

. . .      + τ 2      f0 + f1 f1 + f2 . . . fN−1 + fN      .

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 29 / 38

Heat equation: global (d + 1)-scheme vs. time stepping

As a “sanity” test we consider the heat equation u(0) = g :

∂u ∂t − ∆u = f

in Ω = [0, 1]2, u|∂Ω = 0. (A) The case of a smooth analytic solution (refinement of the grids). Choose f = 0, g(x1, x2) = sin(πx1) sin(πx2). The analytical solution at the time t: u∗(x1, x2) = g(x1, x2) exp(−2π2t). The time interval is [0, 1/2], the residual tolerance for the TT-solve algorithm: 10−6.

Table:

||u−u∗||F ||u∗||F

versus the spatial Nx and time Nt grid sizes at t = 1/2. Nt \ Nx 28 29 210 211 27 4.7598e-03 4.8515e-03 4.8746e-03 4.8803e-03 29 1.8271e-04 2.7475e-04 2.9786e-04 3.0363e-04 211 1.0380e-04 1.1745e-05 1.1363e-05 1.7152e-05 213 1.2171e-04 2.9652e-05 6.5401e-06 7.8656e-07

The convergence is faster than the theoretical bound O(N−2

t

+ N−2

x

) of the Crank-Nicolson - FD scheme. The solution time of the block scheme is almost independent on Nt, Nx (100 − 200 milliseconds for any test in Table) since the QTT ranks are uniformly bounded by a small constant.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 30 / 38

SLIDE 16

Heat equation: non-regular rhs

(B). f (x1, x2) = g = 1, x1, x2 ∈ (0, 1), i.e. discontinuous at the boundary. The time step should be small to resolve transitional processes near Γ, otherwise, as the second-order scheme is not monotonous, the oscillations occur, leading to large tensor ranks. t ∈ [0, 1], Nx = Ny = 256, and the QTT tolerance is 10−6. We compare the time stepping solution scheme with the block approach. The number of time steps varies 28 to 216, the CPU times/sec.

Table: CPU times (sec.), relative residuals, the average QTT rank in the block scheme.

Block solution Time stepping Nt CPU time

||−∆u(1)−f || ||f ||

rank CPU time

||−∆u(1)−f || ||f ||

28 76.96 1.53e+03 37.17 1310.4 1.53e+03 210 83.23 7.34e-03 41.09 4547.1 7.38e-03 212 69.32 1.44e-04 39.40 3845.9 1.69e-03 214 60.40 3.67e-05 35.91 6232.4 1.58e-04 216 61.37 5.49e-05 32.98 12707 7.24e-05

Notice: the solution time even decreases with the number of time steps in the block algorithm since the smaller the time step, the smoother the solution is. The QTT-rank is stable with respect to the number of time steps, and manifests also a slight decrease due to the improving smoothness of the solution. It demonstrates the advantages of the logarithmic scaling of the block scheme, especially if we have to choose extremely small time steps.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 31 / 38

The Fokker-Planck equation

Real-time evolution. The Fokker-Planck equation dψ dt = −Aψ, ψ(0) = ψ0; Aψ = −ε∆ψ + div(ψv), the probability density ψ : Rd → R, and v : Rd → Rd is a given velocity field. Time propagation scheme with tensor truncation ψk+1 = Tε(Sψk), ψk ≈ ψ(tk), tk = τk. The time-propagator S: S = (I + τA)−1 (implicit Euler), S = (I + τ

2 A)−1(I − τ 2 A) (Crank-Nicolson).

The time evolution converges to the null-space of A, ψk → ψ∗, Aψ∗ = 0,

• ψ∗dx1 . . . dxd = 1.

ψ∗ is normalized probability density.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 32 / 38

SLIDE 17

Model problem for the Fokker-Planck eqn.

The dumbbell model discretized on large grids. v = Kq + grad(ϕ), K = β   1   , q =   q1 q2 q3   The potential energy ϕ is given as ϕ = 1 2(q2

1 + q2 2 + q2 3) + 1

2 z p3 e−(q2

1+q2 2+q2 3)/(2p2).

The following functional of the solution is of interest: τ(t) =

• ψ(t) (q ⊗ grad(ϕ)) dq.

In particular, we test η(t) = −τ12 β , Ψ(t) = −τ11 − τ22 β .

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 33 / 38

Numerics for the Fokker-Planck equation: time stepping vs. block solver

The dumbbell model: Complexity vs. Nt and Nx.

Figure: TT-solution time vs. τ,

Nx = 2dx , Nt = 2dt , Time ∼ O(r 2NtNx).

Figure: Block QTT-solution time

• vs. dt, dx, O(r2 log Nt log Nx ),

r = O(log Nx) = O(dx).

[Dolgov, BNK, Oseledets ’SISC-11] Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 34 / 38

SLIDE 18

Chemical master equation (CME): High-dimensional real-time dynamics

◮ QTT-Tucker for chemical master equation: Species S1, ..., Sd react in M channels. The vector of concentration x = (x1, ..., xd), xi ∈ {0, ..., Ni − 1}, i = 1, ..., d. The propensity function w m(x), m = 1, ..., M. The stoichiometric vector zm ∈ Zd.

Jz =              · · · 1 ... ... ... 1 ... . . .              ← row N − z ← row N , if z ≥ 0; Jz = (J−z )⊤, if z < 0.

CME is a deterministic difference equation on the joint probability density P(x, t): dP(t) dt =

M

• m=1

(Jzm − J0) diag(w m)P(t), P(t) ∈ R

d

i=1 Ni

+

, Jz = Jz1 ⊗ · · · ⊗ Jzd , w m = {w m(x)} and P(t) = {P(x, t)}, x ∈

d

• i=1

{0, . . . , Ni − 1}.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 35 / 38

Long-time dynamics: CME for the signaling cascade, d = 20

[Dolgov, BNK, NLAA ’13]

Figure: Cascade signaling network

S1 S2 · · · Sd d = 20, M = 40; for m = 1: w m(x) = 0.7, zm = −δm: generation of the first protein; for m = 2, ..., 20: w m(x) = xm−1 5 + xm−1 , zm = −δm: succeeding creation reactions; for m = 21, ..., 40: w m(x) = 0.07 · xm−20, zm = δm−20: destruction reactions. Ni = 63. δm is the m-th identity vector. Problem size N = 6420.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 36 / 38

SLIDE 19

Convergence history, O(log Nt) CPU time

Two-level t-stepping: numbers of time steps Nt in each interval [(p − 1)T0, pT0], the time interval widths T0 (coarse step): NT = T/T0. logarithmic complexity in Nt

Figure: Mean concentrations

xi (t).

Figure: Closeness to the kernel

||AP|| ||P|| (t)

Figure: CPU time (sec.)

• vs. log2(Nt), T0 = 15.

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 37 / 38

Summary: Rank-structured tensor methods for PDEs in Rd and beyond

Tensor MLA, O(d log N)-approximation in QTT tensor format. New RS tensor formats: approximation tool for complicated data (+). Super-fast computation of oscillatory integrals in O(d log N) cost (±). QTT approximation of elliptic PDEs with “defected“ periodic coefficients (±). Low-rank IGA analysis for d = 2, 3 (±). Multidimensional dynamics by QTT-AMEn solver for global (x, t)-system (∓). Time dependent Fokker-Planck and master eqn. in d + 1 (±). High-dimensional horizons: log-scaling in time and space discretizations This talk is based on the joint works:

[Benner, Khoromskaia, BNK ’16], arXiv: 1606.09218, 2016. [Dolgov, BNK ’13], [Dolgov, BNK, Oseledets ’12] [Mantzaflaris, Jüttler, BNK, Langer ’15-’16], LNCS 9213, 2015; MPI preprint, 2016. [BNK, Repin ’15], RJNAMM 2015. [BNK, Sauter, Veit], CMAM, ’11-’15.

Thank you for your attention !

Boris N. Khoromskij Linz, RICAM, Research Semester, WS2, 10.11.2016 Tensor approximation of parabolic PDEs 38 / 38