Multilevel quasiseparable matrices in PDE-constrained optimization

Multilevel quasiseparable matrices in PDE-constrained optimization

Jacek Gondzio School of Mathematics, The University of Edinburgh, JCMB, The King’s Buildings, Edinburgh, Scotland EH9 3JZ, 11email: J.Gondzio@ed.ac.uk, P.Zhlobich@ed.ac.uk.    Pavel Zhlobich School of Mathematics, The University of Edinburgh, JCMB, The King’s Buildings, Edinburgh, Scotland EH9 3JZ, 11email: J.Gondzio@ed.ac.uk, P.Zhlobich@ed.ac.uk.
Abstract

Optimization problems with constraints in the form of a partial differential 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 Green’s function of a one-dimensional differential 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 preconditioners 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, preconditioning, optimal control, linear systems, quasiseparable matrices
Msc:
49M25 49K20 65F08 65F10 65F50 65N22

1 Introduction

In this paper we introduce a new class of structured matrices that we suggest to call multilevel quasiseparable. It generalizes well-known quasiseparable matrices to the multi-dimensional case. Under the multilevel paradigm, parameters that are used to represent a matrix of a higher hierarchy are themselves multilevel quasiseparable of a lower hierarchy. The usual one-dimensional quasiseparable matrix is the one of the lowest hierarchy.

Quasiseparable matrices found their application in orthogonal polynomials BEGOZ10 (), root finding BDG04 (), system theory BOZ10b () and other branches of applied mathematics and engineering. Consequently, there has been major interest in them in the scientific community, and many fast algorithms for working with them have been developed in the last decade. Due to the similarity in representations 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 Differential Equations (PDEs) and related problems. There have already been recent attempts to use structured matrices similar to quasiseparable ones to solve discretized PDEs XCGL09 (); M09 (). 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 efficient solvers for such systems motivated this paper. In Section 3 we briefly review quasiseparable matrices and related results. Section 4 concerns the extension of quasiseparable matrices to the multi-dimensional case. We give there a formal definition of multilevel quasiseparable matrices, introduce their parametrization and fast arithmetic with them. The last section presents results of preliminary numerical experiments that empirically confirm 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 composed of a cost functional

 minu,f12∥u−ˆu∥22+β∥f∥2 (1)

to be minimized subject to a Poisson’s problem with Dirichlet boundary conditions:

 −∇2u=f in Ω,u=g on ∂Ω. (2)

The problem consists in finding a function that satisfies the PDE constraint and is as close as possible in the L2 norm sense to a given function (“desired state”). The right-hand side of the PDE is not fixed and can be designed in order to achieve the optimality of . In general such a problem 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) numerically is to introduce a weak formulation of it and to discretize the latter by finite elements, see, for instance, RDW10 (); ESW05 (). After these steps, the discrete analog of (1) – (2) reads

 minuh,fh12∥uh−ˆu∥22+β∥fh∥22,uh∈Vhg,fh∈Vh0, (3) subject to ∫Ω∇uh⋅∇vh=∫Ωvhf,∀vh∈Vh0,

where and are finite-dimensional vector spaces spanned by test functions. If and are bases in and , respectively, and and are represented in them as

 uh=n+∂n∑j=1ujϕj,fh=n∑j=1fjϕj,

then the matrix form of (3) is

 min→u,→f12→uTM→u−→uT→b+c+β→fTM→f, (4) subject to K→u=M→f+→d,

where — mass matrix, — stiffness matrix, , , , .

The Lagrangian associated with problem (4) is

 L(→u,→f,→λ)=12→uTM→u−→uTb+c+β→fTM→f+→λT(K→u−M→f−→d),

where is a vector of Lagrange multipliers. The condition for a stationary point of define , and via the solution of the linear system

 ⎡⎢⎣2βM0−M0MKT−MK0⎤⎥⎦⎡⎢ ⎢⎣→f→u→→λ⎤⎥ ⎥⎦=⎡⎢ ⎢⎣→0→b→d⎤⎥ ⎥⎦. (5)

This linear system is of saddle-point type. We refer to BGL05 () for an extensive survey of numerical methods for this type of systems.

Matrix in (5) is symmetric, usually very large but sparse due to finite element discretization. Thus matrix-vector multiplication is cheap and Krylov subspace methods such as Minimal Residual (MINRES) method of Paige and Saunders PS75 () or Projected Preconditioned Conjugate Gradient (PPCG) method of Gould, Hribar and Nocedal GHN02 () are well suited. Their convergence nevertheless depends highly on the choice of a good preconditioner RDW10 ().

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):

 ⎡⎢⎣2βM0−M0MKT−MK0⎤⎥⎦=⎡⎢ ⎢⎣I000I0−12βIKM−1I⎤⎥ ⎥⎦⋅⎡⎢⎣2βM000M000S⎤⎥⎦⋅⎡⎢ ⎢⎣I0−12βI0IM−1KT00I⎤⎥ ⎥⎦,S=−(12βM+KM−1KT). (6)

The hardest part in using this decomposition for solving a system is to compute the Schur complement of the (3,3) block and to solve the corresponding system with it. Since matrix 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 definitions 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 Green’s function of the Sturm–Liouville 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 find in the scientific literature. Among the most well-known are semiseparable GKK85 (), quasiseparable EG99a (), -matrices H99 (), -matrices HKS00 () and mosaic-skeleton matrices T96 (). In the current paper our main tool will be quasiseparable matrices VVM05 () having low rank partitions in their upper and lower parts. The formal definition is given below.

Definition 1 (Rank definition of a block quasiseparable matrix)

Let be a block matrix of block sizes then it is called block -quasiseparable if

 maxKrankA(K+1:N,1:K)≤rl,maxKrankA(1:K,K+1:N)≤ru,K=k∑i=1ni,N=n∑i=1ni,

where () is called lower (upper) order of quasiseparability.

In what follows we will call block -quasiseparable matrices simply quasiseparable for shortness.

A well-known example of a quasiseparable matrix is given by the discretized Green’s function of a one-dimensional differential equation. Consider a regular inhomogeneous Sturm-Liouville boundary value problem:

 (p(x)u′)′−q(x)u=f(x), (7) {α1u(a)+β1u(a)=0,α2u(b)+β2u(b)=0,

where functions , , are continuous on , on and for . It is a classical result in ordinary differential equations that the solution of (7) is given by

 u(x)=b∫ag(x,ξ)f(ξ)dξ,

where is the so-called Green’s function. In this case it has an explicit formula

 g(x,ξ)=1p(a)W(a)={u1(x)u2(ξ),a≤x≤ξ,u1(ξ)u2(x),ξ

with being Wronskian and and being two specific linearly independent solutions of (7). It is easy to see from (8) that discretized Green’s 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 parameters, where 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 . In this paper we will use the so-called generator representation (Definition 2 below) proposed by Eidelman and Gohberg EG99a (). There are several alternative parametrisations such as Givens-weight DV07 () and Hierarchical Semiseparable (HSS) CGP06 ().

Definition 2 (Generator definition of a block quasiseparable matrix)

Let be a block matrix of block sizes then it is called block -quasiseparable if it can be represented in the form

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣d1g1h2g1b2h3⋯g1b2…bn−1hnp2q1d2g2h3⋯g2b3…bn−1hnp3a2q1p3q2d3⋯g3b4…bn−1hn⋮⋮⋮⋱⋮pnan−1…a2q1pnan−1…a3q2pnan−1…a4q3⋯dn⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (9)

where parameters (called generators) , , , , , , are matrices of sizes as in the table below.

Orders of quasiseparability and are maxima over the corresponding sizes of generators:

 rl=max1≤k≤n−1rlk,ru=max1≤k≤n−1ruk.
Remark 1

Generators are not unique, there are infinitely many ways to represent the same quasiseparable matrix with different generators. For the relation between different sets of generators see EG05 ().

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 and of generators are directly related to ranks of submatrices in the lower and upper parts, respectively. Namely, if we let to be the block index: , then

 rankA(K+1:N,1:K)≤rlk,rankA(1:K,K+1:N)≤ruk. (10)

Moreover, for any quasiseparable matrix there exists a set of generators with sizes and that satisfy (10) with exact equalities (such generators are called minimal). For their existence and construction see EG05 ().

One remarkable property of quasiseparable matrices is that this class is closed under inversion EG99a (); SN04 (). For instance, inverse of a banded matrix is not banded but both matrices are quasiseparable. Due to the low-parametric representation in Definition 2 most of the algorithms with quasiseparable matrices have linear complexities. Table 2 lists these algorithms along with the corresponding references.

The ability to solve systems with quasiseparable matrices in operations 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 algorithm for a general block quasiseparable matrix in terms of the generators it is represented with. A restricted version of this algorithm applicable to diagonal plus semiseparable matrices (a subclass of quasiseparable matrices) was derived in GKK85 () and some formulae of the current algorithm are similar to the ones in the inversion algorithm given in EG99a (). 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 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 and factors. The theorem below makes this statement precise and, in addition, relates generators of an original matrix to generators of the factors.

Theorem 3.1

Let be a strongly regular block -quasiseparable matrix given by generators as in (9). Let be its block LU decomposition of the same block sizes. Then

1. Factors and are – and -quasiseparable. Moreover, and for all .

2. and are parametrized by the generators and , , , , , , , where are identity matrices of sizes and new parameters , and can be computed using the following algorithm:

Proof

Statement (i) of the theorem follows from statement (ii), so we only need to prove the latter part.

Denote, as usual, by the block index: and note that quasiseparable representation (9) implies the following recurrence relation between the blocks of :

 A=[A(1:K,1:K)GkHk+1Pk+1QkA(K+1:N,K+1:N)]; (11) Q1=q1,Qk=[akQk−1qk],k=2,…n−1; Pn=pn,Pk=[pk;Pk+1ak],k=n−1,…2; G1=g1,Gk=[Gk−1bk;gk],k=2,…n−1; Hn=hn,Hk=[hkbkHk+1],k=n−1,…2.

The proof will be constructed by induction. We will show that for each

 (12)

For we get from (12):

 d1=˜d1,P2Q1=P2˜Q1˜d1,G1H2=˜G1H2,⟸˜d1=d1,˜q1=q1˜d−11,˜g1=g1.

Let us introduce an auxiliary parameter . It is easy to show by using (11) that satisfies the recurrence relation

 f1=˜q1˜g1,fk=akfk−1bk+˜qk˜gk. (13)

Assume that (12) holds for some fixed , then it holds for if

 dk+1=[pk+1˜Qk1]⋅[˜Gkhk+1;˜dk+1], (14) Pk+2qk+1=Pk+2˜Qk+1⋅[˜Gkhk+1;˜dk+1], (15) gk+1Hk+2=[pk+1˜Qk1]⋅˜Gk+1Hk+2. (16)

The first equality simplifies

 dk+1=pk+1˜Qk˜Gkhk+1+˜dk+1=pk+1fkhk+1+˜dk+1.

The second equality (15) holds if

 qk+1=˜Qk+1⋅[˜Gkhk+1;˜dk+1]==[ak+1˜Qk˜qk+1]⋅[˜Gkhk+1;˜dk+1]=ak+1fkhk+1+˜qk+1˜dk+1.

Finally, the last equality (16) is true if

 gk+1=[pk+1˜Qk1]⋅˜Gk+1=[pk+1˜Qk1]⋅[˜Gkbk+1;˜gk+1]=pk+1fkbk+1+˜gk+1.

Matrix is invertible because, by our assumption, matrix is strongly regular. Hence, we conclude that (12) is true also for index if generators , and are those computed in Algorithm 1. The assertion of the theorem holds by induction.∎

4 Multilevel quasiseparable matrices.

In the previous section we have emphasized the relation between second-order ordinary differential equations and quasiseparable matrices. Let us now look at the problem in two dimensions that extends (7):

for , where with homogeneous Dirichlet boundary conditions. The standard five-point or nine-point discretization of this problem on a uniform grid leads to a system of linear algebraic equations of the form

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣A1B1C1A2B2C2⋱⋱⋱⋱Bm−1Cm−1Am⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣u1u2⋮⋮um⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣f1f2⋮⋮fm⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (17)

where we have assumed that are the discretized unknowns along the ’th column of the grid. In this case each of , , and is an matrix and the whole block matrix has size . Furthermore, each of the , , and 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 we get in the position occupied by the new block

 S2=A2−C1A−11B1.

Observe that even though all the individual matrices on the right-hand side of the expression are tridiagonal, is not, and hence is a dense (non-sparse) matrix. At the next step of Gaussian elimination we would use as the pivot block to eliminate . Now in the position occupied by the block we would get the matrix . Again, since is a dense matrix, in general will be dense, and therefore will also be a dense matrix. What this implies is that during LU factorization of the sparse matrix, we will produce fill-in quickly that causes us to compute inverses (and hence LU factorizations) of dense matrices. If we assume that these dense matrices have no structure, then we would need flops for that operation alone. Therefore, it follows that one would require at least flops to compute block LU factorization of the system matrix.

At first glance it seems that block LU factorization is not a wise approach as it does not use any kind of fill-in minimisation reordering. However, there is a rich structure in successive Schur complements that can be exploited to speed-up computations. In fact, it has been conjectured that if one looks at the off-diagonal blocks of these matrices (, , etc.), then their -rank (number of singular values not greater than ) is going to be small. This conjecture has been justified by the fact that, for example, can be viewed approximately (especially in the limit as becomes large) as a sub-block of the discretized Green’s function of the original PDE. It is known from the theory of elliptic PDEs (see, for instance, F95 ()) that under some mild constraints the Green’s function is smooth away from the diagonal singularity. This, in turn, implies that the numerical ranks of the off-diagonal blocks of are small. There have been several recent attempts to quantify this observation. It has been proved in CDGS10 () that -rank of the off-diagonal blocks of Schur complements in the LU decomposition of the regular 2D-Laplacian are bounded by

 1+8ln4(18ε).

This bound is not effective for any reasonable size of the problem because of the 4’th power in the logarithm (for instance for 1.0e-6 it gives 6.2e+5). However, numerical experiments with Laplacian discretized by the usual five-point stencil suggest much lower bound, see Figure 1.

Careful exploitation of the low-rank structure of Schur complements reduces the complexity of block LU factorization algorithm to linear . 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 complexity and there are steps altogether. First algorithms of this type were recently developed in XCGL09 (); M09 (). These algorithms, although efficient 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 defining the right structured representation of inverses of PDE matrices.

Definition 3 (Multilevel quasiseparable matrix)

Matrix is called -level quasiseparable if it admits the representation

 ⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣D1G1H2G1B2H3⋯G1B2…Bn−1HnP2Q1D2G2H3⋯G2B3…Bn−1HnP3A2Q1P3Q2D3⋯G3B4…Bn−1Hn⋮⋮⋮⋱⋮PnAn−1…A2Q1PnAn−1…A3Q2PnAn−1…A4Q3⋯Dn⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,

where each of the parameters , , , , , , is a block matrix of blocks being -level quasiseparable matrices. -level quasiseparable matrix is a usual quasiseparable matrix, see Definition 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 quasiseparable if blocks are themselves banded matrices. As an example, consider 2D Laplacian on the square discretized by Q1 finite elements:

 K=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣ABB⋱⋱⋱⋱BBA⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,A=13⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣−811⋱⋱⋱⋱11−8⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,B=13⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣111⋱⋱⋱⋱111⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (18)
Example 2 (Tensor product of quasiseparable matrices)

Matrix is 2-level quasiseparable if and are 1-level quasiseparable. Hence, 2D mass matrix on the square discretized by Q1 finite elements:

 M=T⊗T,T=16⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣411⋱⋱⋱⋱114⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦ (19)

is 2-level quasiseparable. Its inverse is also 2-level quasiseparable due to the tensor identity and properties of 1-level quasiseparable matrices.

Quasiseparable orders of multilevel quasiseparable matrices may differ at each level. For instance, it follows immediately from the properties of quasiseparable matrices that inverse of a 2-level quasiseparable matrix has a representation as in Definition 3 for some generators , , , , , , . However, there is no guarantee that entries of this generators are quasiseparable 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 one but as Figure 1 illustrates, generators of LU factors (and, therefore, inverse) have orders larger than one. However, if it is true that numerical orders (ranks) do not differ 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 arithmetic operations with 1-level quasiseparable matrices. Sum or difference of rank- and rank- matrices may have rank . Similarly, any arithmetic operation (except inversion) performed on quasiseparable matrices may add their orders in the worst case. Therefore, to use multilevel quasiseparable structure efficiently we need some kind of order compression algorithm. Such algorithm was proposed in EG05 (), we list it below for the convenience of the reader. This algorithm constitutes the essential part of the proposed method.

In exact arithmetic, new sizes of generators match ranks of the corresponding submatrices in the quasiseparable matrix (see EG05 () for the proof). In floating 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 submatrix. When -ranks are small, complexity of Algorithm 2 is and, hence, the following fact takes place.

Complexity of algorithms listed in Table 2 with 2-level quasiseparable matrices (Definition 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 multilevel quasiseparable matrices approach applied to solving PDEs and PDE-constrained optimization problems.

5 Numerical results.

We illustrate our new method with two different examples: 2D Laplace’s equation with Dirichlet boundary conditions and distributed control problem with a constraint in the form of 2D Poisson’s 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. Although 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 convection–diffusion 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 and consider the problem

 −∇2u=0,in Ω,u=ˆu|∂Ω,on ∂Ω, (20)

where

 ˆu=⎧⎨⎩sin(2πy)if x=0,−sin(2πy)if x=1,0otherwise.

Let us introduce a square reqular grid on of mesh size and discretize equation (20) on this grid by finite elements. Then equation (20) reduces to a matrix equation with matrix 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 different mesh sizes are presented in Tables 3 and 4. Note, that our solver was implemented for 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.

The analysis of results collected in Tables 3 and 4 reveals the linear dependence 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 preconditioned conjugate gradient method. Table 5 bellow illustrates the numerical performance of this approach. Results indicate that the PCG method preconditioned with the inexact LDU decomposition is twice as fast as more accurate LDU decomposition alone. It is enough to truncate quasiseparable order to 2–4 to force PCG to converge in less than 10 iterations.

Solution of problem in Example 3 obtained with PCG method with quasiseparable preconditioner is given in Figure 2 below.

Example 4

Let and consider the problem

 minu,f12∥u−ˆu∥2L2(Ω)+β∥f∥2L2(Ω), (21) s.t.−∇2u=f,in Ω,u=ˆu|∂Ω,on ∂Ω,

where

 ˆu={(2x−1)2(2y−1)2if (x,y)∈[0,12]2,0otherwise.

As it was already discussed in Section 2, problem (21) can be solved by the discretization by Q1 finite elements. Discretized system of equation, in this case, becomes

 ⎡⎢⎣2βM0−M0MKT−MK0⎤⎥⎦⎡⎢ ⎢⎣→f→u→λ⎤⎥ ⎥⎦=⎡⎢ ⎢⎣→0→b→d⎤⎥ ⎥⎦. (22)

After we eliminate variables and in equation (22) it reduces to

 (12βM+KM−1KT)→λ=→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 algorithm. One common approach is to drop either or in the equation and to solve system with the remaining matrix, thus, using it as a preconditioner (see, for instance, RDW10 ()). However, this preconditioner would perform badly if none of the terms is dominant. We propose a completely different approach. Matrices and 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 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 effective 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 and form matrix in quasiseparable matrices arithmetic.

2. Compute an approximate LDU decomposition of using Algorithm 1 and some order reduction tolerance.

3. Use PCG method to solve with approximate LDU decomposition as a preconditioner.

4. Exploit computed , and the factorization (6) to find and .

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 preconditioner 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 problem. Control function and obtained solution computed with the proposed algorithm for different values of the regularization parameter are presented in Figure 3.

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 parameters are needed to store a matrix. Moreover, arithmetic operations and matrix decompositions can be performed in floating-point operations. Multilevel quasiseparable matrices extend the applicability of well-known quasiseparable matrices EG99a () from one-dimensional to multidimensional problems. In particular, we have shown that multilevel quasiseparable structure is well-suited to the description of discretized elliptic operators in 2D.

To demonstrate the usefulness of the new class of matrices we considered distributed control problem with a constraint in the form of a partial differential equation. Such problems were introduced by Lions and Mitter in L71 (). 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 factorization. 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 BGL05 () for details. A large amount of work has been done on developing efficient preconditioners for such systems. In the current paper we have also proposed a new efficient mesh-independent preconditioner 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 differential operators in our examples. However, theory does not require symmetry in the matrix and it may very well be that multilevel quasiseparable matrices are applicable to convection–diffusion 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. Classifications of recurrence relations via subclasses of (H,m)-quasiseparable matrices. In Numerical Linear Algebra in Signals, Systems and Control, volume XV of Lecture Notes in Electrical Engineering, pages 23–53. Springer–Verlag, 2011.
• [2] T. Bella, V. Olshevsky, and P. Zhlobich. Signal flow graph approach to inversion of (H,m)-quasiseparable-Vandermonde matrices and new filter structures. Linear Algebra and its Applications, 432:2032–2051, 2010.
• [3] M. Benzi, G. Golub, and J. Liesen. Numerical solution of saddle point problems. Acta Numerica, 14:1–137, 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:137–152, 2004.
• [5] S. Chandrasekaran, P. Dewilde, M. Gu, and N. Somasunderam. On the numerical rank of the off-diagonal blocks of Schur complements of discretized elliptic PDEs. SIAM Journal on Matrix Analysis and Applications, 31:2261–2290, 2010.
• [6] S. Chandrasekaran, M. Gu, and T. Pals. A Fast Decomposition Solver for Hierarchically Semiseparable Representations. SIAM Journal on Matrix Analysis and Applications, 28(3):603–622, 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):1147–1170, 2007.
• [8] Y. Eidelman and I. Gohberg. Linear complexity inversion algorithms for a class of structured matrices. Integral Equations and Operator Theory, 35:28–52, 1999.
• [9] Y. Eidelman and I. Gohberg. On a new class of structured matrices. Integral Equations and Operator Theory, 34:293–324, 1999.
• [10] Y. Eidelman and I. Gohberg. A modification of the Dewilde-van der Veen method for inversion of finite structured matrices. Linear Algebra and its Applications, 343:419–450, 2002.
• [11] Y. Eidelman and I. Gohberg. On generators of quasiseparable finite block matrices. Calcolo, 42:187–214, 2005.
• [12] H.C. Elman, D.J. Silvester, and A.J. Wathen. Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Oxford University Press, USA, 2005.
• [13] G.B. Folland. Introduction to partial differential 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):780–804, 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 Scientific Computing, 23(4):1376–1395, 2002.
• [16] W. Hackbusch. A sparse matrix arithmetic based on -matrices. Part I: Introduction to -matrices. Computing, 62:89–108, 1999.
• [17] W. Hackbusch, B. Khoromskij, and S. Sauter. On -matrices. In H. Bungartz, R. Hoppe, and C. Zenger, editors, Lectures on Applied Mathematics, pages 9–29. Springer–Verlag, Berlin, 2000.
• [18] J.L. Lions and S.K. Mitter. Optimal control of systems governed by partial differential equations, volume 396. Springer Berlin, 1971.
• [19] P.G. Martinsson. A fast direct solver for a class of elliptic partial differential equations. Journal of Scientific Computing, 38(3):316–330, 2009.
• [20] C.C. Paige and M.A. Saunders. Solution of sparse indefinite systems of linear equations. SIAM Journal on Numerical Analysis, pages 617–629, 1975.
• [21] T. Rees, H.S. Dollar, and A.J. Wathen. Optimal solvers for PDE-constrained optimization. SIAM Journal on Scientific Computing, 32(1):271–298, 2010.
• [22] G. Strang and T. Nguyen. The interplay of ranks of submatrices. SIAM review, 46(4):637–646, 2004.
• [23] E. Tyrtyshnikov. Mosaic-skeleton approximations. Calcolo, 33(1):47–57, 1996.
• [24] R. Vandebril, M. Van Barel, and N. Mastronardi. A note on the representation and definition of semiseparable matrices. Numerical Linear Algebra with Applications, 12(8):839—858, 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):1382–1411, 2009.