Efficient Preconditioners for PDE-Constrained Optimization Problems with a Multi-level Sequentially Semi-Separable Matrix Structure

(1)

REPORT 13-04

Effi ient Pre onditioners for PDE-Constrained

Optimization Problems with a Multi-level

Sequentially SemiSeparable

Matrix Stru ture

YueQiu, Martin B.van Gijzen, Jan-Willem van Wingerden

Mi hel Verhaegen,andCornelis Vuik

ISSN 1389-6520

Reports of the Delft Institute of Applied Mathemati s

(2)

Copyright



2013byDelftInstituteofAppliedMathemati s,Delft,TheNetherlands. Nopartof theJournalmaybereprodu ed, storedinaretrievalsystem,ortransmitted,in

any form orby anymeans, ele troni ,me hani al, photo opying, re ording, orotherwise,

without the prior written permission from Delft Institute of Applied Mathemati s, Delft

(3)

Optimization Problems with a Multi-level Sequentially

SemiSeparable Matrix Stru ture

∗

YueQiu

†

,MartinB.vanGijzen

‡

,Jan-WillemvanWingerden

∗

,

Mi helVerhaegen

∗

,andCornelisVuik

†

Revised on: 10th February 2014

Abstra t

Inthispaper,we onsiderpre onditioningforPDE- onstrainedoptimization

prob-lems. Theunderlyingproblemsyieldalinearsaddle-pointsystem. Westudya lassof

pre onditionersbasedonmultilevelsequentiallysemiseparable(MSSS)matrix

ompu-tations. Thenovelglobalpre onditioneris tomakeuseoftheglobal stru tureof the

saddle-pointsystem, while theblo k pre onditionermakesuse ofthe blo k stru ture

of the saddle-point system. For the novelglobal pre onditioner, itis independent of

theregularizationparameter,whilefortheblo kpre onditioner,thispropertydoesnot

hold. ForthisMSSSmatrix omputationapproa h,modelorderredu tionalgorithms

areessentialto obtainalow omputational omplexity.Westudytwodierentmodel

orderredu tionapproa hes,oneisthenewapproximatebalan edtrun ationalgorithm

with low-rank approximatedGramians and the other is the standard Hankelblo ks

approximation algorithm. We testthe global pre onditionerand theblo k

pre ondi-tionerfortheproblemofoptimal ontrolofthePoissonequationandoptimal ontrolof

the onve tion-diusionequation. Numeri alexperimentsillustratethatboth

pre on-ditionersgivelinear omputational omplexityandtheglobalpre onditioneryieldsthe

fewestnumberofiterationsand omputingtime. Moreover,theapproximatebalan ed

trun ation algorithm onsumesless oating-point operations(ops)thanthe Hankel

blo ksapproximationalgorithm.

Keywords: PDE- onstrainedoptimization,saddle-pointproblem,pre onditioners,

multilevelsequentiallysemiseparablematrix,modelorderredu tion,low-rank

approx-imation

1 Introdu tion

Optimaldesign,optimal ontrolandparameterestimationofsystemsgovernedbypartial

dierentialequations(PDEs)giverisetoa lassofproblemsknownasPDE- onstrained

op-timization. PDE- onstrainedoptimizationproblemshaveawideappli ationsu hasoptimal

ow ontrol[1℄ [2℄, diuseopti al tomography [3℄, and linear(nonlinear) model predi tive

ontrol[4℄. Thesolutionoftheseproblemsisa quiredbysolvingalarge-s alelinearsystem

ofsaddle-pointtype. Mu heorthasbeendedi atedtonde ientiterativesolution

meth-odsforsu hsystems. Someofthemostpopularte hniquesarethe onjugategradient(CG)

[5℄, minimal residual (MINRES) [6℄, generalizedminimal residual(GMRES) and indu ed

∗

Thisresear hissupportedbytheNWOVeniGrant

#

11930"Re ongurableFloatingWindFarms".

†

DelftCenter forSystemandControl,DelftUniversityofTe hnology,2628 CDDelft,theNetherlands,

Y.Qiutudelft.nl, YueCiouGmail. om.

‡

DelftInstituteofAppliedMathemati s,DelftUniversityofTe hnology,2628CDDelft,theNetherlands,

(4)

hoi e of pre onditioners. In this paper,westudy a lass of pre onditionersthat exploits

themultilevelsequentiallysemiseparable(MSSS)stru tureoftheblo ksofthesaddle-point

system.

Semiseparablematri esappearin severaltypesofappli ations,e.g. theeld ofintegral

equations [8℄, Gauss-Markov pro esses [9℄, boundary value problems [10℄ and rational

in-terpolation [11℄. Semiseparablematri es are matri es of whi h allthe sub-matri estaken

from the lower-triangularorthe upper-triangular partare of rank at most 1by [12℄.

Se-quentially semiseparable (SSS) matri es of whi h the o-diagonal blo ks are of low-rank,

notlimited to 1, named by Dewilde et. al. in [13℄ generalize thesemiseparable matri es.

Multi-levelsequentiallysemiseparablegeneralizethesequentiallysemiseparablematri esto

themulti-dimensional ases. Systemsthatarisefromthedis retizationof1Dpartial

dier-entialequations typi allyhaveanSSS stru ture. Dis retizationofhigher dimensional(2D

or3D)partial dierentialequations giveriseto matri esthathaveanMSSSstru ture[14℄

[15℄. Under the multilevel paradigm, generators that are used to represent a matrixof a

higherhierar hyare themselves multilevel sequentiallysemiseparable of alowerhierar hy.

The usual one-level sequentially semiseparable matrix is the one of the lowest hierar hy.

Operations likethe matrix inversionand the matrix-matrix multipli ation are still losed

underthisstru ture. The

LU

fa torization analsobeperformedin astru turepreserving way. Thisfa torizationresultsinagrowthoftherankoftheo-diagonalblo ksoftheS hur

omplement. As aresult, the

LU

fa torization is notof linear omputational omplexity. Themodelorderredu tionplaysakeyrolein redu ingtherankoftheo-diagonalblo ks.

Be auseofthemodelorderredu tionoperationbeingperformed,itis possibleto ompute

aninexa t

LU

de ompositionofanMSSSmatrixthat anbeusedasapre onditioner. In [14℄, Gondzio et. al. rst introdu ed the pre onditioningof PDE- onstrained

opti-mizationproblemsbyMSSSmatrix omputations. TheyexploitedtheMSSSmatrix

stru -tureof theblo ks of the saddle-point systemand performed an

LU

fa torization method forMSSSmatri estoapproximatetheS hur omplementofthesaddle-pointsystem. With

theapproximateS hur omplement, onjugategradientmethodwasperformedtosolvethe

pre onditioned saddle-point system blo k-by-blo k. As aforementioned, the model order

redu tionplaysavitalroleinobtainingalinear omputational omplexityofthe

LU

fa tor-ization. In[14℄,Gondzioet. al. usedastandardmodelorderredu tionalgorithm[16℄[13℄to

redu ethe omputational omplexity. Inthispaper,ourwork extends[14℄ inthefollowing

ways. 1)Weproposeanewmodelorderredu tionalgorithmforSSSmatrix omputations

basedonthe orresponden ebetweenlineartime-varying(LTV)systemsandblo ksofSSS

matri es. Thenewmodel orderredu tionalgorithm is motivated by[17℄. In [17℄, the

ap-proximate balan ed trun ationwasaddressedfor themodel orderredu tionof lineartime

invariant(LTI) systems. In thispaper,weextendthat method to thelineartime varying

(LTV) systems. Be ause of the orresponden e between MSSSmatrix and LTV systems,

it is suitable for model order redu tion for MSSS matrix omputations. Compared with

the onventional model order redu tionalgorithms in [13℄ [16℄, theapproximate balan ed

trun ationneedslessoating-pointoperations(ops). 2)Withthesemodelorderredu tion

algorithms,we an omputeaninexa t

LU

fa torizationoftheMSSSmatrixblo ksofthe saddle-pointsystemin linear omputational omplexity. This yieldsblo kpre onditioners

forthesaddle-pointsystemsofthetypethataredes ribedin[18℄ whileonlysingle

pre on-ditioner forthe last blo k of the saddle-point system is studied in [14℄. 3) By permuting

the blo ks of the saddle-point system, we an also ompute an inexa t

LU

fa torization ofthe globalsystem,whi h givesanovelglobal pre onditioner. 4)Besides theproblem of

optimal ontrolofthePoissonequation,wealsostudytheproblemofoptimal ontrolofthe

onve tion-diusionequation. 5)Wealsoextendthesepre onditioningte hniquetothe3D

ases.

Notethatthestandardblo kpre onditionersdependontheregularizationparameter

β

forthePDE- onstrainedoptimizationproblem[19℄. Bypermutingthesaddle-pointsystem

(5)

LU

fa torization of the global system in linear omputational omplexity, whi h is alled the global pre onditioner. Numeri al experiments for the optimal ontrol of the Poisson

equationandthe onve tion-diusionequationillustratethattheperforman eoftheglobal

pre onditionerisindependentoftheregularizationparameter

β

and isalsoindependentof themeshsize.

The stru ture of this paper is as follows: we start with formulating adistributed

op-timal ontrol problem onstrained by PDEs. This problem yields a linear system of the

saddle point type. Demand for e ient pre onditioners to solve this typeof system with

iterative solutionmethods motivatesthis paper. InSe tion 3, we briey givean overview

of somedenitions and the widelyused omputations of MSSSmatri es and then dis uss

theMSSSpre onditioningte hnique. Thenovelmodelorderredu tionalgorithmisalso

de-s ribed. BasedontheMSSSmatrix omputations,weproposethreepre onditionersforthis

saddle-pointproblem,theyarethenovelglobalpre onditioner,thestandardblo k-diagonal

pre onditionerandthestandardblo klower-triangularpre onditioner. InSe tion4,weuse

thedistributed optimal ontrolof thePoisson equation andthe onve tion-diusion

equa-tion asnumeri alexperiments to illustrate the performan e of ourmethod. InSe tion 5,

weextendthispre onditioningte hniquetotheoptimal ontrol of3Dproblems. Se tion 6

drawsthe on lusion anddes ribesfuture work.

2 Problem Formulation

2.1 PDE-Constrained Optimization Problem

ConsiderthePDE- onstrainedoptimizationproblemdes ribedby

min

u, f

1

2 ku − ˆuk

2 _{+ βkfk}

2 s.t. Lu = f

in

Ω

u = u

D

on

Γ

D

,

(1)

where

L

isanoperator,

u

is thesystemstate,

f

isthe systeminput,

u

ˆ

isthedesiredstate of the system,

β

is the weight of the system input in the ost fun tion orregularization parameter and

β > 0

. In this paper, we onsider

L = −∇

2

for optimal ontrol of the

Poisson equation and

L = −ǫ∇

2 _{+ −}

→

_{w · ∇}

for optimal ontrol of the onve tion-diusion

equation. Here

−

→

w

is ave torin

Ω

,

∇

is thegradient operator,and

ǫ

is apositive s alar. Ifwewanttosolvesu haproblem numeri ally, itis learthat weneedto dis retizethese

quantities involved at some point. There are two kinds of approa hes, one is to derive

theoptimality onditionsrstandthendis retizefromthere (optimize-then-dis retize),the

otheristodis retizethe ostfun tionandthePDErstandthenoptimizethat

(dis retize-then-optimize). Fortheproblemofoptimal ontrolofthePoissonequation,bothapproa hes

leadtoequivalentsolutionswhiledierentanswersarerea hed fortheproblem ofoptimal

ontrolofthe onve tion-diusionequation[19℄. Sin eourfo usisonmultilevelsequentially

semiseparablepre onditioners,thedis retize-then-optimizeapproa his hoseninthispaper.

By introdu ing the weak formulation and dis retizing (1) using the Galerkin method,

thedis reteanalogueoftheminimizationproblem(1)istherefore,

min

u, f

1

2 u

T

M u − u

T

b + c + βf

T

M f

s.t. Ku = M f + d,

(2) where

K = [K

i,j

] ∈ R

N ×N

isthe stinessmatrix,

M = [M

i,j

] ∈ R

N ×N

_{, M}

ij

=

Z

Ω

φ

i

φ

j

dΩ

is themassmatrixand is symmetri positivedenite,

b = [b

i

] ∈ R

N

_{, b}

i

=

Z

Ω

ˆ

u

i

φ

i

dΩ

,

c ∈

(6)

R, c =

Z

Ω

ˆ

u

2 dΩ

,

d = [d

i

] ∈ R

N

_{, d}

i

= −

N

+∂N

X

j=N +1

u

j

Z

Ω

∇φ

j

· ∇φ

i

dΩ

. The

φ

i

(i = 1, 2, . . . N )

and

φ

j

(j = 1, 2, . . . N, N + 1, . . . N + ∂N )

form abasisof

V

h

0

and

V

h

g

,respe tively. Considerthe ostfun tionin(2)andasso iatewiththeequality onstrain,weintrodu e

theLagrangianfun tion

J (u, f, λ) =

1

2 u

T

_{M u − u}

T

_{b + c + βf}

T

_{M f + λ}

T

_{(Ku − Mf − d),}

where

λ

istheLagrangemultiplier. Thenitiswell-knownthattheoptimalsolutionisgiven bynding

u

,

f

and

λ

su hthat

∇

u

J (u, f, λ) = Mu − b + K

T

λ = 0,

∇

f

J (u, f, λ) = 2βMf − Mλ = 0,

∇

λ

J (u, f, λ) = Ku − Mf − d = 0.

Thisyieldsthelinearsystem





2βM

0 −M

0 M

K

T

−M

K

0 







f

u

λ





=





0 b

d





.

(3)

The system (3) is of the saddle-point system type [18℄, i.e., the system matrix, whi h is

denotedas

A

, hasthefollowingstru ture

A =

A B

T

B

0 ,

(4) where

A ∈ R

n×n

,

B ∈ R

n×m

. Forsystem(3),wehave

A =

2βM

0

0 M

and

B =

−M K

. Thesystemmatrixofthesaddle-pointsystem(3)islargeandsparse. Thusitisamenable

tosolvesu hsystemsbypre onditionedKrylovsolvers,su hasMINRES[6℄andIDR(s)[7℄.

2.2 Pre onditioning of Saddle-Point Systems

The performan e of iterative solution methods highly depends on the quality of the

pre onditioners [20℄. For numeri almethods to solve system(3)and onstru tionof

pre- onditioners,wereferto[18℄ foranextensivesurveyof numeri almethodsfor thistypeof

systems. In this paper, we study three types of pre onditioners. The rst two types

ex-ploittheMSSSstru tureoftheblo ksofthesaddle-pointsystem,whereasthese ond type

exploitstheMSSSstru ture ofthepermutedsaddle-pointsystem.

2.2.1 Blo kPre onditioners

Re allfrom(4),if

A

isnonsingular,then

A

admitsthefollowing

LDU

fa torizationgiven by





2βM

0 −M

0 M

K

T

−M

K

0 



=





I

0 I

−

1 2β

I

KM

−

1 _I









2βM

M

S









I

0 −

1 2β

I

M

−

1 _K

T

I





,

where

S = −

1 2β

M + KM

−

1 _K

T

istheS hur omplement.

The most di ult part for this fa torization is to ompute the S hur omplement

S

be ause of omputing theinverseof alargesparse matrix. Meanwhile,solvingthesystem

(7)

introdu edlater. ThentheS hur omplement

S

alsohastheMSSSstru ture. Ifweexploit theMSSS stru tureof

S

, we an both ompute

S

and solvethepre onditioned system in linear omplexity.

In this paper, we rst study two types of blo k pre onditioners for the saddle-point

system. Theyaretheblo k-diagonalpre onditioner

P

1

andtheblo klower-triangular pre- onditioner

P

2

,where

P

1 =





2β ˆ

M

ˆ

M

− ˆ

S





, P

2 =





2β ˆ

M

0 M

ˆ

−M

K

S

ˆ





,

(5)

where

M

ˆ

isanapproximationofthemassmatrixand

S

ˆ

isanapproximationoftheS hur om-plement. For

M

ˆ

and

S

ˆ

withoutapproximation,i.e.,

M = M

ˆ

and

S = S

ˆ

,thepre onditioned system

P

−

1

1 A

hasthreedistin teigenvaluesandGMRESappliedtothepre onditioned sys-tem delivers the solution in at most three steps, while the pre onditioned system

P

−

1

2 A

hastwodistin teigenvaluesandGMRESappliedtothepre onditionedsystemdeliversthe

solutionin atmosttwosteps[18℄. Forthegeneralpropertiesof

P

1

and

P

2

,wereferto[18℄ foranextensivestudy.

Aspointedoutin[19℄,toapproximatetheS hur omplement

S = −

1 2β

M + KM

−

1 _K

T

,

ˆ

S = −KM

−

1 _K

T

ould be used for big to middle range of

β

while

S = −

ˆ

1 2β

M

ould be hosenforsmall

β

. Thustheblo k-diagonalpre onditioneris

P

1 =





2β ˆ

M

ˆ

M

ˆ

KM

−

1 _K

_ˆ

T





,

(6)

forbigormiddlerangeof

β

,and

P

1 =







2β ˆ

M

ˆ

M

1 2β

M

ˆ







,

(7)

for small

β

, where

M

ˆ

and

K

ˆ

are approximated byMSSS matrix omputation. Note that thesub-blo ksof

P

1

and

P

2

allhaveanMSSSmatrixstru turesu hthat thelinearsystem

P

1 y = r

or

P

2 y = r

anbesolvedwithlinear omputational omplexity.

2.2.2 Global Pre onditioners

Sin e the blo ks of the saddle-point system(3) keep the MSSS matrixstru ture, it is

possibletopermutethesaddle-pointsystem(3)withMSSSmatrixblo kstoalinearsystem

withglobalMSSSmatrixstru ture,wherethedetailswillbeintrodu edinthenextse tion.

Thuswehavethepermutedsaddle-pointsystemdes ribedby

˜

A˜x = ˜g,

(8)

where

A

˜

,

x

˜

and

g

˜

are permutations of

A

,

f

T

_u

T

_λ

T

and

0 T

_b

T

_d

T

in (3)and

(4),respe tively. Sin etheglobalmatrix

A

˜

ofthepermutedsaddle-pointsystemisanMSSS matrix,we andoaninexa t

LU

fa torizationof

A

˜

inlinear omputational omplexitywith MSSSmatrix omputations,i.e.,

˜

A ≈ ˜

L ˜

U ,

(9)

andusethis inexa tfa torizationasapre onditioner. We allthis fa torizationin (9)the

globalpre onditioner. Sin enoinformationof

β

islostduring thepermutationand fa tor-ization,theglobalpre onditionerisindependentof

β

whileforstandardblo kpre onditioner

P

1

and

P

2

in(5)thisdoesnothold. Thisisabigadvantageoftheglobalpre onditionerover thestandardblo kpre onditioner. Numeri alexamplesin Se tion4verifythisstatement.

(8)

arable Matrix Computations

Matri es in this paperwill alwaysbereal and theirdimensions are ompatible forthe

matrix-matrixoperations and the matrix-ve toroperationswhen their sizes are not

men-tioned.

3.1 Multi-level Sequentially Semiseparable Matri es

Thegeneratorsrepresentationofsequentiallysemiseparablematri esaredenedby

Def-inition3.1[21℄.

Denition 3.1. Let

A

be an

N × N

matrix with SSS matrix stru ture and let

n

positive integers

m

1 , m

2 , · · · m

n

with

N = m

1 + m

2 + · · · + m

n

su hthat

A

anbe writtenin the following blo k-partitionedform

A

ij

=







U

i

W

i+1

· · · W

j−1

V

j

T

,

if

i < j;

D

i

,

if

i = j;

P

i

R

i−1

· · · R

j+1

Q

T

j

,

if

i > j.

wherethe supers ript

′

T

′

denotes the transposeof the matrix.

Table1: GeneratorsizefortheSSSmatrix

A

in Denition3.1

Generators

U

i

W

i

V

i

D

i

P

i

R

i

Q

i

Sizes

m

i

× k

i

k

i−1

× k

i

m

i

× k

i−1

m

i

× m

i

m

i

× l

i

l

i−1

× l

i

m

i

× l

i+1

The sequen es

{U

i

}

n−1

i=1

,

{W

i

}

n−1

i=2

,

{V

i

}

n

i=2

,

{D

i

}

n

i=1

,

{P

i

}

n

i=2

,

{R

i

}

n−1

i=2

,

{Q

i

}

n−1

i=1

are

matri eswhosesizesare listedinTable1andtheyare alled generatorsof theSSSmatrix

A

. Withthegeneratorsrepresentation,theSSSmatrix

A

isdenotedas

A = SSS(P

s

, R

s

, Q

s

, D

s

, U

s

, W

s

, V

s

).

Take

n = 5

forexample,theSSSmatrix

A

isshownby(10),

A =







D

1 U

1 V

2 T

U

1 W

2 V

3 T

U

1 W

2 W

3 V

4 T

U

1 W

2 W

3 W

4 V

5 T

P

2 Q

T

1 D

2 U

2 V

3 T

U

2 W

3 V

4 T

U

2 W

3 W

4 V

5 T

P

3 R

2 Q

T

1 P

3 Q

T

2 D

3 U

3 V

4 T

U

3 W

4 V

5 T

P

4 R

3 R

2 Q

T

1 P

4 R

3 Q

T

2 P

4 Q

T

3 D

4 U

4 V

5 T

P

5 R

4 R

3 R

2 Q

T

1 P

5 R

4 R

3 Q

T

2 P

5 R

4 Q

T

3 P

5 Q

T

4 D

5 





.

(10)

Remark 3.1. The generatorsof a SSS matrixis not unique, there existsa seriesof

non-singular transformations between two dierent sets of generators of the same SSS matrix

A

.

With the generatorsrepresentation ofSSSmatri es, basi operationsof theunderlying

matri essu hasaddition,multipli ationandinversionare losedunderSSSmatrixstru ture

and anbedoneinlinear omputational omplexity. Moreover,de omposition/fa torization

su h as

QR

fa torization [22℄ [23℄,

LU/LDU

de omposition [24℄ [14℄, and

U LV

de ompo-sition [25℄ an also be omputed in a stru ture preserving way. Many operationson SSS

matri es anbeperformedwithlinear omputational omplexity. Examplesarethe

matrix-matrixmultipli ation[21℄, thematrix-ve tormultipli ation [21℄, thematrixinversion[24℄,

the

QR

[22℄,

LU

[12℄, and

U LV

fa torization [26℄. To keep a lear stru ture of this pa-per, Table2liststhemost widelyusedoperationsfor SSSmatri esand the orresponding

(9)

operations

Ax

A ± B

AB

A

−

1 _LU

_{Lx = b}

*

referen es [13℄[24℄[21℄ [13℄[24℄ [21℄ [13℄[24℄[21℄ [23℄[27℄[28℄ [24℄[21℄ [21℄

*

L

isalower-triangularSSSmatrix.

Similarto Denition 3.1for SSSmatri es, thegeneratorsrepresentationforMSSS

ma-tri es,spe i allythe

k

-levelSSSmatri es,aredened byDenition3.2.

Denition3.2. Thematrix

A

issaidtobea

k

-levelSSSmatrixifallitsgeneratorsare

(k −

1)

-levelSSSmatri es. The

1

-level SSSmatrixistheSSSmatrixthatsatisesDenition3.1. OperationslistedinTable2fortheSSSmatri es anbeextendedtotheMSSSmatri es,

whi hyieldslinear omputational omplexityforMSSSmatri es. MSSSmatri eshavemany

appli ations,oneofthem isthedis retizedpartial dierentialequations(PDEs) [15℄[14℄.

Example 3.1. For the

P

1

nite-elementdis retization of the 2D Lapla ian equation with homogeneousDiri hlet boundary ondition, the stiness matrix

K

isgiven by

K =







A

B

A

B

. . . . . . . . . . . .

B

A







,where

A =







4 −1

−1

4 −1

−1

. . . . . . . . . . . .

−1

4 





,

B = −I

,and

I

is

the identitymatrix. The matrix

A

and

B

areboth SSSmatri es bedenoted by

A =

SSS(1, 0, −1, 4, 1, 0, −1),

B

=

SSS(0, 0, 0, −1, 0, 0, 0).

Thematrix

K

has the MSSS(2-level SSS)matrixstru tureandisdenotedby

K = MSSS(I, 0, B

T

_{, A, I, 0, B}

T

_).

Remark3.2. Similar withSSS matri es, for MSSSmatrix, itsgeneratorsarenot unique.

There exists a set of nonsingular transformations between two dierent sets of generators

for aspe iedMSSSmatrix.

Remark 3.3. For SSS or MSSS matri es, it is not ne essary for their diagonals,

sub-diagonalorsuper-diagonals tobe onstantlikethatinExample 3.1. Theirsizes anevenbe

dierentaslongasthe blo k-partitionedrepresentation in Denition3.1issatised.

Note that for a saddle-point system from the PDE- onstrainedoptimization problem,

allits blo ks areMSSS matri es, whi h enablesus to omputethe

LU

fa torization ofall its blo ks with MSSS matrix omputations in linear omputational omplexity. However,

we fail to omputethe

LU

fa torization ofthe whole saddle-point systemmatrix be ause thesaddle-point systemmatrix is notan MSSSmatrix but just has MSSSmatrix blo ks.

Thefollowinglemmatellsus how topermuteamatrixwithSSS matrixblo ksto asingle

SSSmatrix. We aneasilyextendthislemmatotheMSSSmatrix ases,whi hallowsusto

permuteamatrixwithMSSSmatrixblo ksto asingleMSSSmatrix.

Lemma3.1. [29℄Let

A

,

B

,

C

and

D

beSSS matri eswiththe generatorsrepresentations

A =

SSS(P

s

a

, R

a

s

, Q

a

s

, D

s

a

, U

s

a

, W

s

a

, V

s

a

),

B

=

SSS(P

b

s

, R

b

s

, Q

b

s

, D

s

b

, U

s

b

, W

s

b

, V

s

b

),

C

=

SSS(P

c

s

, R

c

s

, Q

c

s

, D

s

c

, U

s

c

, W

s

c

, V

s

c

),

D

=

SSS(P

d

s

, R

d

s

, Q

d

s

, D

d

s

, U

s

d

, W

s

d

, V

s

d

).

(10)

f

g

=

A B

C

D

a

b

,

and

f

g

= T

a

_b

areequivalentwithrowand olumnpermutationsofthematrixblo ks. Theve tors

f

g

and

a

b

arepermutationsof

f

g

and

a

b

,respe tively. Thematrix

T

isanSSSmatrixandhas the generatorsrepresentation

T = SSS(P

t

s

, R

t

s

, Q

t

s

, D

s

t

, U

s

t

, W

s

t

, V

s

t

),

where

P

t

s

=

P

a

s

P

s

b

0

0 P

c

s

P

s

d

, R

t

s

=







R

a

s

R

b

s

R

c

s

R

d

s







, Q

t

s

=

Q

a

s

0 Q

c

s

0

0 Q

b

s

0 Q

d

s

, D

s

t

=

D

a

s

D

s

b

D

c

s

D

d

s

, U

t

s

=

U

a

s

U

s

b

0

0 U

c

s

U

s

d

, W

t

s

=







W

a

s

W

b

s

W

c

s

W

d

s







, V

t

s

=

V

a

s

0 V

s

c

0

0 V

b

s

0 V

s

d

.

Remark 3.4. Lemma 3.1 is for a

2 × 2

blo k matrix, but it an be extended to matri es withdierent numberofblo ksaswell.

Remark 3.5. Extending Lemma 3.1 tothe

k

-level SSS matrix ase isalso possible. If

A

,

B

,

C

,and

D

are

k

-level SSSmatri es, thentheir generatorsare

(k − 1)

-levelSSS matri es. Forthe permuted

k

-levelSSSmatrix

T

,its

(k − 1)

-levelSSSmatrixgeneratorswith

(k − 1)

-levelSSS matrixblo ksarederivedfromthe permutationsof rowsandblo kstogetasingle

(k − 1)

-level SSS matrixbyLemma3.1.

Forthesaddle-pointsystem(3)derivedfromthe2DPDE- onstrainedoptimization

prob-lem,dis retizingusing

P

1

niteelementmethodyieldsasaddle-pointsystemthathasMSSS (

2

-levelSSS)matrixblo ks. Thestru ture ofthesaddle-pointsystemmatrixformesh size

h = 2

−

3

is shown in Figure 1(a). Permuting the saddle-point system using Lemma 3.1

givessystem (8). The saddle-point systemmatrix stru ture before and after permutation

areshownin Figure1.

0

50

100

150

0

20

40

60

80

100

120

140

160

180 nz = 2120

(a)Beforepermutation.

0

50

100

150

0

20

40

60

80

100

120

140

160

180 nz = 2120

(b)Afterpermutation.

Figure 1: Stru tureofsystemmatrixof(3)beforeand afterpermutationfor

h = 2

−

3

.

3.2 Multi-level Sequentially Semiseparable Pre onditioners

TheabilitytosolvealinearsystemwithMSSSmatrixstru tureinlinear omputational

(11)

fa -rst introdu e the

LU

fa torization of MSSS matri esand then givea novel model order redu tionalgorithm for SSS matri es that is requiredin omputing the

LU

fa torization. For omparison,the onventionalmodelorderredu tionalgorithmisalsodis ussed.

3.2.1

LU

Fa torization ofMultilevelSequentiallySemiseparable Matri es Thesemiseparable orderdened inDenition 3.3 playsan importantrule in theMSSS

matrix omputations. Note that Dewilde et. al. and Golberg et. al. studied this kind

of stru tured matri es separately, SSS matri es named in [21℄ are alled quasiseparable

matri esin[24℄. HereweusetheMATLABstyleofnotationformatri es,i.e.,foramatrix

A

,

A(i : j, s : t)

sele tsrowsofblo ksfrom

i

to

j

and olumnsofblo ksfrom

s

to

t

of

A

.

Denition3.3. [16℄Let

rank

A(k + 1 : n, 1 : k) = l

k

, k = 1, 2, · · · , n − 1.

The numbers

l

k

(k = 1, 2, · · · , n − 1)

are alled the lower order numbers of the matrix

A

. Let

rank

A(1 : k, k + 1 : n) = u

k

, k = 1, 2, · · · , n − 1.

The numbers

u

k

(k = 1, 2, · · · , n − 1)

are alled the upper ordernumbers of the matrix

A

. Set

r

l

_{= max l}

k

and

r

u

_{= max u}

k

,where

r

l

and

r

u

are alledthelowerquasi-separableorder

andthe upper quasi-separable orderof

A

,respe tively.

Denition 3.4. [30℄ The SSS matrix

A

with lower and upper semiseparable order

r

l

and

r

u

is alledblo k

(r

l

_{, r}

u

₎

semiseparable.

Thedenitionsin Denition3.3and3.4ofSSSmatri es anbedire tlyextendedtothe

MSSSmatri es,whi hleadstoDenition3.5and 3.6.

Denition3.5. Letthematrix

A

bean

N × N

blo k

k

-levelSSS matrixwithitsgenerators be

M × M

blo k

(k − 1)

-level SSS matri es. Let

rank

A(k + 1 : N, 1 : k) = l

k

, k = 1, 2, · · · , N − 1.

Thenumbers

l

k

(k = 1, 2, · · · , N −1)

are alledthe

k

-levellowerordernumbersofthematrix

A

. Let

rank

A(1 : k, k + 1 : N ) = u

k

, k = 1, 2, · · · , N − 1.

The numbers

u

k

(k = 1, 2, · · · , N − 1)

are alled the

k

-level upper order numbers of the matrix

A

. Set

r

l

_{= max l}

k

and

r

u

_{= max u}

k

, where

r

l

and

r

u

are alled the

k

-level lower semiseparableorderandthe

k

-level uppersemiseparable order of the

k

-level SSS matrix

A

, respe tively.

Denition3.6. The

k

-levelSSSmatrix

A

with

k

-level loweranduppersemiseparableorder

r

l

and

r

u

is alled

k

-level blo k

(r

l

_{, r}

u

₎

semiseparable.

Withthesedenitions,wehavethefollowingalgorithmto omputethe

LU

fa torization ofa

k

-levelSSSmatrix.

Lemma3.2. [12℄[14℄Let

A

beastronglyregular

N ×N

blo k

k

-levelsequentially semisepara-blematrixof

k

-levelblo k

(r

l

_{, r}

u

₎

semiseparableanddenotedbyitsgeneratorsrepresentation

A = MSSS(P

s

, R

s

, Q

s

, D

s

, U

s

, W

s

, V

s

)

. Let

A = LU

be its blo k LU fa torization. Then,

1. The fa tor

L

is a

k

-level sequentially semiseparable matrix of

k

-level blo k

(r

L

_{, 0)}

semiseparable and

U

is a

k

-level sequentially semiseparable matrix of

k

-level blo k

(0, r

U

₎

semiseparable. Moreover,

r

L

_{= r}

l

and

r

U

_{= r}

u

.

(12)

2. The fa tors

L

and

U

anbe denotedby the generatorsrepresentation

L

= MSSS(P

s

, R

s

, ˆ

Q

s

, D

L

s

, 0, 0, 0),

U

= MSSS(0, 0, 0, D

U

s

, ˆ

U

s

, W

s

, V

s

).

where

Q

ˆ

s

,

D

L

s

,

D

U

s

and

U

ˆ

s

are

(k − 1)

-level sequentially semiseparable matri es and omputedbythe followingalgorithm:

Algorithm1

LU

fa torizationofa

k

-levelSSSmatrix

A

Initialize:

M

1 ← 0 ∈ R

r

l

×

_r

u

bea

(k − 1)

-levelSSSmatrix Computethe

LU

fa torizationofthe

(k − 1)

-levelSSSmatrix

D

1 = D

L

1 D

U

1

,let

U

ˆ

1 = (D

L

1 )

−

1 _U

1

and

Q

ˆ

1 = (D

L

1 )

−

T

_Q

1

for

i = 2 : N − 1

do

M

i

= ˆ

Q

T

i−1

U

ˆ

i−1

+ R

i

M

i−1

W

i

,

Computethe

LU

fa torizationofthe

(k − 1)

-levelSSSmatrix

(D

i

− P

i

M

i

V

i

) = D

L

i

D

U

i

, Let,

U

ˆ

i

= (D

L

i

)

−

1 _(U

i

− P

i

M

i−1

W

i

)

,

Q

ˆ

i

= (D

U

i

)

−

T

_(Q

i

− V

i

M

i−1

T

R

T

i

)

. endfor

Computethe

LU

fa torizationofthe

(k − 1)

-levelSSSmatrix

D

n

− P

n

M

n−1

V

n

T

= D

L

n

D

n

U

Output:

D

L

i

,

D

U

i

,

Q

ˆ

i

,

U

ˆ

i

Proof. Fortheproofofthelemma,wereferto[12℄and[14℄.

Remark3.6. InAlgorithm1,the

LU

fa torization of a

0

-level SSS matrixisjustthe

LU

fa torizationof anordinarymatrixwithoutSSS stru ture.

ForMSSSmatri es, matrix-matrixoperationssu h asaddition and multipli ation will

leadtoagrowthofthe semiseparableorder,whi h anbeveriedfrom thematrix-matrix

operationsintrodu edin[21℄[24℄. Thisresultsinthegrowthofthe omputational

omplex-ity. Takethe1-levelSSSmatrix

A

forexample,theopsneededfor omputing

A

2

is

40n

3 _N

where

n

isthesemiseparableorder[21℄and

N

isthenumberofblo ksof

A

. Tobespe i , thefollowinglemmaisintrodu ed.

Lemma3.3. [24℄Let

A

1

,

A

2

beSSSmatri esofsizes

N ×N

whi harelowersemiseparable of orders

m

1

,

n

1

respe tively. Then the produ t

A

1 A

2

is lower semiseparable of order at most

m

1 + n

1

. Let

A

1

,

A

2

beSSS matri esof sizes

N × N

whi hareuppersemiseparableof orders

m

2

,

n

2

respe tively. Thenthe produ t

A

1 A

2

isupper semiseparable of orderatmost

m

2 + n

2

.

Remark3.7. For

k

-level SSSmatri es,sin esemiseparableordervaries atdierentlevels, result of Lemma 3.3 holds for the

k

-level semiseparable order. But we do not know the

(k − 1)

-level semiseparable order of the

(k − 1)

-level SSS generators exa tly, wejust know the

(k − 1)

-level semiseparable orderalsoin reases.

Lemma3.3givesrisetothequestionwhetherthereexistsaminimalsemiseparableorder

foraSSSmatrixsu h that theSSSmatrixwithabiggersemiseparableorder isequivalent

toanSSSmatrixwithminimalsemiseparableorder. Denition3.7andLemma3.4givethe

answerto theaforementionedquestion.

Denition3.7. [16℄Wesaythatthe lowergenerators

P

i

(i = 2, . . . , N )

,

Q

j

(j = 1, . . . , N −

1)

,

R

k

(k = 2, . . . , N −1)

ofanSSSmatrix

A

areminimalifalltheirorders

l

k

(k = 1, . . . , N −

1)

areassmallaspossibleamongall lower generatorsofthe samematrix

A

,i.e., for lower generatorsof the matrix

A

with orders

l

′

(13)

l

k

≤ l

′

k

, k = 1, . . . , N − 1

hold. We also saythat the orders

l

k

(l = 1, . . . , N − 1)

are the minimal ordersof the lower generatorsof

A

.

Lemma 3.4. [16℄Let

A = {A

ij

}

N

i,j=1

be a blo k matrix with lower rank numbers

r

k

(k =

1, . . . , N − 1)

. Then

A

has lower generators with orders equal to the orresponding rank numbers. Moreover, for any matri es, the rank numbers are the minimal orders of the

generators.

Remark3.8. Lemma3.4 an beextendedtothe

k

-level SSSmatri es dire tly.

Remark3.9. Lemma3.4showsthatthereexistsaminimalsemiseparableorderforanSSS

matrix. Thus, for an SSS matrixof semiseparableorderbigger than the minimal separable

order,the semiseparableorder anberedu edtomake theredu edsemiseparable orderequal

toorsmaller thanthe minimalsemiseparable ordersu hthatthe resultingSSSmatrixwith

redu ed semiseparable order is equal to or equivalent with the SSS matri es without order

redu tionuptoasmall toleran e.

Theaimofmodelorderredu tionofa

k

-levelSSSmatrix

A

withitsgenerators represen-tation

A = MSSS(P

s

, R

s

, Q

s

, D

s

, U

s

, W

s

, V

s

)

istond

(k−1)

-levelSSSmatri es

P

ˆ

s

,

R

ˆ

s

,

ˆ

Q

s

,

U

ˆ

s

,

W

ˆ

s

,

V

ˆ

s

ofsmallersize omparedwith

P

s

,

R

s

,

Q

s

,

U

s

,

W

s

,

V

s

,respe tivelysu hthat

ˆ

A = MSSS( ˆ

P

s

, ˆ

R

s

, ˆ

Q

s

, D

s

, ˆ

U

s

, ˆ

W

s

, ˆ

V

s

)

isof

k

-levelsemiseparableordersmallerthanor equalto theminimal

k

-levelsemiseparable orderof

A

. Meanwhile,

A

ˆ

isanapproximation of

A

uptoasmalltoleran e

ǫ

, i.e.,

k ˆ

A − Ak < ǫ

.

Remark 3.10. In Algorithm 1, for omputing the

LU

fa torization of a

k

-level SSS ma-trix, matrix-matrix operations are performedon its generatorswhi h are

(k − 1)

-level SSS matri es. Su hoperations lead tothegrowthof semiseparable orderofthe

(k − 1)

-level SSS matri es, whi h indu esgrowthof omputational omplexity. Model order redu tion is

ne -essary toredu ethe semiseparable orderor keep the semiseparable order under athreshold

duringthe

LU

fa torization, su has omputationof the re urren eof

M

i

inAlgorithm1. Remark 3.11. Sin e the

LU

fa torization of a

k

-level SSS matrixneeds the model order redu tionfor

(k − 1)

-levelSSSmatri es, the

LU

fa torizationinLemma3.2isanexa t fa -torizationforSSSmatri esbe ausenomodelorderredu tionisneededforordinarymatri es

(

0

-level SSS matri es). Itis an inexa t fa torizationfor the

k

-level

(k ≥ 2)

SSS matri es. Therefore, for dis retizedone-dimensional PDEs onaregular grid,this fa torization ould

beperformedasadire t solverandasan e ientpre onditionerfor the dis retizedtwo-or

higher-dimensional PDEsonaregulargrid.

Remark 3.12. The model orderredu tion algorithm for SSS matri es has been studiedin

[13℄[16℄,whilefor

2 −level

orevenhigher-levelSSSmatri es, itisstillabig hallengesin e modelorderredu tionof

k

-levelSSSmatri eswhere

k ≥ 2

needsthe redu edgeneratorsstill be

(k − 1)

-level SSS matri es. The model order redu tion algorithms in [13℄ [16℄ applied to the

k

-level SSS matri es will not return stru ture preserving

(k − 1)

-level SSS matrix generators.

3.2.2 Approximate Balan ed Trun ation

In this paper, we design a novel model order redu tion algorithm for SSS matri es.

Withthisalgorithm,we an onstru tane ientpre onditionerfortwo-dimensional

PDE- onstrained optimization problem, whi h will be studied in the next se tion. The

orre-sponden e betweenSSSmatri es and thelinear time-varying (LTV)systemsmotivates us

(14)

Dewildeet. al. in [27℄. Consideramixed- ausalsystemthat isdes ribed bythefollowing state-spa emodel

x

c

i+1

x

a

i−1

=

R

i

W

i

x

c

i

x

a

i

+

Q

i

V

i

u

i

y

i

=

P

i

U

i

x

c

i

x

a

i

+ D

i

u

i

,

(11) where

x

c

denotes the ausal systemstates,

x

a

representsthe anti- ausal systemstates,

u

i

isthesysteminput, and

y

i

isthesystemoutput. Withzeroinitialsystemstatesandsta k all theinput and output as

u = u

¯

T

1 , u

T

2 , . . .

u

T

N

T

,

y = y

¯

T

1 , y

2 T

, . . .

y

N

T

, the

matrix

H

thatdes ribestheinput-outputbehaviorofthismixed- ausalsystem,i.e.,

y = Hu

, indu esanSSSmatrixstru ture. Take,

N = 4

forexample,thematrix

H

is,

H =







D

i

U

1 V

2 U

1 W

2 V

3 U

1 W

2 W

3 V

4 P

2 Q

1 D

2 U

2 V

3 U

2 W

3 V

4 P

3 R

2 Q

1 P

3 Q

2 D

3 U

3 V

4 P

4 R

3 R

2 Q

1 P

4 R

3 Q

2 P

4 Q

3 D

4 





.

(12)

Remark3.13. Toredu e thesemiseparable orderof the SSS matrix

H

in (12),the orders of

P

s

,

R

s

,

Q

s

,

U

s

,

W

s

and

V

s

needtoberedu ed. This orrespondstoredu etheorderofthe mixed- ausalLTVsystem(11). Model redu tion forLTVsystem(11) ouldbeperformedto

redu ethe semiseparable orderof

H

.

Model orderredu tion forLTV systemsis studied in [31℄ [32℄. In [32℄, alinear matrix

inequality(LMI)wasintrodu edtosolvetheLyapunovinequalitiesforthe ontrollabilityand

observabilityGramians. In[31℄,thelow-rankSmithmethodwaspresentedto approximate

thesquare-rootofthe ontrollabilityandobservabilityGramians.

Sin ethe ausalLTVsystemandtheanti- ausalLTVsystemhavesimilarsystem

stru -turethat orrespondto thestri tly lower-triangularpartand thestri tlyupper-triangular

partof thematrix

H

,respe tively. Herewejust onsiderthe ausalLTVsystemdes ribed bythefollowingstate-spa emodel,

(

x

k+1

= R

k

x

k

+ Q

k

u

k

y

k

= P

k

x

k

,

(13)

overthetimeinterval

[k

o

, k

f

]

withzeroinitialstates.The ontrollabilityGramian

G

c

(k)

and observabilityGramian

G

o

(k)

are omputedfromthefollowingSteinre urren eformulas:

G

c

(k + 1) = R

k

G

c

(k)R

T

k

+ Q

k

Q

T

k

,

(14)

G

o

(k)

= R

T

k

G

o

(k + 1)R

k

+ P

k

T

P

k

,

(15)

withinitial onditions

G

c

(k

o

) = 0

and

G

o

(k

f

+ 1) = 0

.

Note that the ontrollability Gramian

G

c

(k)

andobservability Gramian

G

o

(k)

are pos-itive deniteif thesystemis ompletely ontrollable andobservableorsemi-denite ifthe

systemispartly ontrollableandobservable,thustheireigenvaluesarenon-negative. Their

eigenvaluesoften havealarge jump at anearly stage as pointedout in [17℄ [33℄ [34℄ [35℄,

whi h suggeststo approximatethese twoGramiansatea hstepbyalow-rank

approxima-tion. Belowweshowhowtoobtainsu happroximations. Sin ethe ontrollabilityGramian

G

c

(k)

and observability Gramian

G

o

(k)

have similar stru ture, we will only fo us on the ontrollabilityGramian

G

c

(k)

.

Thekeypointofthelow-rankapproximationistosubstitutetheCholeskyfa torization

ofthe ontrollabilityGramian

G

c

(k)

(15)

where

L

k

∈ R

N ×N

inea hstep

k

byitsapproximateCholeskyfa torization,

˜

G

c

(k) = ˜

L

k

L

˜

T

k

,

(17)

with

L

˜

k

∈ R

N ×n

k

where

n

k

is the numeri al rank of

G

c

(k)

and

N > n

k

at ea h step

k

. Typi ally,

n

k

issettobe onstant,i.e.,

n

k

= n

atea hstep. Sin e

G

c

(k)

isoflownumeri al rank,itisreasonabletousetherank

n

k

fa tor

L

˜

k

toapproximate

G

c

(k)

.

In[17℄,are ursivelow-rankGramianmethodwasintrodu edtoapproximatethe

Grami-ansofalineartime-invariantsystem. Here,weextendthatmethodtothelineartime-varying

systems,whi h is similar withthe method in [36℄. This method is shown in Algorithm2.

From [17℄, weknow that

G

˜

c

(i) = ˜

L

c

(i) ˜

L

c

(i)

T

and

G

˜

o

(i) = ˜

L

o

(i) ˜

L

o

(i)

T

in Algorithm 2are

thebest rank

n

approximationsto

G

c

(i)

and

G

o

(i)

. Algorithm2Low-rankapproximationoftheGramians

Initialize:

G

˜

c

(1) ← 0 ∈ R

M ×n

,

G

˜

o

(N + 1) ← 0 ∈ R

M ×n

,

N

isthenumberoftimesteps,

M

istheunredu ed order,

n

isthenumeri alrank.

fori=2: Ndo

Computethesingularvaluede ompositions

Q

i−1

R

i−1

G

˜

c

(i − 1)

= U

c

Σ

c

V

c

T

,

P

T

i

R

T

i

G

˜

o

(i + 1)

= U

o

Σ

o

V

o

T

. Let

U

c

=

U

c1

U

c2

,

Σ

c

=

Σ

c1

Σ

c2

with

U

c1

∈ R

M ×n

and

Σ

c1

∈ R

n×n

.

U

o

=

U

o1

U

o2

,

Σ

o

=

Σ

o1

Σ

o2

with

U

o1

∈ R

M ×n

and

Σ

o1

∈ R

n×n

. Make

˜

L

c

(i) = U

c1

Σ

c1

,

L

˜

o

(i) = U

o1

Σ

o1

. endfor Output:

L

˜

c

(i) ∈ R

M ×n

and

L

˜

o

(i) ∈ R

M ×n

.

With the approximate ontrollability Gramian

G

c

(i)

and observability Gramian

G

o

(i)

, the balan ed trun ation ould beperformed to redu e the order of theLTV system. For

theapproximate balan ed trun ation, thekeyis to usethe low-rankapproximationofthe

fa torsof Gramianstoprovideanapproximationtothebalan edtrun ation.

FortheLTVsystem(13),todoabalan edtrun ation,rstthesystemstatesare

trans-formedbythenonsingulartransformation

x

k

= T

k

x

¯

k

togeta"balan ed"system,

( ¯x

k+1

= T

−

1 k+1

R

k

T

k

x

k

+ T

−

1 k+1

Q

k

u

k

y

k

= P

k

T

k

x

k

,

(18)

wherethestates

x

¯

k

= ˜

x

T

k

x

ˆ

T

k

T

. Thekeptsystemstatesare

x

˜

k

=

I

n

0 ¯

x

k

where

n

is thesystemorderafterredu tion. Theredu edLTVsystemof(13)is

(

˜

x

k+1

= Π

l

(k + 1)R

k

Π

r

(k)˜

x

k

+ Π

l

(k + 1)Q

k

u

k

y

k

= P

k

Π

r

(k)˜

x

k

,

(19) where

Π

l

(k + 1) =

I

n

0 T

−

1 k+1

and

Π

r

(k) = T

k

I

n

0

.

Next,weextendthebalan edtrun ationalgorithmtothelineartime-varying ase. This

methodisdes ribedin Algorithm3.

Remark3.14. These ondloopofAlgorithm3ensuresthat

Π

l

(i)

and

Π

r

(i)

are"balan ed". This is vital sin e we approximate the ontrollability andobservability Gramians

(16)

(19) bythe low-rankapproximatebalan edtrun ation.

Remark 3.16. For an SSS matrix

A = SSS(P

s

, R

s

, Q

s

, D

s

, U

s

, W

s

, V

s

)

with lower semiseparableorder

M

,Algorithm2andAlgorithm3 ouldbeperformedtothestri tly lower-triangular partof

A

toredu e the lower semiseparable orderto

n

,yielding the approximate SSS matrix

A = SSS( ˜

˜

P

s

, ˜

R

s

, ˜

Q

s

, D

s

, U

s

, W

s

, V

s

)

. For the stri tly upper-triangular part of

A

, rst transpose it to be stri tly lower-triangular then perform Algorithm 2 and Algorithm3. After the redu tion,transpose the redu edstri tly lower-triangular part tobe

stri tlyupper-triangular.

Algorithm3Approximate balan edtrun ationforLTVsystems

Pro edure: Setthenumeri alrank

n

.

Usethe low-rankapproximationAlgorithm 2to omputethe rank

n

approximationsto the ontrollabilityGramian

G

c

(i)

andobservabilityGramian

G

o

(i)

,denotedby

G

˜

c

(i)

and

˜

G

o

(i)

,respe tively. loop

Computethesingularvaluede omposition

˜

G

T

c

(i) ˜

G

o

(i) = U

i

Σ

i

V

i

T

. endloop loop Let

Π

l

(i) = ˜

G

o

(i)V

i

Σ

−

1

2 i

,

Π

r

(i) = ˜

G

c

(i)U

i

Σ

−

1

2 i

. endloop End Pro edure Output:

Π

l

(i) ∈ R

n×M

and

Π

r

(i) ∈ R

M ×n

.

3.2.3 Hankel Blo ksApproximation

Themodelorderredu tionalgorithmsforSSSmatri esin [13℄[21℄[27℄approximatethe

Hankelblo ksoftheSSSmatri es,wheretheHankelblo ksofanSSSmatrix

A

aredened byDenition3.8.

Denition 3.8. [13℄Hankel blo ks denote the o-diagonal blo ks that extendfromthe

di-agonaltothe northeast orner(forthe upper ase)ortothe southwest orner (forthe lower

ase).

Takea

4×4

SSSmatrix

A

forexample,theHankelblo ksforthestri tlyuppertriangular partareshownin Figure2by

H

1

,

H

2

and

H

3

.

1,2

A

1,3

A

1,4

2,3

A

2,4

3,4

A

1 H

2 H

3 H

Figure 2: Hankelblo ksofaSSSmatrix

A

Themodelorder redu tionalgorithmsin [13℄[21℄[27℄ are

Hankel norm

optimalorder redu tion[29℄algorithms. Thatis,givenanSSSmatrix

A

withalowersemiseparableorder

r

L

andanuppersemiseparableorder

r

U

,we angetanapproximate SSSmatrix

A

ˆ

with a lowersemiseparableorder

˜

r

L

andanuppersemiseparableorder

r

˜

U

where

r

L

> ˜

r

L

,

r

U

> ˜

r

U

a hieves

(17)

where

kAk

H

= max

i

kH

i

(A)k

2

and

H

i

(A)

aretheHankelblo ksof

A

denedinDenition3.8. For omparison, this model order redu tionalgorithm to the stri tly upper-triangular

partofSSSmatri esislistedinAlgorithm4[13℄.

Algorithm4Hankelblo ksapproximationforSSSmatri es

Initialize:

H ← 0 ∈ R

M ×M

,

G ← 0 ∈ R

M ×M

,

M

is theuppersemiseparableorderbefore redu tion,set theredu eduppersemiseparableorder

m

andthenumberofblo ks

N

. performtheforwardre ursion

for

i = 1 : N − 1

do

Computethesingularvaluede omposition(SVD)

H

U

i

= U ΣV

T

and partition

U =

U

T

U

K

=

U

a

(·)

(·) (·)

with

U

K

∈ R

n

u

×

(·)

,

U

a

∈

R

M ×M

,where

n

u

isthenumberofrowsof

U

i

. Let,

U

i

= U

K

,

W

i

= U

a

,

V

i+1

= ΣV

T

_V

i+1

and

H = ΣV

T

_W

i+1

. endfor

performtheba kwardre ursion

for

i = N : −1 : 2

do

Computethesingularvaluede omposition(SVD)

V

i

G

T

= U ΣV

T

and partition

U =

U

a

(·)

U

b

(·)

with

U

a

∈ R

n

v

×

_M

,

U

b

∈ R

M ×M

,

Σ =

Σ

a

(·)

with

Σ

a

∈ R

M ×M

,where

n

v

isthenumberofrowsof

V

i

. Let,

V

i

= U

a

,

U

i−1

= U

i−1

V Σ

T

a

,

W

i

= U

T

b

,

G = W

i−1

V Σ

T

a

. endfor dothetrun ation for

i = 1 : N

do Partition

U

i

=

U

i1

(·)

with

U

i1

∈ R

(·)×m

,

W

i

=

W

i1

(·)

with

W

i1

∈ R

m×m

,

V

i

=

V

i1

(·)

with

V

i1

∈ R

m×(·)

. Let,

˜

U

i

= U

i1

,

W

˜

i

= W

i1

,

V

˜

i

= V

i1

. endfor Output:

U

˜

i

∈ R

(·)×m

,

W

˜

i

∈ R

m×m

and

V

˜

i

∈ R

m×(·)

.

Remark 3.17. With the Hankel blo ks approximation, we an also onstru t an e ient

pre onditioner for two-dimensional PDE- onstrained optimization problem, whi h will be

studiedin thenext se tion.

Remark3.18. Foran SSS matrix

A

withlower and upper semiseparable order

r

l

and

r

u

, respe tively. The bigger the semiseparable order

ˆ

r

l

and

ˆ

r

u

after model order redu tion by Algorithm2-3or4is,the losertheredu edSSSmatrix

A

ˆ

isto

A

.Forapropersemiseparable order set, the model order redu tion is a urate enough. This makes the

LU

fa torization of the

2

-level SSS matrixbyAlgorithm1 a urate enough that an beperformedasadire t solver. Numeri al experimentsin thenext se tion illustratethis.

GivenanSSSmatrix

A = SSS(P

S

, R

s

, Q

s

, D

s

, U

s

, W

s

, V

s

)

,to omparetheopsof theapproximatebalan edtrun ationinAlgorithm2-3andtheHankelblo ksapproximation

Algorithm4,weassumethat thegeneratorssizesin Table3.1are

m

i

= n

and

k

i

= l

i

= M

where

N

isthe number of SSSblo ks and

N ≫ M ≫ n

. This is easy to verify from the

(18)

SSSmatrix

A = SSS( ˜

˜

P

S

, ˜

R

s

, ˜

Q

s

, D

s

, ˜

U

s

, ˜

W

s

, ˜

V

s

)

,where

k

˜

i

= ˜

l

i

= m

,

m

istheredu ed semiseparableorderand

m ≪ M

. ForAlgorithm2-3,theops ount

F

N

are

F

N

= O (3m

2 + 4mn + n

2 )M N + (m + n)M

2 N

,

(20)

whiletheops ount

F

C

forAlgorithm4is

F

C

= O M

3 N + (2m + n)M

2 N + 2mnM N .

(21)

Sin e

N ≫ M ≫ m, n

,wedes ribethat

F

N

= O M

2 N

,

(22)

and

F

C

= O M

3 N .

(23)

Remark3.19. From(22)and(23),itisobviousthatbothmodelorderredu tionalgorithm

forSSSmatri eshavelinear omputational omplexity

O(N)

,whiletheapproximatebalan ed trun ation (ops ount denoted by

F

N

) is omputationally heaper thanthe Hankel blo ks approximation (ops ountdenotedby

F

C

)for largeenough

M

. Thiswillalsobeillustrated bynumeri al experiments inthe next se tion.

Remark3.20. Asstatedin[32℄,thebalan edtrun ationyieldsanoptimalindu ed

L

2

-norm approximation. Thus forthe approximatebalan edtrun ation,theredu ed ontrol systemis

losetotheoptimal

L

2

-normapproximation. TheAlgorithm4returnsthe optimal

Hankel

-normapproximation. Thus,both algorithmsfor model orderredu tion ofSSS matri eswill

yield an a urate approximation. Sin e the inequality

kZk

H

6 kZk

2

6 √

N kZk

H

for all

Z ∈ R

n×n

holds [29℄, the Hankel blo ks approximation Algorithm4 yields amorea urate

approximation thanthe approximate balan ed trun ation Algorithm 2-3 in theory. But the

a ura yofAlgorithm4andAlgorithm2-3are omparable,whi hwillbeshownbynumeri al

experimentsinthe next se tion.

4 Numeri al Experiments

We study twotest examplesfor optimal ontrol of2D PDEs in this se tion, i.e.,

opti-mal ontrol ofthe onve tion-diusion equationin Example4.1and optimal ontrolofthe

Poisson equation in Example A.1 in the appendix. We apply the blo k-diagonal

pre on-ditioner

P

1

in (5) forthe MINRES method and thelower-triangularpre onditioner

P

2

in (5) for the IDR(s) method to both examples. The global pre onditioner

A

ˆ

in (9) is also performed for the two test examplesto show its superior performan e over the standard

blo kpre onditionersforsaddle-pointsystem.

Example4.1. [19℄Let

Ω = {(x, y)|0 ≤ x ≤ 1, 0 ≤ y ≤ 1}

and onsider the problem

min

u,f

1

2 ku − ˆuk +

β

2 kfk

2 s.t. − ǫ∇

2 _{u + −}

→

_{ω .∇u = f}

in

Ω

u = u

D

on

Γ

D

,

where

Γ

D

= ∂Ω

and

u

D

=

(2x − 1)

2 _{(2y − 1)}

2

if

0 ≤ x ≤

1

2 ,

and

0 ≤ y ≤

1

2 ,

0

otherwise

.

ǫ

isa positive s alar,

−

→

ω

is the unit dire tional ve tor that

−

→

ω = (cos(θ), sin(θ))

T

and the

(19)

of 2.53GHzand 4Gb memorywith Matlab R2010b. Thestop toleran eof the 2-normof

therelativeresidualisset tobe

10 −

6

forallthenumeri alexperiments. Theproblemsizes

3.07e+03, 1.23e+04, 4.92e+04 and 1.97e+05 orrespond to the mesh sizes

h = 2

−

5

,

2 −

6

,

2 −

7

, and

2 −

8

,respe tively. Themaximumsemiseparableorder isin thebra ketsfollowing

theproblemsize. Thetimeto omputethepre onditioners anditerativesolutionmethods

timeismeasuredin se onds.

4.1 Comparison of Two Model Order Redu tion Algorithms

Inthispart,wetesttheperforman eofthetwomodelorderredu tionalgorithms.

Con-siderthe pre onditioningofoptimal ontrol ofthe onve tion-diusionequation des ribed

inExample4.1. Withtheblo k-diagonalpre onditioner

P

1

byapproximatebalan ed trun- ation and theHankel blo ks approximation methods, the results for dierentvalues of