Noname manuscript No.

(will be inserted by the editor)

Multilevel quasiseparable matrices in PDE-constrained
optimization
Jacek Gondzio Pavel Zhlobich
the date of receipt and acceptance should be inserted later
Abstract Optimization problems with constraints in the form of a partial
dierential equation arise frequently in the process of engineering design. The
discretization of PDE-constrained optimization problems results in large-scale
linear systems of saddle-point type. In this paper we propose and develop a
novel approach to solving such systems by exploiting so-called quasiseparable
matrices. One may think of a usual quasiseparable matrix as of a discrete
analog of the Greens function of a one-dimensional dierential operator. Nice
feature of such matrices is that almost every algorithm which employs them
has linear complexity. We extend the application of quasiseparable matrices
to problems in higher dimensions. Namely, we construct a class of precon-
ditioners which can be computed and applied at a linear computational cost.
Their use with appropriate Krylov methods leads to algorithms of nearly linear
complexity.
Keywords saddle-point problems, PDE-constrained optimization, precondi-
tioning, optimal control, linear systems, quasiseparable matrices
Mathematics Subject Classication (2000) 49M25 49K20 65F08
65F10 65F50 65N22
1 Introduction
In this paper we introduce a new class of structured matrices that we sug-
gest to call multilevel quasiseparable. It generalizes well-known quasiseparable
matrices to the multi-dimensional case. Under the multilevel paradigm, param-
eters that are used to represent a matrix of a higher hierarchy are themselves
multilevel quasiseparable of a lower hierarchy. The usual one-dimensional qua-
siseparable matrix is the one of the lowest hierarchy.
School of Mathematics, The University of Edinburgh, JCMB, The Kings Buildings, Edin-
burgh, Scotland EH9 3JZ, E-mail: [email protected], [email protected].
2 Jacek Gondzio, Pavel Zhlobich
Quasiseparable matrices found their application in orthogonal polynomials
[1], root nding [4], system theory [2] and other branches of applied mathe-
matics and engineering. Consequently, there has been major interest in them
in the scientic community, and many fast algorithms for working with them
have been developed in the last decade. Due to the similarity in represen-
tations of multilevel quasiseparable and usual quasiseparable matrices, it is
possible to extend the applicability of almost all the algorithms developed for
quasiseparable matrices to the new class. It is a celebrated fact that operations
with quasiseparable matrices can be performed in linear time with respect to
their size. We show that under some conditions this is also true for multilevel
quasiseparable matrices.
Natural applications of multilevel quasiseparable matrices are discretized
Partial Dierential Equations (PDEs) and related problems. There have al-
ready been recent attempts to use structured matrices similar to quasisepa-
rable ones to solve discretized PDEs [25, 19]. The very distinct feature of our
approach is the ability not only to solve systems in linear time, but also to
invert matrices and to perform arithmetic operations with them in linear time.
Therefore, we can solve, for instance, saddle-point systems with PDE matrices
that former approaches did not allow.
The structure of the paper is as follows: we start with a motivation for
this paper by considering a simple PDE-constrained optimization problem in
Section 2. This type of problems leads to systems in a saddle-point form. The
need in ecient solvers for such systems motivated this paper. In Section 3 we
briey review quasiseparable matrices and related results. Section 4 concerns
the extension of quasiseparable matrices to the multi-dimensional case. We
give there a formal denition of multilevel quasiseparable matrices, introduce
their parametrization and fast arithmetic with them. The last section presents
results of preliminary numerical experiments that empirically conrm linear
time complexity of algorithms with multilevel quasiseparable matrices.
2 Model PDE-constrained optimization problem.
As a simple model let us consider a distributed control problem which is com-
posed of a cost functional
min
u,f
1
2
|u u|
2
2
+|f|
2
(1)
to be minimized subject to a Poissons problem with Dirichlet boundary con-
ditions:

2
u = f in ,
u = g on .
(2)
The problem consists in nding a function u that satises the PDE con-
straint and is as close as possible in the L2 norm sense to a given function u
(desired state). The right-hand side of the PDE f is not xed and can be
designed in order to achieve the optimality of u. In general such a problem
Multilevel quasiseparable matrices in PDE-constrained optimization 3
is ill-posed and thus needs Tikhonov regularization that is introduced in the
form of the second term in (1).
The usual approach to solving the minimization problem (1) (2) numer-
ically is to introduce a weak formulation of it and to discretize the latter by
nite elements, see, for instance, [21, 12]. After these steps, the discrete analog
of (1) (2) reads
min
u
h
,f
h
1
2
|u
h
u|
2
2
+|f
h
|
2
2
, u
h
V
h
g
, f
h
V
h
0
,
subject to
_

u
h
v
h
=
_

v
h
f, v
h
V
h
0
,
(3)
where V
h
0
and V
h
g
are nite-dimensional vector spaces spanned by test func-
tions. If
1
, . . .
n
and
1
, . . .
n
,
n+1
, . . .
n+n
are bases in V
h
0
and V
h
g
,
respectively, and u
h
and v
h
are represented in them as
u
h
=
n+n

j=1
u
j

j
, f
h
=
n

j=1
f
j

j
,
then the matrix form of (3) is
min
u,f
1
2
u
T
Mu u
T
b +c +f
T
Mf ,
subject to Ku = Mf +d,
(4)
where M
ij
=
_

j
mass matrix, K
ij
=
_

j
stiness matrix,
b
j
=
_

j
u, c =
1
2
_

u
2
, d
j
=
n+n

k=n+1
u
k
_

j
, i, j = 1, . . . , n.
The Lagrangian associated with problem (4) is
/(u, f , ) =
1
2
u
T
Mu u
T
b +c +f
T
Mf +
T
(Ku Mf d),
where is a vector of Lagrange multipliers. The condition for a stationary
point of / dene u, f and via the solution of the linear system
_
_
2M 0 M
0 M K
T
M K 0
_
_
_
_
f
u

_
_
=
_
_
0
b
d
_
_
. (5)
This linear system is of saddle-point type. We refer to [3] for an extensive
survey of numerical methods for this type of systems.
Matrix in (5) is symmetric, usually very large but sparse due to nite ele-
ment discretization. Thus matrix-vector multiplication is cheap and Krylov
subspace methods such as Minimal Residual (MINRES) method of Paige
and Saunders [20] or Projected Preconditioned Conjugate Gradient (PPCG)
method of Gould, Hribar and Nocedal [15] are well suited. Their convergence
nevertheless depends highly on the choice of a good preconditioner [21].
4 Jacek Gondzio, Pavel Zhlobich
In this paper we propose a structured matrices approach to the solution
of systems of the type described above. First, let us write the block LDU
decomposition of (5):
_
_
2M 0 M
0 M K
T
M K 0
_
_
=
_
_
I 0 0
0 I 0
1
2
I KM
1
I
_
_

_
_
2M 0 0
0 M 0
0 0 S
_
_

_
_
I 0
1
2
I
0 I M
1
K
T
0 0 I
_
_
,
S =
_
1
2
M +KM
1
K
T
_
.
(6)
The hardest part in using this decomposition for solving a system is to compute
the Schur complement S of the (3,3) block and to solve the corresponding
system with it. Since matrix S is usually dense, this task is untractable if we
use its entrywise representation but it is feasible if we use a proper structured
representation. In Section 4 we will show that Schur complement indeed has a
structure and this structure is the one called multilevel quasiseparable. We will
give all necessary denitions and show that the use of multilevel quasiseparable
matrices leads to asymptotically (with respect to the size of the system) linear
storage and linear complexity algorithm for solving (5).
3 Quasiseparable matrices.
Matrix of a rank much smaller than its size is a discrete analog of a separable
function. More generally, we may think of a matrix with certain blocks being
low rank rather than the whole matrix. This case corresponds to a function
being separable in some subdomains of the whole domain. One simple example
of such a function is Greens function of the SturmLiouville equation (this
example will be considered in some more details later in this section). There
is no surprise that matrices with low rank blocks found their applications in
many branches of applied mathematics, especially in integral equations. There
are several related classes of rank-structured matrices that one can nd in the
scientic literature. Among the most well-known are semiseparable [14], quasi-
separable [9], H-matrices [16], H
2
-matrices [17] and mosaic-skeleton matrices
[23]. In the current paper our main tool will be quasiseparable matrices [24]
having low rank partitions in their upper and lower parts. The formal deni-
tion is given below.
Denition 1 (Rank denition of a block quasiseparable matrix) Let
A be a block matrix of block sizes n
k

n
k=1
then it is called block (r
l
, r
u
)-
quasiseparable if
max
K
rank A(K + 1 : N, 1 : K) r
l
, max
K
rank A(1 : K, K + 1 : N) r
u
,
K =
k

i=1
n
i
, N =
n

i=1
n
i
,
where r
l
(r
u
) is called lower (upper) order of quasiseparability.
Multilevel quasiseparable matrices in PDE-constrained optimization 5
In what follows we will call block (r
l
, r
u
)-quasiseparable matrices simply
quasiseparable for shortness.
A well-known example of a quasiseparable matrix is given by the discretized
Greens function of a one-dimensional dierential equation. Consider a regular
inhomogeneous Sturm-Liouville boundary value problem:
(p(x)u

q(x)u = f(x),
_

1
u(a) +
1
u(a) = 0,

2
u(b) +
2
u(b) = 0,
(7)
where functions p, p

, q are continuous on [a, b], p(x) > 0 on [a, b] and [

i
[ +
[
i
[ , = 0 for i = 1, 2. It is a classical result in ordinary dierential equations
that the solution f(x) of (7) is given by
u(x) =
b
_
a
g(x, )f() d,
where g(x, ) is the so-called Greens function. In this case it has an explicit
formula
g(x, ) =
1
p(a)W(a)
=
_
u
1
(x)u
2
(), a x ,
u
1
()u
2
(x), < x b,
(8)
with W(x) being Wronskian and u
1
(x) and u
2
(x) being two specic linearly
independent solutions of (7). It is easy to see from (8) that discretized Greens
function is a quasiseparable matrix of order one.
In order to exploit the quasiseparability of matrices in practice one must
have a low-parametric representation of them. There are many alternative
parametrisations of quasiseparable matrices all of which use O(N) parameters,
where N is the size of the matrix. Having such a parametrisation at hand one
can write most of the algorithms, e.g., inversion, LU or QR decomposition,
matrix-vector multiplication in terms of these parameters and, therefore, the
complexity of these algorithms is O(N). In this paper we will use the so-
called generator representation (Denition 2 below) proposed by Eidelman and
Gohberg [9]. There are several alternative parametrisations such as Givens-
weight [7] and Hierarchical Semiseparable (HSS) [6].
Denition 2 (Generator denition of a block quasiseparable matrix)
Let A be a block matrix of block sizes n
k

n
k=1
then it is called block (r
l
, r
u
)-
quasiseparable if it can be represented in the form
_

_
d
1
g
1
h
2
g
1
b
2
h
3
g
1
b
2
. . . b
n1
h
n
p
2
q
1
d
2
g
2
h
3
g
2
b
3
. . . b
n1
h
n
p
3
a
2
q
1
p
3
q
2
d
3
g
3
b
4
. . . b
n1
h
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
p
n
a
n1
. . . a
2
q
1
p
n
a
n1
. . . a
3
q
2
p
n
a
n1
. . . a
4
q
3
d
n
_

_
, (9)
6 Jacek Gondzio, Pavel Zhlobich
Table 1 Sizes of generators of a block quasiseparable matrix.
d
k
q
k
a
k
p
k
g
k
b
k
h
k
# of rows n
k
r
l
k
r
l
k
n
k
n
k
r
u
k1
r
u
k1
# of cols n
k
n
k
r
l
k1
r
l
k1
r
u
k
r
u
k
n
u
k
where parameters (called generators) d
k
, q
k
, a
k
, p
k
, g
k
, b
k
, h
k
are matrices
of sizes as in the table below.
Orders of quasiseparability r
l
and r
u
are maxima over the corresponding
sizes of generators:
r
l
= max
1kn1
r
l
k
, r
u
= max
1kn1
r
u
k
.
Remark 1 Generators are not unique, there are innitely many ways to repre-
sent the same quasiseparable matrix with dierent generators. For the relation
between dierent sets of generators see [11].
Remark 2 It is obvious that any scalar quasiseparable matrix can be converted
to the block form by simple aggregation of its generators.
Remark 3 Sizes r
l
k
and r
u
k
of generators are directly related to ranks of sub-
matrices in the lower and upper parts, respectively. Namely, if we let K to be
the block index: K =
k

i=1
n
i
, then
rank A(K + 1 : N, 1 : K) r
l
k
, rank A(1 : K, K + 1 : N) r
u
k
. (10)
Moreover, for any quasiseparable matrix there exists a set of generators with
sizes r
l
k
and r
u
k
that satisfy (10) with exact equalities (such generators are
called minimal ). For their existence and construction see [11].
One remarkable property of quasiseparable matrices is that this class is
closed under inversion [9, 22]. For instance, inverse of a banded matrix is not
banded but both matrices are quasiseparable. Due to the low-parametric rep-
resentation in Denition 2 most of the algorithms with quasiseparable matrices
have linear complexities. Table 2 lists these algorithms along with the corre-
sponding references.
Table 2 O(n) algorithms for quasiseparable matrices
A v LU QR A
1
A B AB A
F
[8] Theorem 1 [10] [9] [9] [11]
The ability to solve systems with quasiseparable matrices in O(n) oper-
ations is essential for the purpose of this paper. One of the ways to do it is
through the use of fast LU decomposition. We next derive LU decomposition
Multilevel quasiseparable matrices in PDE-constrained optimization 7
algorithm for a general block quasiseparable matrix in terms of the generators
it is represented with. A restricted version of this algorithm applicable to diag-
onal plus semiseparable matrices (a subclass of quasiseparable matrices) was
derived in [14] and some formulae of the current algorithm are similar to the
ones in the inversion algorithm given in [9]. Still, to the best of our knowledge,
current LU decomposition algorithm has never been published and may be
useful to those who need a fast system solver for quasiseparable matrices. The
complexity of the algorithm is O(N) and it is valid in the strongly regular case
(i.e. all of the principal leading minors are non-vanishing). First, let us note
that quasiseparable structure of the original matrix implies the quasiseparable
structure of L and U factors. The theorem below makes this statement precise
and, in addition, relates generators of an original matrix to generators of the
factors.
Theorem 1 Let A be a strongly regular N N block (r
l
, r
u
)-quasiseparable
matrix given by generators d
k
, q
k
, a
k
, p
k
, g
k
, b
k
, h
k
as in (9). Let A = LU be
its block LU decomposition of the same block sizes. Then
(i) Factors L and U are (r
l
, 0) and (0, r
u
)-quasiseparable. Moreover, r
l
k
(L) =
r
l
k
(A) and r
u
k
(U) = r
u
k
(A) for all k = 1, . . . , n 1.
(ii) L and U are parametrized by the generators I
k
, q
k
, a
k
, p
k
, 0, 0, 0 and

d
k
,
0, 0, 0, g
k
, b
k
, h
k
, where I
k
are identity matrices of sizes n
k
n
k
and new
parameters q
k
,

d
k
and g
k
can be computed using the following algorithm:
Algorithm 1 Fast quasiseparable LU decomposition.
Input: d
k
, q
k
, a
k
, p
k
, g
k
, b
k
, h
k
1:

d
1
= d
1
, q
1
= q
1

d
1
1
, g
1
= g
1
, f
1
= q
1
g
1
2: for k = 2 to n 1 do
3:

d
k
= d
k
p
k
f
k1
h
k
4: q
k
= (q
k
a
k
f
k1
h
k
)

d
1
k
5: g
k
= g
k
p
k
f
k1
b
k
6: f
k
= a
k
f
k1
b
k
+ q
k
g
k
7: end for
8:

d
n
= d
n
p
n
f
n1
h
n
.
Output:

d
k
, q
k
, g
k
Proof Statement (i) of the theorem follows from statement (ii), so we only
need to prove the latter part.
Denote, as usual, by K the block index: K =
k

i=1
n
i
and note that quasi-
separable representation (9) implies the following recurrence relation between
8 Jacek Gondzio, Pavel Zhlobich
the blocks of A:
A =
_
A(1 : K, 1 : K) G
k
H
k+1
P
k+1
Q
k
A(K + 1 : N, K + 1 : N)
_
;
Q
1
= q
1
, Q
k
= [a
k
Q
k1
q
k
], k = 2, . . . n 1;
P
n
= p
n
, P
k
= [p
k
; P
k+1
a
k
], k = n 1, . . . 2;
G
1
= g
1
, G
k
= [G
k1
b
k
; g
k
], k = 2, . . . n 1;
H
n
= h
n
, H
k
= [h
k
b
k
H
k+1
], k = n 1, . . . 2.
(11)
The proof will be constructed by induction. We will show that for each k
_
A
k
11
G
k
H
k+1
P
k+1
Q
k

_
=
_
L
k
11
0
P
k+1

Q
k

_

_
U
k
11

G
k
H
k+1
0
_
. (12)
For k = 1 we get from (12):
d
1
=

d
1
,
P
2
Q
1
= P
2

Q
1

d
1
,
G
1
H
2
=

G
1
H
2
,
=

d
1
= d
1
,
q
1
= q
1

d
1
1
,
g
1
= g
1
.
Let us introduce an auxiliary parameter f
k
=

Q
k

G
k
. It is easy to show by
using (11) that f
k
satises the recurrence relation
f
1
= q
1
g
1
, f
k
= a
k
f
k1
b
k
+ q
k
g
k
. (13)
Assume that (12) holds for some xed k, then it holds for k + 1 if
d
k+1
= [p
k+1

Q
k
1] [

G
k
h
k+1
;

d
k+1
], (14)
P
k+2
q
k+1
= P
k+2

Q
k+1
[

G
k
h
k+1
;

d
k+1
], (15)
g
k+1
H
k+2
= [p
k+1

Q
k
1]

G
k+1
H
k+2
. (16)
The rst equality simplies
d
k+1
= p
k+1

Q
k

G
k
h
k+1
+

d
k+1
= p
k+1
f
k
h
k+1
+

d
k+1
.
The second equality (15) holds if
q
k+1
=

Q
k+1
[

G
k
h
k+1
;

d
k+1
] =
= [a
k+1

Q
k
q
k+1
] [

G
k
h
k+1
;

d
k+1
] = a
k+1
f
k
h
k+1
+ q
k+1

d
k+1
.
Finally, the last equality (16) is true if
g
k+1
= [p
k+1

Q
k
1]

G
k+1
= [p
k+1

Q
k
1][

G
k
b
k+1
; g
k+1
] = p
k+1
f
k
b
k+1
+ g
k+1
.
Matrix

d
k+1
is invertible because, by our assumption, matrix A is strongly
regular. Hence, we conclude that (12) is true also for index k +1 if generators
q
k+1
,

d
k+1
and g
k+1
are those computed in Algorithm 1. The assertion of the
theorem holds by induction. .
Multilevel quasiseparable matrices in PDE-constrained optimization 9
4 Multilevel quasiseparable matrices.
In the previous section we have emphasized the relation between second-order
ordinary dierential equations and quasiseparable matrices. Let us now look
at the problem in two dimensions that extends (7):

x
_
p(x, y)

x
u(x, y)
_
+

y
_
q(x, y)

y
u(x, y)
_
r(x, y)u(x, y) = f(x, y)
for (x, y) , where = [0, 1] [0, 1] with homogeneous Dirichlet boundary
conditions. The standard ve-point or nine-point discretization of this problem
on a n m uniform grid leads to a system of linear algebraic equations of the
form
_

_
A
1
B
1
C
1
A
2
B
2
C
2
.
.
.
.
.
.
.
.
.
.
.
.
B
m1
C
m1
A
m
_

_
_

_
u
1
u
2
.
.
.
.
.
.
u
m
_

_
=
_

_
f
1
f
2
.
.
.
.
.
.
f
m
_

_
, (17)
where we have assumed that u
i
are the discretized unknowns along the ith
column of the nm grid. In this case each of A
i
, B
i
, and C
i
is an nn matrix
and the whole block matrix has size nm nm. Furthermore, each of the A
i
,
B
i
, and C
i
is a tridiagonal matrix.
One way to solve system (17) is to do block LU factorization assuming
that it exists. When we eliminate the block entry C
1
we get in the position
occupied by A
2
the new block
S
2
= A
2
C
1
A
1
1
B
1
.
Observe that even though all the individual matrices on the right-hand side
of the expression are tridiagonal, A
1
1
is not, and hence S
2
is a dense (non-
sparse) matrix. At the next step of Gaussian elimination we would use S
1
as the
pivot block to eliminate C
2
. Now in the position occupied by the block A
3
we
would get the matrix S
3
= A
3
C
2
S
1
2
B
2
. Again, since S
2
is a dense matrix,
in general S
1
2
will be dense, and therefore S
3
will also be a dense matrix.
What this implies is that during LU factorization of the sparse matrix, we
will produce ll-in quickly that causes us to compute inverses (and hence LU
factorizations) of dense nn matrices. If we assume that these dense matrices
have no structure, then we would need O(n
3
) ops for that operation alone.
Therefore, it follows that one would require at least O(mn
3
) ops to compute
block LU factorization of the system matrix.
At rst glance it seems that block LU factorization is not a wise approach
as it does not use any kind of ll-in minimisation reordering. However, there is
a rich structure in successive Schur complements S
k
that can be exploited to
speed-up computations. In fact, it has been conjectured that if one looks at the
o-diagonal blocks of these matrices (S
2
, S
3
, etc.), then their -rank (number
of singular values not greater than ) is going to be small. This conjecture has
10 Jacek Gondzio, Pavel Zhlobich
Fig. 1 Growth of -rank (quasiseparable order) with size for dierent .
been justied by the fact that, for example, S
k
can be viewed approximately
(especially in the limit as n becomes large) as a sub-block of the discretized
Greens function of the original PDE. It is known from the theory of elliptic
PDEs (see, for instance, [13]) that under some mild constraints the Greens
function is smooth away from the diagonal singularity. This, in turn, implies
that the numerical ranks of the o-diagonal blocks of S
1
are small. There have
been several recent attempts to quantify this observation. It has been proved
in [5] that -rank of the o-diagonal blocks of Schur complements in the LU
decomposition of the regular 2D-Laplacian are bounded by
1 + 8 ln
4
_
18

_
.
This bound is not eective for any reasonable size of the problem because of
the 4th power in the logarithm (for instance for =1.0e-6 it gives 6.2e+5).
However, numerical experiments with Laplacian discretized by the usual ve-
point stencil suggest much lower bound, see Figure 1.
Careful exploitation of the low-rank structure of Schur complements re-
duces the complexity of block LU factorization algorithm to linear O(mn).
Indeed, as it has been mentioned in the previous section, all the operations
with quasiseparable matrices that are needed at each step of the algorithm
have O(n) complexity and there are m steps altogether. First algorithms of
this type were recently developed in [25, 19]. These algorithms, although ef-
cient in the case of simple PDEs, are not applicable to PDE-constrained
optimization problems, such as the one described in Section 2. Saddle-point
matrices arising there are dense block matrices with blocks themselves being
discretized PDE matrices. To use block LU factorization in this case one would
need to invert PDE matrices of type (17). We next show that it is possible by
dening the right structured representation of inverses of PDE matrices.
Multilevel quasiseparable matrices in PDE-constrained optimization 11
Denition 3 (Multilevel quasiseparable matrix) Matrix A is called d-
level quasiseparable if it admits the representation
_

_
D
1
G
1
H
2
G
1
B
2
H
3
G
1
B
2
. . . B
n1
H
n
P
2
Q
1
D
2
G
2
H
3
G
2
B
3
. . . B
n1
H
n
P
3
A
2
Q
1
P
3
Q
2
D
3
G
3
B
4
. . . B
n1
H
n
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
P
n
A
n1
. . . A
2
Q
1
P
n
A
n1
. . . A
3
Q
2
P
n
A
n1
. . . A
4
Q
3
D
n
_

_
,
where each of the parameters Q
k
, A
k
, P
k
, D
k
, G
k
, B
k
, H
k
is a block matrix
of blocks being (d 1)-level quasiseparable matrices. 1-level quasiseparable
matrix is a usual quasiseparable matrix, see Denition 1.
Let us give some simple examples of multilevel quasiseparable matrices to
justify the usefulness of the proposed structure.
Example 1 (Block banded matrices with banded blocks) Any banded matrix is
automatically quasiseparable. Similarly, block banded matrix is 2-level quasi-
separable if blocks are themselves banded matrices. As an example, consider
2D Laplacian on the square discretized by Q1 nite elements:
K =
_

_
A B
B
.
.
.
.
.
.
.
.
.
.
.
.
B
B A
_

_
, A =
1
3
_

_
8 1
1
.
.
.
.
.
.
.
.
.
.
.
.
1
1 8
_

_
, B =
1
3
_

_
1 1
1
.
.
.
.
.
.
.
.
.
.
.
.
1
1 1
_

_
. (18)
Example 2 (Tensor product of quasiseparable matrices) Matrix AB is 2-level
quasiseparable if A and B are 1-level quasiseparable. Hence, 2D mass matrix
M on the square discretized by Q1 nite elements:
M = T T, T =
1
6
_

_
4 1
1
.
.
.
.
.
.
.
.
.
.
.
.
1
1 4
_

_
(19)
is 2-level quasiseparable. Its inverse M
1
is also 2-level quasiseparable due to
the tensor identity (AB)
1
= A
1
B
1
and properties of 1-level quasise-
parable matrices.
Quasiseparable orders of multilevel quasiseparable matrices may dier at
each level. For instance, it follows immediately from the properties of quasise-
parable matrices that inverse of a 2-level quasiseparable matrix has a represen-
tation as in Denition 3 for some generators Q
k
, A
k
, P
k
, D
k
, G
k
, B
k
, H
k
.
However, there is no guarantee that entries of this generators are quasisepa-
rable matrices of the same order as for the original matrix. 2D Laplace matrix
in Example 1 is 2-level quasiseparable with generators of quasiseparable order
12 Jacek Gondzio, Pavel Zhlobich
one but as Figure 1 illustrates, generators of LU factors (and, therefore, in-
verse) have orders larger than one. However, if it is true that numerical orders
(ranks) do not dier much from the original orders we may say that inverses
of multilevel quasiseparable matrices retain the structure approximately. This
observation, in its turn, leads to algorithms of linear complexity.
To compute LU factorization of a 2-level quasiseparable matrix we may
use Algorithm 1. It consists of approximately 12n arithmetic operations with
1-level quasiseparable matrices. Sum or dierence of rank-r
1
and rank-r
2
ma-
trices may have rank r
1
+r
2
. Similarly, any arithmetic operation (except inver-
sion) performed on quasiseparable matrices may add their orders in the worst
case. Therefore, to use multilevel quasiseparable structure eciently we need
some kind of order compression algorithm. Such algorithm was proposed in
[11], we list it below for the convenience of the reader. This algorithm consti-
tutes the essential part of the proposed method.
Algorithm 2 Quasiseparable order reduction.
Input: q
k
, a
k
, p
k
of sizes r
k
n
k
, r
k
r
k1
and n
k
r
k1
, respectively.
1: Using LQ or SVD decomposition determine matrices L
1
and q

1
of sizes r
1
r

1
and
r

1
n
1
, respectively, such that q
1
= L
1
q

1
, where q

1
(q

1
)

= I
r

1
and r

1
= rank q
1
.
2: for k = 2 to n 1 do
3: Using LQ or SVD decomposition determine matrices L
k
and V
k
of sizes r
k
r

k
and
r

k
(r

k1
+ n
k
), respectively, such that [a
k
L
k1
q
k
] = L
k
V
k
, where V

k
(V

k
)

= I
r

k
and r

k
= rank[a
k
L
k1
q
k
].
4: Split V
k
into the new generators a

k
, q

k
of sizes r

k
r

k1
, r

k
n
k
, respectively:
V
k
= [a

k
q

k
].
5: p

k
= p
k
L
k1
6: end for
7: p

n
= p
n
L
n1
8: Using QR or SVD decomposition determine matrices p

n
and S
n1
of sizes n
n
r

n1
and r

n1
r

n1
, respectively, such that p

n
= p

n
S
n1
, where (p

n
)

n
= I
r

n1
and
r

n1
= rank p

n
.
9: for k = n 1 to 2 do
10: Using QR or SVD decomposition determine matrices U
k
and S
k1
of sizes (n
k
+
r

k
) r

k1
and r

k1
r

k1
, respectively, such that [p

k
; S
k
a

k
] = U
k
S
k1
, where
(U
k
)

U
k
= I
r

k1
and r

k1
= rank p

k
.
11: Split U
k
into the new generators p

k
, a

k
of sizes n
k
r

k1
, r

k
r

k1
, respectively:
U
k
= [p

k
; ; a

k
].
12: q

k
= S
k
q

k
13: end for
14: q

1
= S
1
q

1
Output: New generators q

k
, a

k
, p

k
of minimal possible sizes r

k
n
k
, r

k
r

k1
and
n
k
r

k1
, respectively.
In exact arithmetic, new sizes r

k
of generators match ranks of the corre-
sponding submatrices in the quasiseparable matrix (see [11] for the proof). In
oating point arithmetic we may use truncated SVD or rank-revealing LQ/QR
to decrease sizes of generators to the value of the -rank of the corresponding
Multilevel quasiseparable matrices in PDE-constrained optimization 13
submatrix. When -ranks are small, complexity of Algorithm 2 is O(n) and,
hence, the following fact takes place.
Complexity of algorithms listed in Table 2 with 2-level quasiseparable ma-
trices (Denition 3) is linear in the size if quasiseparable orders of 1-level
quasiseparable matrices stay small during computations.
In the next section we will demonstrate practical properties of multi-
level quasiseparable matrices approach applied to solving PDEs and PDE-
constrained optimization problems.
5 Numerical results.
We illustrate our new method with two dierent examples: 2D Laplaces equa-
tion with Dirichlet boundary conditions and distributed control problem with
a constraint in the form of 2D Poissons equation. We show that the proposed
method can be used either directly or as a preconditioner for an iterative
method depending on the level of order reduction used in Algorithm 2. Al-
though we consider symmetric problems only, our approach does not require
symmetry and only depends on low rank properties of operators. Therefore,
it can also be applied to convectiondiusion and Helmholtz problems. We do
not consider 3D problems as it is not clear yet how to adapt SVD and LQ/QR
decomposition used in Algorithm 2 to the block case.
Our software is written in C++ and is compiled with GCC compiler v.4.2.1.
It is freely available for download at http://www.maths.ed.ac.uk/ERGO/
software.html. All tests were done on a 2.4 GHz Intel Core 2 Duo with
4 Gb of RAM.
Example 3 Let = [0, 1]
2
and consider the problem

2
u = 0, in ,
u = u[

, on ,
(20)
where
u =
_

_
sin(2y) if x = 0,
sin(2y) if x = 1,
0 otherwise.
Let us introduce a square reqular grid on [0, 1]
2
of mesh size h = 1/(n+1)
and discretize equation (20) on this grid by Q1 nite elements. Then equation
(20) reduces to a matrix equation with matrix K in (18). This matrix, as it
was noted before, is 2-level quasiseparable and, therefore, we can apply fast
LDU decomposition Algorithm 1 with fast order-reduction Algorithm 2 at
each step. Maximal quasiseparable order can be chosen adaptively in SVD
but we set it manually to 4 and 8 for simplicity. Larger order corresponds
to better approximation and, hence, more accurate result although it makes
algorithm slower. Results of the computational tests for dierent mesh sizes
are presented in Tables 3 and 4. Note, that our solver was implemented for
14 Jacek Gondzio, Pavel Zhlobich
unsymmetric problems and did not use symmetry of the matrix in the current
example. Exploiting the symmetry would roughly halve the time and memory
used.
Table 3 Time and memory used by direct quasiseparable LDU decomposition based solver
applied to the problem in Example 3, quasiseparable order (maximal o-diagonal rank) is
set to 4 during computations.
n LDU mem, Mb solve
Auf
2
f
2
2
12
0.11 3 0.003 8.22e-05
2
14
0.51 12 0.011 1.85e-04
2
16
2.18 51 0.043 3.93e-04
2
18
9.00 210 0.169 6.91e-04
2
20
36.88 847 0.677 8.81e-04
Table 4 Time and memory used by direct quasiseparable LDU decomposition based solver
applied to the problem in Example 3, quasiseparable order (maximal o-diagonal rank) is
set to 8 during computations.
n LDU mem, Mb solve
Auf
2
f
2
2
12
0.19 4 0.003 3.31e-09
2
14
0.91 19 0.013 6.19e-08
2
16
4.09 83 0.052 5.72e-07
2
18
17.44 342 0.209 2.33e-06
2
20
72.83 1388 0.852 5.41e-06
The analysis of results collected in Tables 3 and 4 reveals the linear depen-
dence of time and memory used by the algorithm on the size of the problem.
The computational experience supports our claim about asymptotically linear
complexity of algorithms with multilevel quasiseparable matrices. Note also the
linear dependence of time and the sub-linear dependence of memory on the
maximal quasiseparable order used during computations.
If we truncate the order of quasiseparable matrices in the quasiseparable
LDU decomposition to a very small number, the decomposition becomes less
accurate. However, the algorithm becomes much faster and it is tempting to
try this approximate LDU decomposition as a preconditioner in the precon-
ditioned conjugate gradient method. Table 5 bellow illustrates the numerical
performance of this approach. Results indicate that the PCG method precon-
ditioned with the inexact LDU decomposition is twice as fast as more accurate
LDU decomposition alone. It is enough to truncate quasiseparable order to 24
to force PCG to converge in less than 10 iterations.
Solution of problem in Example 3 obtained with PCG method with quasi-
separable preconditioner is given in Figure 2 below.
Multilevel quasiseparable matrices in PDE-constrained optimization 15
Table 5 Time and number of iterations (in brackets) used by PCG when applied to the
problem in Example 3, preconditioner is based on approximate quasiseparable LDU decom-
position with order tolerance r, system is solved to the 1.0e-08 accuracy.
n r LDU PCG total time
2
12
1 0.05 0.0359 (9) 0.09
2 0.07 0.0252 (6) 0.09
2
14
1 0.21 0.212 (14) 0.42
2 0.30 0.144 (9) 0.45
2
16
3 1.76 0.482 (7) 2.24
4 2.18 0.282 (4) 2.47
2
18
3 7.17 2.96 (11) 10.13
4 9.05 1.98 (7) 11.03
2
20
4 37.07 10.1 (9) 47.22
5 45.63 8.19 (7) 53.83
Fig. 2 Solution of problem in Example 3.
Example 4 Let = [0, 1]
2
and consider the problem
min
u,f
1
2
|u u|
2
L
2
()
+|f|
2
L
2
()
,
s.t.
2
u = f, in ,
u = u[

, on ,
(21)
where
u =
_
(2x 1)
2
(2y 1)
2
if (x, y) [0,
1
2
]
2
,
0 otherwise.
16 Jacek Gondzio, Pavel Zhlobich
As it was already discussed in Section 2, problem (21) can be solved by the
discretization by Q1 nite elements. Discretized system of equation, in this
case, becomes
_
_
2M 0 M
0 M K
T
M K 0
_
_
_
_
f
u

_
_
=
_
_
0
b
d
_
_
. (22)
After we eliminate variables f and u in equation (22) it reduces to
_
1
2
M +KM
1
K
T
_
= y. (23)
This system is sometimes called normal equation in optimization community.
Computing matrix on the left side of (23) is prohibitively expensive as it is
usually very large and dense. Even more prohibitively expensive is solving the
system with this matrix as it is O(n
3
) algorithm. One common approach is
to drop either
1
2
M or KM
1
K
T
in the equation and to solve system with
the remaining matrix, thus, using it as a preconditioner (see, for instance,
[21]). However, this preconditioner would perform badly if none of the terms
is dominant. We propose a completely dierent approach. Matrices M and
K in (23) are multilevel quasiseparable (see Examples 1 and 2). Thus, we
can compute the Schur complement matrix in the left hand side of (23) in
quasiseparable arithmetic using O(n) arithmetic operations if quasiseparable
orders stay small during computations. The last condition is true in practice
because Schur complement is itself a discretized elliptic operator similar to
Laplacian. In the previous example we have shown that it is more eective
to use the proposed approach as a preconditioner rather than a direct solver.
Thus, our approach to solving equation (22) consists of the following steps:
1. Invert matrix M and form matrix S =
1
2
M+KM
1
K
T
in quasiseparable
matrices arithmetic.
2. Compute an approximate LDU decomposition of S using Algorithm 1 and
some order reduction tolerance.
3. Use PCG method to solve S = y with approximate LDU decomposition
as a preconditioner.
4. Exploit computed , M
1
and the factorization (6) to nd f and u.
We have realized the proposed approach in practice. Tables 6 and 7 gather
the numerical results we have obtained while solving the problem in Example
4.
Analyzing the results presented in Tables 6 and 7 we conclude that our pre-
conditioner is mesh-independent and, as in the case of simple PDE problem,
CPU time for the overall algorithm grows linearly with the size of the prob-
lem. Control function f and obtained solution u computed with the proposed
algorithm for dierent values of the regularization parameter are presented
in Figure 3.
Multilevel quasiseparable matrices in PDE-constrained optimization 17
Table 6 Time required to construct the normal equation matrix S in the left hand side of
(23) with = 10
2
, compute its approximate LDU decomposition and solve system with
it to a specied tolerance using PCG method (number of iterations is given in brackets).
Quasiseparable order used in order reduction algorithm was set to 1.
= 10
2
, tol = 10
4
n S LDU(S) PCG total time
12 0.73 1.12 0.0839 (3) 1.93
14 3.23 5.01 0.337 (3) 8.57
16 13.57 21.12 1.37 (3) 36.05
18 55.75 86.85 3.91 (2) 146.51
= 10
2
, tol = 10
8
n S LDU(S) PCG total time
12 0.71 1.11 0.13 (5) 1.95
14 3.23 5.00 0.534 (5) 8.76
16 13.58 21.11 2.17 (5) 36.86
18 55.79 86.83 8.76 (5) 151.38
Table 7 Time required to construct the normal equation matrix S in the left hand side of
(23) with = 10
5
, compute its approximate LDU decomposition and solve system with
it to a specied tolerance using PCG method (number of iterations is given in brackets).
Quasiseparable order used in order reduction algorithm was set to 1.
= 10
5
, tol = 10
4
n S LDU(S) PCG total time
12 0.71 1.12 0.0344 (1) 1.86
14 3.23 4.98 0.14 (1) 8.35
16 13.57 21.05 0.566 (1) 35.19
18 55.83 86.58 2.28 (1) 144.70
= 10
5
, tol = 10
8
n S LDU(S) PCG total time
12 0.71 1.11 0.058 (2) 1.88
14 3.23 4.98 0.238 (2) 8.45
16 13.57 21.05 0.966 (2) 35.59
18 55.79 86.71 3.9 (2) 146.40
6 Conclusion.
In this paper we have introduced a new class of rank-structured matrices called
multilevel quasiseparable. This matrices are low-parametric in the sense that
only O(n) parameters are needed to store a matrix. Moreover, arithmetic op-
erations and matrix decompositions can be performed in O(n) oating-point
operations. Multilevel quasiseparable matrices extend the applicability of well-
known quasiseparable matrices [9] from one-dimensional to multidimensional
problems. In particular, we have shown that multilevel quasiseparable struc-
ture is well-suited to the description of discretized elliptic operators in 2D.
18 Jacek Gondzio, Pavel Zhlobich
Fig. 3 Plots of the state, u, and the control, f, for = 10
2
and 10
2
.
(a) u, = 10
2
,
u
h
u
h

2
u
h

2
=0.03 (b) f, = 10
2
(c) u, = 10
2
,
u
h
u
h

2
u
h

2
=1.9e-4 (d) f, = 10
2
To demonstrate the usefulness of the new class of matrices we considered
distributed control problem with a constraint in the form of a partial dieren-
tial equation. Such problems were introduced by Lions and Mitter in [18]. In
the course of solving them large-scale block matrices of saddle-point type arise.
A straightforward way of solving these systems is to use block LU factoriza-
tion. This is impractical as it requires direct inversion of large PDE matrices
and further manipulations with them. However, we have shown that inverses
of PDE matrices can be approximated by multilevel quasiseparable matrices
with any desired accuracy and, hence, what was impossible in dense linear
algebra became possible in structured linear algebra.
Development of numerical methods for solving systems of saddle-point type
is an important subject of modern numerical linear algebra, see an excellent
survey paper by Benzi, Golub, and Liesen [3] for details. A large amount of
work has been done on developing ecient preconditioners for such systems. In
the current paper we have also proposed a new ecient mesh-independent pre-
conditioner for saddle-point systems arising in PDE-constrained optimization.
Performance of the new preconditioner is studied numerically.
There are several open questions left for future research. In particular, it is
not clear what the range of applicability of multilevel quasiseparable matrices
is. We have only considered problems with symmetric dierential operators in
Multilevel quasiseparable matrices in PDE-constrained optimization 19
our examples. However, theory does not require symmetry in the matrix and
it may very well be that multilevel quasiseparable matrices are applicable to
convectiondiusion and other unsymmetric operators as well. Next, robust
techniques are needed to make the approach work on non-tensor grids. The
biggest question though is how to extend order reduction Algorithm 2 to the
block case to be able to use multilevel quasiseparable matrices in 3D and higher
dimensions.
References
1. T. Bella, Y. Eidelman, I. Gohberg, V. Olshevsky, and P. Zhlobich. Classications of
recurrence relations via subclasses of (H,m)-quasiseparable matrices. In Numerical Lin-
ear Algebra in Signals, Systems and Control, volume XV of Lecture Notes in Electrical
Engineering, pages 2353. SpringerVerlag, 2011.
2. T. Bella, V. Olshevsky, and P. Zhlobich. Signal ow graph approach to inversion of
(H,m)-quasiseparable-Vandermonde matrices and new lter structures. Linear Algebra
and its Applications, 432:20322051, 2010.
3. M. Benzi, G. Golub, and J. Liesen. Numerical solution of saddle point problems. Acta
Numerica, 14:1137, May 2005.
4. D. Bini, F. Daddi, and L. Gemignani. On the shifted QR iteration applied to companion
matrices. Electronic Transactions on Numerical Analysis, 18:137152, 2004.
5. S. Chandrasekaran, P. Dewilde, M. Gu, and N. Somasunderam. On the numerical rank
of the o-diagonal blocks of Schur complements of discretized elliptic PDEs. SIAM
Journal on Matrix Analysis and Applications, 31:22612290, 2010.
6. S. Chandrasekaran, M. Gu, and T. Pals. A Fast ULV Decomposition Solver for Hi-
erarchically Semiseparable Representations. SIAM Journal on Matrix Analysis and
Applications, 28(3):603622, 2006.
7. S. Delvaux and M. Van Barel. A Givens-weight representation for rank structured
matrices. SIAM Journal on Matrix Analysis and Applications, 29(4):11471170, 2007.
8. Y. Eidelman and I. Gohberg. Linear complexity inversion algorithms for a class of
structured matrices. Integral Equations and Operator Theory, 35:2852, 1999.
9. Y. Eidelman and I. Gohberg. On a new class of structured matrices. Integral Equations
and Operator Theory, 34:293324, 1999.
10. Y. Eidelman and I. Gohberg. A modication of the Dewilde-van der Veen method for
inversion of nite structured matrices. Linear Algebra and its Applications, 343:419
450, 2002.
11. Y. Eidelman and I. Gohberg. On generators of quasiseparable nite block matrices.
Calcolo, 42:187214, 2005.
12. H.C. Elman, D.J. Silvester, and A.J. Wathen. Finite elements and fast iterative solvers:
with applications in incompressible uid dynamics. Oxford University Press, USA, 2005.
13. G.B. Folland. Introduction to partial dierential equations. Princeton University Press,
Princeton, NJ, 1995.
14. I. Gohberg, T. Kailath, and I. Koltracht. Linear complexity algorithms for semiseparable
matrices. Integral Equations and Operator Theory, 8(6):780804, 1985.
15. N.I.M. Gould, M.E. Hribar, and J. Nocedal. On the solution of equality constrained
quadratic programming problems arising in optimization. SIAM Journal on Scientic
Computing, 23(4):13761395, 2002.
16. W. Hackbusch. A sparse matrix arithmetic based on H-matrices. Part I: Introduction
to H-matrices. Computing, 62:89108, 1999.
17. W. Hackbusch, B. Khoromskij, and S. Sauter. On H
2
-matrices. In H. Bungartz,
R. Hoppe, and C. Zenger, editors, Lectures on Applied Mathematics, pages 929.
SpringerVerlag, Berlin, 2000.
18. J.L. Lions and S.K. Mitter. Optimal control of systems governed by partial dierential
equations, volume 396. Springer Berlin, 1971.
20 Jacek Gondzio, Pavel Zhlobich
19. P.G. Martinsson. A fast direct solver for a class of elliptic partial dierential equations.
Journal of Scientic Computing, 38(3):316330, 2009.
20. C.C. Paige and M.A. Saunders. Solution of sparse indenite systems of linear equations.
SIAM Journal on Numerical Analysis, pages 617629, 1975.
21. T. Rees, H.S. Dollar, and A.J. Wathen. Optimal solvers for PDE-constrained optimiza-
tion. SIAM Journal on Scientic Computing, 32(1):271298, 2010.
22. G. Strang and T. Nguyen. The interplay of ranks of submatrices. SIAM review,
46(4):637646, 2004.
23. E. Tyrtyshnikov. Mosaic-skeleton approximations. Calcolo, 33(1):4757, 1996.
24. R. Vandebril, M. Van Barel, and N. Mastronardi. A note on the representation and
denition of semiseparable matrices. Numerical Linear Algebra with Applications,
12(8):839858, 2005.
25. J. Xia, S. Chandrasekaran, M. Gu, and X.S. Li. Superfast multifrontal method for
large structured linear systems of equations. SIAM Journal on Matrix Analysis and
Applications, 31(3):13821411, 2009.