Risoluzione di equazioni matriciali con

matrici infinite quasi-Toeplitz

Leonardo Robol, ISTI-CNR, Pisa, Italy

<leonardo.robol@isti.cnr.it>

joint work with:

D. A. Bini, B. Meini (UniPi)

S. Massei (EPFL)

Montecatini, 16 Feb 2018

1 / 26

GNCS Project / Participants

Project: “Metodi numerici avanzati per equazioni e funzioni di

matrici con struttura”, coordinated by B. Meini.

• Dario Bini

• Gianna Del Corso

• Dario Fasino

• Luca Gemignani

• Bruno Iannazzo

• Nicola Mastronardi

• Federico Poloni

• Valeria Simoncini

• Massimiliano Fasi

• Stefano Massei

• Davide Palitta

• Leonardo Robol

2 / 26

GNCS Project / Participants

Project: “Metodi numerici avanzati per equazioni e funzioni di

matrici con struttura”, coordinated by B. Meini.

• Dario Bini

• Gianna Del Corso

• Dario Fasino

• Luca Gemignani

• Bruno Iannazzo

• Nicola Mastronardi

• Federico Poloni

• Valeria Simoncini

• Massimiliano Fasi

• Stefano Massei

• Davide Palitta

• Leonardo Robol

2 / 26

Random walks on a quadrant

We consider random walks on N× N:

• Allowed moves only to adjacent states.

• Probabilities of going up/down/left/right are eventually

independent of position (i , j).

• The problem can be stated also on the infinite strip

{1, . . . ,m} × N.

3 / 26

The link with Toeplitz matrices

Several queuing problems can be recasted in this framework.

The simpler 1D case considers movements on a line:

The probability transition matrix is as follows:

P =

ρ0 + ρ−1 ρ1ρ−1 ρ0 ρ1

. . .. . .

. . .

,where ρ−1 is the probability of moving left, ρ1 right, and

ρ0 = 1− ρ−1 − ρ1: a quasi-Toeplitz semi-infinite matrix.4 / 26

Toeplitz matrices

• Toeplitz matrices have constant diagonals.

• In our case, the entry (1, 1) is different: we have a

top-left correction (quasi-Toeplitz structure).

• In the 2D case, the matrix P is (quasi-)TBBT

(Toeplitz-block block-Toeplitz):

P =

A0 A1

A−1 A0 A1

. . .. . .

. . .

,with Ai , A0 quasi-Toeplitz (and semi-infinite).

5 / 26

Computing the invariant vector

The steady state probability vector solving πTP = πT can be

recovered by finding the minimal non-negative solutions to

A−1 + A0G + A1G2 = G , R2A−1 + RA0 + A1 = R

• Cyclic reduction (CR) is a matrix iteration that converges

to the correct solutions.

• Needs to perform multiplications, sum, and inversions.

• Complexity O(m3) if the matrices Ai are m ×m — for

instance a random walk on {1, . . . ,m} × N.

6 / 26

A quick look at the matrix iteration

A(h+1)0 = A

(h)0 + A

(h)−1S

(h)A(h)1 + A

(h)1 S (h)A

(h)−1

A(h+1)−1 = A

(h)−1S

(h)A(h)−1

A(h+1)1 = A

(h)1 S (h)A

(h)1

S (h) = (I − A(h)0 )−1

The initial matrices are A(0)i := Ai , so they are Toeplitz. Is the

structure preserved?

• Experimentally, approximately, up to a top-left correction.

• Non-trivial to exploit in practice. In fact, typical

implementation ignore the structure and require O(m3)

flops =⇒ m =∞ is unsolved.

7 / 26

A quick look at the matrix iteration

A(h+1)0 = A

(h)0 + A

(h)−1S

(h)A(h)1 + A

(h)1 S (h)A

(h)−1

A(h+1)−1 = A

(h)−1S

(h)A(h)−1

A(h+1)1 = A

(h)1 S (h)A

(h)1

S (h) = (I − A(h)0 )−1

The initial matrices are A(0)i := Ai , so they are Toeplitz. Is the

structure preserved?

• Experimentally, approximately, up to a top-left correction.

• Non-trivial to exploit in practice. In fact, typical

implementation ignore the structure and require O(m3)

flops =⇒ m =∞ is unsolved.

7 / 26

Basic facts about Toeplitz matrices

T (a(z)) =

a0 a1 . . . . . .

a−1 a0 a1... a−1

. . .. . .

.... . .

. . .

, a(z) =∑j∈Z

ajzj

• a(z) is in the Wiener algebra W if∑j∈Z

|aj | <∞

• a(z) ∈ W1 ⇐⇒ a′(z) ∈ W ⇐⇒∑

j∈Z |i · ai | <∞.

8 / 26

Hankel matrices

Given a power series f (z) =∑∞

j=0 fjzj , we have

H(f (z)) =

f1 f2 f3 . . .

f2 f3

f3. . .

...

the Hankel matrix defined by f (z). Often, we use the notation

a+(z) =∑j≥1

ajzj , a−(z) =

∑j≥1

a−jzj .

9 / 26

Toeplitz matrices are not an algebra

The power series in W and W1 form an algebra. However,

T (a(z)) · T (b(z)) 6= T (c(z)), c(z) = a(z)b(z).

The link between Toeplitz matrices and Laurent series is

indeed an isomorphism for bi-infinite matrices.

Theorem (Gohberg–Feldman)

Let a(z), b(z) ∈ W . Then,

T (a)T (b) = T (c)− H(a−)H(b+),

H(f ) Hankel matrix. The correction H(a−)H(b+) is a

compact operator on `2.

Toeplitz matrices are an algebra up to compact corrections.

10 / 26

Toeplitz matrices are not an algebra

The power series in W and W1 form an algebra. However,

T (a(z)) · T (b(z)) 6= T (c(z)), c(z) = a(z)b(z).

The link between Toeplitz matrices and Laurent series is

indeed an isomorphism for bi-infinite matrices.

Theorem (Gohberg–Feldman)

Let a(z), b(z) ∈ W . Then,

T (a)T (b) = T (c)− H(a−)H(b+),

H(f ) Hankel matrix. The correction H(a−)H(b+) is a

compact operator on `2.

Toeplitz matrices are an algebra up to compact corrections.

10 / 26

Toeplitz matrices are not an algebra

The power series in W and W1 form an algebra. However,

T (a(z)) · T (b(z)) 6= T (c(z)), c(z) = a(z)b(z).

The link between Toeplitz matrices and Laurent series is

indeed an isomorphism for bi-infinite matrices.

Theorem (Gohberg–Feldman)

Let a(z), b(z) ∈ W . Then,

T (a)T (b) = T (c)− H(a−)H(b+),

H(f ) Hankel matrix. The correction H(a−)H(b+) is a

compact operator on `2.

Toeplitz matrices are an algebra up to compact corrections.10 / 26

The Hankel correction

Since a(z), b(z) are in W , then aj , bj → 0 for j →∞, so:

• T (a), T (b) are numerically banded.

• H(a−),H(b+) have the support in the top-left corner.

Therefore, T (a) · T (b) = T (c)− H(a−)H(b+) is numerically:

= +

11 / 26

A new class of matrices

Consider the norm ‖E‖F :=∑

i ,j |Eij |, and

F := {E ∈ RN×N | ‖E‖F <∞}.

Definition

A = T (a(z)) + Ea is quasi-Toeplitz (QT), if:

• a(z) ∈ W1 ⊆ W .

• Ea satisfies ‖Ea‖F <∞.

We define ‖A‖QT := ‖a(z)‖W + ‖a′(z)‖W + ‖Ea‖F .

• QT-matrices form a (computationally-friendly) algebra.

• ‖Ea‖F <∞ =⇒ Ea compact operator on `2.

12 / 26

A new class of matrices

Consider the norm ‖E‖F :=∑

i ,j |Eij |, and

F := {E ∈ RN×N | ‖E‖F <∞}.

Definition

A = T (a(z)) + Ea is quasi-Toeplitz (QT), if:

• a(z) ∈ W1 ⊆ W .

• Ea satisfies ‖Ea‖F <∞.

We define ‖A‖QT := ‖a(z)‖W + ‖a′(z)‖W + ‖Ea‖F .

• QT-matrices form a (computationally-friendly) algebra.

• ‖Ea‖F <∞ =⇒ Ea compact operator on `2.

12 / 26

A new class of matrices

Consider the norm ‖E‖F :=∑

i ,j |Eij |, and

F := {E ∈ RN×N | ‖E‖F <∞}.

Definition

A = T (a(z)) + Ea is quasi-Toeplitz (QT), if:

• a(z) ∈ W1 ⊆ W .

• Ea satisfies ‖Ea‖F <∞.

We define ‖A‖QT := ‖a(z)‖W + ‖a′(z)‖W + ‖Ea‖F .

• QT-matrices form a (computationally-friendly) algebra.

• ‖Ea‖F <∞ =⇒ Ea compact operator on `2.

12 / 26

Numerical representation of a QT matrix

We can represent QT matrices with a finite number of

parameters and arbitrary accuracy.

If A = T (a) + Ea:

• T (a) can be stored by a (truncated) approximation of its

symbol, since∑|i |>j |ai | <∞;

• Ea ≈ UV T by means of SVD; U ,V have compact

support and σj(Ea)→ 0, since Ea is compact and so

“numerically low-rank”.

In our case we use ‖Ea‖F <∞ to pick a top-left block as

approximation to Ea.

13 / 26

QT is an algebra

We can check that:

• A,B ∈ QT =⇒ A + B ∈ QT :

T (a) + Ea + T (b) + Eb = T (a + b) + (Ea + Eb).

• A,B ∈ QT =⇒ AB ∈ QT :

(T (a) + Ea) · (T (b) + Eb) = T (ab)+

(T (a)Eb + EaT (b) + EaEb − H(a−)H(b+))

• A = T (a) + Ea implies A−1 = T (a−1) + ?

14 / 26

QT is an algebra

We can check that:

• A,B ∈ QT =⇒ A + B ∈ QT :

T (a) + Ea + T (b) + Eb = T (a + b) + (Ea + Eb).

• A,B ∈ QT =⇒ AB ∈ QT :

(T (a) + Ea) · (T (b) + Eb) = T (ab)+

(T (a)Eb + EaT (b) + EaEb − H(a−)H(b+))

• A = T (a) + Ea implies A−1 = T (a−1) + ?

14 / 26

QT is an algebra

We can check that:

• A,B ∈ QT =⇒ A + B ∈ QT :

T (a) + Ea + T (b) + Eb = T (a + b) + (Ea + Eb).

• A,B ∈ QT =⇒ AB ∈ QT :

(T (a) + Ea) · (T (b) + Eb) = T (ab)+

(T (a)Eb + EaT (b) + EaEb − H(a−)H(b+))

• A = T (a) + Ea implies A−1 = T (a−1) + ?

14 / 26

QT is an algebra (cont’d)

If a(z), b(z) ∈ W1 and Ea,Eb ∈ F , then the corrections

obtained by sum and multiplication are in F .

Proof: some (not so interesting) computations. a(z) ∈ W1

plays an important role.

And the inverse?

15 / 26

QT is an algebra (cont’d)

If a(z), b(z) ∈ W1 and Ea,Eb ∈ F , then the corrections

obtained by sum and multiplication are in F .

Proof: some (not so interesting) computations. a(z) ∈ W1

plays an important role.

And the inverse?

15 / 26

QT is an algebra (cont’d)

If a(z), b(z) ∈ W1 and Ea,Eb ∈ F , then the corrections

obtained by sum and multiplication are in F .

Proof: some (not so interesting) computations. a(z) ∈ W1

plays an important role.

And the inverse?

15 / 26

Computing the inverse

For the inverse of T (a), we need a few steps:

• Factor a(z) = u(z)`(z−1), where u, ` are power series

(Wiener-Hopf factorization).

• This corresponds to the UL factorization

T (a) = T (u(z)) · T (`(z−1))

• Triangular Toeplitz matrices have Toeplitz inverses:

T (a(z))−1 = T (`−1(z−1))T (u−1(z)),

and we can use what we know about multiplications.

16 / 26

Computational remarks

We can operate on quasi-Toeplitz matrices combining:

• Operations/Factorizations on power series (FFT-based);

• Toeplitz-vector multiplications (again FFT-based);

• Compression of low-rank matrices of the form

H(a−)H(b+).

For the last item, fast matvec is available, so we can run either

Lanczos or random sampling to compress the products of

Hankel matrices.

Ranks of the corrections are small in practice: 10–20 is a

typical value in applications.

17 / 26

Computational remarks

We can operate on quasi-Toeplitz matrices combining:

• Operations/Factorizations on power series (FFT-based);

• Toeplitz-vector multiplications (again FFT-based);

• Compression of low-rank matrices of the form

H(a−)H(b+).

For the last item, fast matvec is available, so we can run either

Lanczos or random sampling to compress the products of

Hankel matrices.

Ranks of the corrections are small in practice: 10–20 is a

typical value in applications.

17 / 26

Keeping the rank bounded

When we perform arithmetic operations, the rank of the

correction increases in the representation. This can be kept

low by recompression.

Assume Ea = UaVTa , Ua,Va with k columns.

• [QU ,RU ] = qr(Ua), and [QV ,RV ] = qr(Va).

• [U1,Σ1,V1] ≈ svd(RURTV ) (truncated SVD).

• Ea ≈ QUU1

√Σ1 · (QVV1

√Σ1).

The new representation has k ′ < k columns.

The cost of the compression is O(nk2 + k3) flops, where n is

the support of Ea.

18 / 26

Keeping the rank bounded

When we perform arithmetic operations, the rank of the

correction increases in the representation. This can be kept

low by recompression.

Assume Ea = UaVTa , Ua,Va with k columns.

• [QU ,RU ] = qr(Ua), and [QV ,RV ] = qr(Va).

• [U1,Σ1,V1] ≈ svd(RURTV ) (truncated SVD).

• Ea ≈ QUU1

√Σ1 · (QVV1

√Σ1).

The new representation has k ′ < k columns.

The cost of the compression is O(nk2 + k3) flops, where n is

the support of Ea.

18 / 26

Keeping the rank bounded

When we perform arithmetic operations, the rank of the

correction increases in the representation. This can be kept

low by recompression.

Assume Ea = UaVTa , Ua,Va with k columns.

• [QU ,RU ] = qr(Ua), and [QV ,RV ] = qr(Va).

• [U1,Σ1,V1] ≈ svd(RURTV ) (truncated SVD).

• Ea ≈ QUU1

√Σ1 · (QVV1

√Σ1).

The new representation has k ′ < k columns.

The cost of the compression is O(nk2 + k3) flops, where n is

the support of Ea.

18 / 26

Keeping the rank bounded

When we perform arithmetic operations, the rank of the

correction increases in the representation. This can be kept

low by recompression.

Assume Ea = UaVTa , Ua,Va with k columns.

• [QU ,RU ] = qr(Ua), and [QV ,RV ] = qr(Va).

• [U1,Σ1,V1] ≈ svd(RURTV ) (truncated SVD).

• Ea ≈ QUU1

√Σ1 · (QVV1

√Σ1).

The new representation has k ′ < k columns.

The cost of the compression is O(nk2 + k3) flops, where n is

the support of Ea.

18 / 26

Keeping the rank bounded

When we perform arithmetic operations, the rank of the

correction increases in the representation. This can be kept

low by recompression.

Assume Ea = UaVTa , Ua,Va with k columns.

• [QU ,RU ] = qr(Ua), and [QV ,RV ] = qr(Va).

• [U1,Σ1,V1] ≈ svd(RURTV ) (truncated SVD).

• Ea ≈ QUU1

√Σ1 · (QVV1

√Σ1).

The new representation has k ′ < k columns.

The cost of the compression is O(nk2 + k3) flops, where n is

the support of Ea.

18 / 26

Keeping the rank low (pictorial version)

= +

=

19 / 26

Keeping the rank low (pictorial version)

= +

=

19 / 26

A Banach algebra

We can say even more about QT .

Theorem

The set of QT-matrices, endowed with the norm

‖T (a) + Ea‖ := ‖a‖W + ‖a′‖W + ‖Ea‖F ,

is a Banach algebra.

This means that it is an algebra, it is a Banach space, and that

‖AB‖ ≤ ‖A‖ · ‖B‖.

20 / 26

Back to our matrix iteration

Recall we wanted to use cyclic reduction to solve

A−1 + A0G + A1G2 = G , R2A−1 + RA0 + A1 = R

Theorem

Let Ai ∈ B, a Banach algebra. If

ϕ(z) = z−1A1 + (A0 − I ) + zA1

is invertible on {t−1 < |z | < t} for t > 1, then CR converges

quadratically and G , R can be obtained from the limits of A(h)i .

21 / 26

Good news

• The hypotheses of the previous results are always satisfied

in the Markov chain setting.

• QT-matrices are a Banach algebra so cyclic reduction

works in this class, and we know how to compute the

iteration.

• The semi-infinite solutions G ,R , will also be represented

in the QT format.

• We have a MATLAB toolbox that makes it easy to play

with these matrices: You can google “cqt-toolbox”.

22 / 26

A simple example

Consider a tandem Jackson queue. It has been shown that

finite truncation may not approximate the infinite-dimensional

case behavior1.

1Motyer, A.J. and Taylor, P.G., 2006. Decay rates for

quasi-birth-and-death processes with countably many phases and

tridiagonal block generators. Advances in applied probability, 38(2),

pp.522-544.23 / 26

Explicit solution

We can compute the steady state vector π using CR:

Problem Time (s) Residue ‖·‖∞ Band Support Rank

1 2.61 2.02 · 10−13 561 541 8

2 2.91 9.09 · 10−13 561 555 8

3 0.28 2.02 · 10−13 143 89 8

4 2.32 1.77 · 10−13 463 481 9

5 0.47 1.93 · 10−13 233 148 9

6 7.96 1.16 · 10−12 455 462 10...

......

......

...

The residue is ‖πTP − πT‖∞, and each problem corresponds

to different rates and parameters. Ranks and support refer to

the corrections in G ,R , and the band to their Toeplitz part.

24 / 26

Conclusions and future outlook

• Computing with infinite matrices might be fun, if you

have the right structure.

• Some hypotheses on the norms can be relaxed — this is

work in progress.

• We can do much more fancy things: compute matrix

functions, solve other kinds of matrix equations, . . . . This

has applications to other fields, such as finance.

• Finite case handled as well.

• A MATLAB toolbox is available for you to try: feedback

is very welcome!

25 / 26

Thank you for your attention!

26 / 26