V. Simoncini - Dipartimento di Matematica

On the decay of the inverse of matrices
that are sum of Kronecker products
V. Simoncini
Dipartimento di Matematica, Universit`a di Bologna
[email protected]
Joint work with
C. Canuto (Politecnico Torino) and M. Verani (Politecnico Milano)
1
The application. I
Adaptive Legendre-Galerkin discretizations for PDEs:
H01 Tensorized Babuska-Shen basis in Ω = (0, 1) × (0, 1):
ηk (x1 , x2 ) = ηk1 (x1 )ηk2 (x2 ),
k1 , k2 ≥ 2,
k = (k1 , k2 )
{ηki }: ki -order Legendre polyn (1D BS basis)
Stiffness matrix:
(ηk , ηm )H 1 (Ω) = (ηk1 , ηm1 )H 1 (I) (ηk2 , ηm2 )L2 (I) +(ηk1 , ηm1 )L2 (I) (ηk2 , ηm2 )H 1 (I)
0
0
0
Kronecker structure: Sηp = Mp ⊗ Ip + Ip ⊗ Mp (max p polyn degree)
Note: If higher order polynomial used, then Sηp simply expands
(augmented Mp )
2
The application. I
Adaptive Legendre-Galerkin discretizations for PDEs:
H01 Tensorized Babuska-Shen basis in Ω = (0, 1) × (0, 1):
ηk (x1 , x2 ) = ηk1 (x1 )ηk2 (x2 ),
k1 , k2 ≥ 2,
k = (k1 , k2 )
{ηki }: ki -order Legendre polyn (1D BS basis)
Stiffness matrix:
(ηk , ηm )H 1 (Ω) = (ηk1 , ηm1 )H 1 (I) (ηk2 , ηm2 )L2 (I) +(ηk1 , ηm1 )L2 (I) (ηk2 , ηm2 )H 1 (I)
0
0
0
Kronecker structure: Sηp = Mp ⊗ Ip + Ip ⊗ Mp (max p polyn degree)
Note: If higher order polynomial used, then Sηp simply expands
(augmented Mp )
3
The application. II
Adaptive Legendre-Galerkin discretizations for PDEs:
• Inner product:
v=
X
kvk2H 1 = vˆT Sη vˆ
vˆk ηk ,
0
• (Full!) Orthonormalization:
{Φk } orth basis,
X
v=
v˜k Φk ,
kvk2H 1 = v˜T GT Sη (G˜
v ) = v˜T v˜
0
with G = L−1 where Sη = LLT
• (Cheap!) Quasi-orthonormalization:
{Ψk } quasi-orth basis,
X
ˇ T Sη (G˜
ˇ v ) ≈ v˜T D˜
v=
vˇk Ψk ,
kvk2H 1 = v˜T G
v
0
ˇ very sparse version of G,
G
D diagonal
ˇ exist? ...Analyze sparsity of S −1
Q: Does such a G
η
4
The application. II
Adaptive Legendre-Galerkin discretizations for PDEs:
• Inner product:
v=
X
kvk2H 1 = vˆT Sη vˆ
vˆk ηk ,
0
• (Full!) Orthonormalization:
{Φk } orth basis,
X
v ) = v˜T v˜
v=
v˜k Φk ,
kvk2H 1 = v˜T GT Sη (G˜
0
with G = L−1 where Sη = LLT
• (Cheap!) Quasi-orthonormalization:
{Ψk } quasi-orth basis,
X
ˇ T Sη (G˜
ˇ v ) ≈ v˜T D˜
v
v=
vˇk Ψk ,
kvk2H 1 = v˜T G
0
ˇ very sparse version of G,
G
D diagonal
ˇ exist? ...Analyze sparsity of S −1
Q: Does such a G
η
5
The application. II
Adaptive Legendre-Galerkin discretizations for PDEs:
• Inner product:
v=
X
kvk2H 1 = vˆT Sη vˆ
vˆk ηk ,
0
• (Full!) Orthonormalization:
{Φk } orth basis,
X
v ) = v˜T v˜
v=
v˜k Φk ,
kvk2H 1 = v˜T GT Sη (G˜
0
with G = L−1 where Sη = LLT
• (Cheap!) Quasi-orthonormalization:
{Ψk } quasi-orth basis,
X
ˇ T Sη (G˜
ˇ v ) ≈ v˜T D˜
v
v=
vˇk Ψk ,
kvk2H 1 = v˜T G
0
ˇ very sparse version of G,
G
D diagonal
ˇ exist? ...Analyze sparsity of S −1
Q: Does such a G
η
6
The application. II
Adaptive Legendre-Galerkin discretizations for PDEs:
• Inner product:
v=
X
kvk2H 1 = vˆT Sη vˆ
vˆk ηk ,
0
• (Full!) Orthonormalization:
{Φk } orth basis,
X
v ) = v˜T v˜
v=
v˜k Φk ,
kvk2H 1 = v˜T GT Sη (G˜
0
with G = L−1 where Sη = LLT
• (Cheap!) Quasi-orthonormalization:
{Ψk } quasi-orth basis,
X
ˇ T Sη (G˜
ˇ v ) ≈ v˜T D˜
v
v=
vˇk Ψk ,
kvk2H 1 = v˜T G
0
ˇ very sparse version of G,
G
D diagonal
ˇ exist? ...Analyze sparsity of S −1
Q: Does such a G
η
7
The stiffness matrix
S := M ⊗ In + In ⊗ M,
with M symmetric and positive definite, banded with bandwidth b
• Finite differences: M is second order operator in one space
dimension (b = 1)
⇒
for instance, S: 2D Laplace operator
• Legendre Spectral methods: M spd, nonconstant (b = 1)
• ...
More generally,
Sg := M1 ⊗ In + In ⊗ M2 ,
with M1 6= M2 , banded, with not necessarily the same dimensions
8
The inverse of the 2D Laplace matrix on the unit square
S := M ⊗ In + In ⊗ M,
M = tridiag(−1, 2, −1)
Sparsity pattern:
0
0
10
10
20
20
30
30
40
40
50
50
60
60
70
70
80
80
90
90
100
100
0
10
20
30
40
50
60
nz = 460
70
80
90
100
0
10
20
30
40
50
60
nz = 9380
S −1
Matrix S
9
70
80
90
100
The inverse of the 2D Laplace matrix on the unit square
S := M ⊗ In + In ⊗ M,
M = tridiag(−1, 2, −1)
Sparsity pattern:
0
10
0.7
20
0.6
30
0.5
40
0.4
50
0.3
60
0.2
0.1
70
0
100
80
80
90
100
60
80
60
40
100
0
10
20
30
40
50
60
nz = 460
70
80
90
40
20
100
20
0
0
|(S −1 )ij |
S
10
The exponential decay of the entries of S −1
The classical bound (Demko, Moss & Smith):
If S spd is banded with bandwidth b, then
|(S
−1
)ij | ≤ γq
|i−j|
b
where
κ: condition number of S
√
κ−1
q := √
<1
κ+1
γ := max{λmin (S)
−1
√
(1 + κ)2
, γˆ }, and γˆ =
2λmax (S)
(λmin (·), λmax (·) smallest and largest eigenvalues of the given symmetric matrix)
Many contributions: Bebendorf, Hackbusch, Benzi, Boito, Razouk, Golub, Tuma,
Concus, Meurant, Mastronardi, Ng, Tyrtyshnikov, Nabben, ...
11
The true decay
1
column 35 of S−1
0.9
0.7
0.8
0.6
0.7
magnitude
0.5
0.4
0.3
0.2
0.1
0.6
0.5
0.4
0.3
0
100
0.2
80
100
60
80
0.1
60
40
0
0
40
20
20
0
... a very peculiar pattern
20
40
60
row index
0
⇒ much higher sparsity
12
80
100
Where do the repeated peaks come from?
For S = M ⊗ In + In ⊗ M ∈ Rn
2
×n2
xt := (S −1 ):,t = S −1 et
:
⇔
Solve : Sxt = et
Let
Xt ∈ Rn×n be such that xt = vec(Xt )
Et ∈ Rn×n be such that et = vec(Et )
Then
Sxt = et
⇔
M X t + X t M = Et
13
Where do the repeated peaks come from?
For S = M ⊗ In + In ⊗ M ∈ Rn
2
×n2
xt := (S −1 ):,t = S −1 et
:
⇔
Solve : Sxt = et
Let
Xt ∈ Rn×n be such that xt = vec(Xt )
Et ∈ Rn×n be such that et = vec(Et )
Then
Sxt = et
⇔
M X t + X t M = Et
14
For S the 2D Laplace operator, t = 1, . . . , n2
t = 35,
Sxt = et
⇔
1
M X t + X t M = Et
0.7
0.6
0.8
0.5
0.6
0.4
0.3
0.4
0.2
0.2
0.1
0
10
0
10
8
10
6
6
4
2
0
0
matrix Et
6
4
2
2
0
10
8
4
4
2
8
6
8
and
0
matrix Xt
Et has only one nonzero element
Lexicographic order: (Et )ij , j = ⌊(t − 1)/n⌋ + 1, i = tn⌊(t − 1)/n⌋
15
For S the 2D Laplace operator, t = 1, . . . , n2
t = 35,
Sxt = et
⇔
1
M X t + X t M = Et
0.7
0.6
0.8
0.5
0.6
0.4
0.3
0.4
0.2
0.2
0.1
0
10
0
10
8
10
6
6
4
2
0
0
matrix Et
6
4
2
2
0
10
8
4
4
2
8
6
8
and
0
matrix Xt
Et has only one nonzero element
Lexicographic order: (Et )ij , j = ⌊(t − 1)/n⌋ + 1, i = tn⌊(t − 1)/n⌋
16
0.7
0.6
0.7
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.1
0.2
0
10
8
0.1
10
6
8
6
4
0
4
2
0
10
20
30
40
Left: Row of S −1
50
component
60
70
80
90
100
j
2
0
0
i
Right: same row on the grid
17
0.7
0.6
0.7
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.1
0.2
0
10
8
0.1
10
6
8
6
4
0
4
2
0
10
20
30
40
Left: Row of S −1
50
component
60
70
80
90
100
j
2
0
0
i
Right: same row on the grid
18
0.7
0.6
0.7
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.1
0.2
0
10
8
0.1
10
6
8
6
4
0
4
2
0
10
20
30
40
Left: Row of S −1
50
component
60
70
80
90
100
j
2
0
0
i
Right: same row on the grid
19
0.7
0.6
0.7
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.1
0.2
0
10
8
0.1
10
6
8
6
4
0
4
2
0
10
20
30
40
Left: Row of S −1
50
component
60
70
80
90
100
j
2
0
0
i
Right: same row on the grid
20
0.7
0.6
0.7
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.1
0.2
0
10
8
0.1
10
6
8
6
4
0
4
2
0
10
20
30
40
Left: Row of S −1
50
component
60
70
80
90
100
j
2
0
0
i
Right: same row on the grid
21
0.7
0.6
0.7
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.1
0.2
0
10
8
0.1
10
6
8
6
4
0
4
2
0
10
20
30
40
Left: Row of S −1
50
component
60
70
80
90
100
j
2
0
0
i
Right: same row on the grid
22
Resolving the entry indexing using M Xt + Xt M = Et
(S −1 )k,t = (S −1 )ℓ+n(m−1),t = e⊤
ℓ Xt e m ,
ℓ, m ∈ {1, . . . , n}
⇒ All the elements of the t-th column, (S −1 ):,t , are obtained by
varying m, ℓ ∈ {1, . . . , n}
From the Lyapunov equation theory,
Z ∞
1
Xt =
(ıωI + M )−1 Et (ıωI + M )−∗ dω
2π −∞
with Et = ei e⊤
j , j = ⌊(t − 1)/n⌋ + 1, i = t − n⌊(t − 1)/n⌋
Therefore,
e⊤
ℓ Xt e m
1
=
2π
Z
∞
−∞
−1
−∗
e⊤
ei e⊤
em dω
ℓ (ıωI + M )
j (ıωI + M )
23
Resolving the entry indexing using M Xt + Xt M = Et
(S −1 )k,t = (S −1 )ℓ+n(m−1),t = e⊤
ℓ Xt e m ,
ℓ, m ∈ {1, . . . , n}
⇒ All the elements of the t-th column, (S −1 ):,t , are obtained by
varying m, ℓ ∈ {1, . . . , n}
From the Lyapunov equation theory,
Z ∞
1
Xt =
(ıωI + M )−1 Et (ıωI + M )−∗ dω
2π −∞
with Et = ei ej⊤ , j = ⌊(t − 1)/n⌋ + 1, i = t − n⌊(t − 1)/n⌋
Therefore,
e⊤
ℓ Xt e m
1
=
2π
Z
∞
−∞
−1
−∗
e⊤
ei e⊤
em dω
ℓ (ıωI + M )
j (ıωI + M )
24
Qualitative bounds
Let κ = λmax /λmin = cond(M )
i)Assume ℓ, i, m, j : ℓ 6= i, m 6= j. n2 := |ℓ − i| + |m − j| − 2 > 0
√
κ2 + 1 1
−1
√ .
|(S )k,t | ≤
2λmin
n2
ii)Assume ℓ, i, m, j : ℓ = i or m = j. n1 := |ℓ − i| + |m − j| − 1 > 0
√
κ κ2 + 1 1
−1
√ .
|(S )k,t | ≤
2
n1
1
0.7
0.6
0.8
0.5
0.6
0.4
0.3
0.4
0.2
0.2
0.1
0
10
0
10
8
8
10
6
4
2
0
0
(i, j)
4
2
2
0
8
6
4
4
2
10
6
8
6
and
25
0
(ℓ, m)
Examples. Symmetric positive definite matrix
M = tridiag(−0.5, 2, −0.5) ∈ R10×10
1.4
−1
column id of S
new estimate
Demko etal
1.2
magnitude
1
0.8
0.6
0.4
0.2
0
0
20
40
60
row index
26
80
100
Examples. Legendre stiffness matrix (scaled to have peak equal to 1)
M = tridiag(δk , γk , δk )
1
column 35 of S−1
new estimate
0.9
γk
=
0.8
magnitude
0.7
δk
0.6
k = 1, . . . , n, and
−1
p
(4k + 1) (4k − 1)(4k + 3)
k = 1, . . . , n − 1
0.5
0.4
0.3
0.2
0.1
0
0
=
2
(4k − 3)(4k + 1)
20
40
60
80
100
row index
27
Connections to point-wise estimates for discrete Laplacian
For the discrete Green function Gh on the discrete d-dimensional grid
Rh , there exist constants h0 and C such that for h ≤ h0 , x, y ∈ Rh ,

C
 C log
if d = 2
|x−y|+h
Gh (x, y) ≤
C

if d ≥ 3
d−2
(|x−y|+h)
(Bramble & Thomee, ’69)
Our estimate: entries depend on inverse square root of the distance!
28
Explored generalizations
• M spd of bandwidth b > 1
• S = M1 ⊗ I + I ⊗ M2 , M1 6= M2
• M1 , M2 of different bandwidth
• LLT = S, then L−1 (lower triang.) has same sparsity pattern
references:
C. Canuto, V. Simoncini and M. Verani, LAA, v.452, 2014.
C. Canuto, V. Simoncini and M. Verani, Adaptive Legendre-Galerkin methods, in
preparation, 2014.
29

Download Report