Sherman-Morrison-Woodbury formula for Sylvester and T -Sylvester - - PowerPoint PPT Presentation

Sherman-Morrison-Woodbury formula for Sylvester and T-Sylvester equation

Ivana Kuzmanović

Department of Mathematics, University of Osijek, Croatia

I.Kuzmanović (University of Osijek, Croatia) 1 / 9

Problem

We consider the Sylvester equation (A0 + U1C1V1)X + X(B0 + U2C2V2) = E, (1) where A0 ∈ Rm×m, U1, V T

1 ∈ Rm×r1, C1 ∈ Rr1×r1 , B0 ∈ Rn×n,

U2, V T

2 ∈ Rn×r2, C2 ∈ Rr2×r2, X, E ∈ Rm×n.

I.Kuzmanović (University of Osijek, Croatia) 2 / 9

Problem

We consider the Sylvester equation (A0 + U1C1V1)X + X(B0 + U2C2V2) = E, (1) where A0 ∈ Rm×m, U1, V T

1 ∈ Rm×r1, C1 ∈ Rr1×r1 , B0 ∈ Rn×n,

U2, V T

2 ∈ Rn×r2, C2 ∈ Rr2×r2, X, E ∈ Rm×n.

We assume r1, r2 ≪ m, n and that C1, C2 are regular matrices

I.Kuzmanović (University of Osijek, Croatia) 2 / 9

Problem

We consider the Sylvester equation (A0 + U1C1V1)X + X(B0 + U2C2V2) = E, (1) where A0 ∈ Rm×m, U1, V T

1 ∈ Rm×r1, C1 ∈ Rr1×r1 , B0 ∈ Rn×n,

U2, V T

2 ∈ Rn×r2, C2 ∈ Rr2×r2, X, E ∈ Rm×n.

We assume r1, r2 ≪ m, n and that C1, C2 are regular matrices Goal: efficient method for solving equation (1) many times with different update matrices U1, C1, V1, U2, C2, V2.

I.Kuzmanović (University of Osijek, Croatia) 2 / 9

Basic idea

Using Kronecker product and vectorization operator, as well as equality vec(ABC) = (CT ⊗ A) vec(B), equation (1) can be rewritten as mn × mn matrix equation (In ⊗ (A0 + U1C1V1) + (B0 + U2C2V2)T ⊗ Im) vec(X) = vec(E)

I.Kuzmanović (University of Osijek, Croatia) 3 / 9

Basic idea

Using Kronecker product and vectorization operator, as well as equality vec(ABC) = (CT ⊗ A) vec(B), equation (1) can be rewritten as mn × mn matrix equation (In ⊗ (A0 + U1C1V1) + (B0 + U2C2V2)T ⊗ Im) vec(X) = vec(E) and further as equation (L0 + UCV) vec(X) = vec(E)

I.Kuzmanović (University of Osijek, Croatia) 3 / 9

Basic idea

Using Kronecker product and vectorization operator, as well as equality vec(ABC) = (CT ⊗ A) vec(B), equation (1) can be rewritten as mn × mn matrix equation (In ⊗ (A0 + U1C1V1) + (B0 + U2C2V2)T ⊗ Im) vec(X) = vec(E) and further as equation (L0 + UCV) vec(X) = vec(E) where

L0 = In ⊗ A0 + BT

0 ⊗ Im ∈ Rmn×mn,

U =

• In ⊗ U1 V T

2 ⊗ Im

• ∈ Rmn×(r1n+r2m),

C = In ⊗ C1 C T

2 ⊗ Im

• ∈ R(r1n+r2m)×(r1n+r2m),

V = In ⊗ V1 UT

2 ⊗ Im

• ∈ R(r1n+r2m)×mn.

I.Kuzmanović (University of Osijek, Croatia) 3 / 9

Basic idea

(A0 + U1C1V1)X + X(B0 + U2C2V2) = E

• (L0 + UCV) vec(X) = vec(E)

I.Kuzmanović (University of Osijek, Croatia) 4 / 9

Basic idea

(A0 + U1C1V1)X + X(B0 + U2C2V2) = E

• (L0 + UCV) vec(X) = vec(E)

After applying the Sherman-Morrison-Woodbury formula on the right hand side of the previous equation, we obtain vec(X) = (L−1 − L−1

0 U(C−1 + VL−1 0 U)−1VL−1 0 ) vec (E).

L0 = In ⊗ A0 + BT

0 ⊗ Im ∈ Rmn×mn,

U =

• In ⊗ U1 V T

2 ⊗ Im

• ∈ Rmn×(r1n+r2m),

C = In ⊗ C1 C T

2 ⊗ Im

• ∈ R(r1n+r2m)×(r1n+r2m),

V = In ⊗ V1 UT

2 ⊗ Im

• ∈ R(r1n+r2m)×mn.

I.Kuzmanović (University of Osijek, Croatia) 4 / 9

Special case

Special case that we have considered: optimization of the solution of the parameter-dependent Sylvester equation (A0 − vU1V1)X(v) + X(v)(B0 − vU2V2) = E. We want to optimize the solution X(v) of the above Sylvester equation (i.e. to determine optimal parameter v ∈ R) with respect to some

• ptimization criterion.

I.Kuzmanović (University of Osijek, Croatia) 5 / 9

Special case

Special case that we have considered: optimization of the solution of the parameter-dependent Sylvester equation (A0 − vU1V1)X(v) + X(v)(B0 − vU2V2) = E. We want to optimize the solution X(v) of the above Sylvester equation (i.e. to determine optimal parameter v ∈ R) with respect to some

• ptimization criterion.

Optimization criteria: tr(X(v)) → min for m = n, X(v)F → min

I.Kuzmanović (University of Osijek, Croatia) 5 / 9

Special case

vec(X) = (L−1 − vL−1

0 U(I + v∆)−1VL−1 0 ) vec (E)

∆ = VL−1

0 U

We need efficient methods for solving system L0x = y

L−1 given or A0 ˆ X + ˆ XB0 = ˆ Y where vec(ˆ X) = x,vec( ˆ Y ) = y.

matrix - vector product with U and V

Vy =

• In ⊗ V1

UT

2 ⊗ Im

• vec( ˆ

Y ) =

• (In ⊗ V1) vec( ˆ

Y ) (UT

2 ⊗ Im) vec( ˆ

Y )

• =
• vec(V1 ˆ

Y ) vec( ˆ Y U2)

• ,

solving system with I + v∆

Krylov subspace methods (e.g. Full Orthogonalization Method )

I.Kuzmanović (University of Osijek, Croatia) 6 / 9

Number of operations

Let: O(r1) = O(r2) =: r, m = n Calculation vec(X(v)) for the first value of v: 8(rk + r + k + 2)m2 + O(rmk).

I.Kuzmanović (University of Osijek, Croatia) 7 / 9

Number of operations

Let: O(r1) = O(r2) =: r, m = n Calculation vec(X(v)) for the first value of v: 8(rk + r + k + 2)m2 + O(rmk). For every additional evaluation of vec(X(v)) (4r + 8)m2 + O(rmk).

I.Kuzmanović (University of Osijek, Croatia) 7 / 9

Number of operations

Let: O(r1) = O(r2) =: r, m = n Calculation vec(X(v)) for the first value of v: 8(rk + r + k + 2)m2 + O(rmk). For every additional evaluation of vec(X(v)) (4r + 8)m2 + O(rmk). Additional number of flops for vec(X ′(v)) and vec(X ′′(v)) is for each (4r + 8)m2 + O(rmk) where k ≪ m

I.Kuzmanović (University of Osijek, Croatia) 7 / 9

Flops number growth for algorithms lyap, old, and new for the case of r=1,2,3,4

1000 2000 3000 3500 10

5

10

10

m Number of flops r=1 1000 2000 3000 3500 10

5

10

10

m Number of flops r=2 1000 2000 3000 3500 10

5

10

10

m Number of flops r=3 1000 2000 3000 3500 10

5

10

10

m Number of flops r=4 lyap

• ld

new lyap

• ld

new lyap

• ld

new lyap

• ld

new

I.Kuzmanović (University of Osijek, Croatia) 8 / 9

T-Sylvester equation

Similarly as for Sylvester equation, Sherman-Morrison-Woodbury formula can be used for updating the solution of equation (A0 + U1C1V1)X + X T(B0 + U2C2V2) = E. Corresponding matrices are L′ = I ⊗ A0 + (BT

0 ⊗ I)Π

U′ =

• I ⊗ U1 V T

2 ⊗ I

• C′

=

• I ⊗ C1

CT

2 ⊗ I

• V′

=

• I ⊗ V1

(UT

2 ⊗ I)Π

• ,

where Π vec (X) = vec (X T).

I.Kuzmanović (University of Osijek, Croatia) 9 / 9