Fast Parallel-in-Time Quasi-Boundary Value Methods for Backward Heat Conduction Problems

by   Jun Liu, et al.

In this paper we proposed two new quasi-boundary value methods for regularizing the ill-posed backward heat conduction problems. With a standard finite difference discretization in space and time, the obtained all-at-once nonsymmetric sparse linear systems have the desired block ω-circulant structure, which can be utilized to design an efficient parallel-in-time (PinT) direct solver that built upon an explicit FFT-based diagonalization of the time discretization matrix. Convergence analysis is presented to justify the optimal choice of the regularization parameter. Numerical examples are reported to validate our analysis and illustrate the superior computational efficiency of our proposed PinT methods.



There are no comments yet.


page 1

page 2

page 3

page 4


A Direct Parallel-in-Time Quasi-Boundary Value Method for Inverse Space-Dependent Source Problems

Inverse source problems arise often in real-world applications, such as ...

Numerical Analysis of Backward Subdiffusion Problems

The aim of this paper is to develop and analyze numerical schemes for ap...

Efficient Direct Space-Time Finite Element Solvers for Parabolic Initial-Boundary Value Problems in Anisotropic Sobolev Spaces

We consider a space-time variational formulation of parabolic initial-bo...

On the identification of source term in the heat equation from sparse data

We consider the recovery of a source term f(x,t)=p(x)q(t) for the nonhom...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Let and be an open and bounded domain in with a piecewise smooth boundary . We consider the classical backward heat conduction problem (BHCP) of reconstructing the initial data from the final time condition according to a homogeneous heat equation


This gives a severely ill-posed linear inverse problem that requires effective regularization techniques for stable and accurate numerical approximations (13; 22). Our proposed methods also work for other more general boundary conditions and spatial differential operators. In practice, the exact final condition is always unknown and we only have a noisy measurement , which is assumed to satisfy

with an estimated noise level


Due to the wide applications of BHCP, there appeared many different regularization methods in literature, such as quasi-boundary value methods (8; 11; 36; 18), optimal filtering method (37; 39), kernel–based method (1; 20; 41), Tikhonov regularization method (38; 43; 7; 12; 5), optimization method (21; 33), optimal control method (29; 24), total–variation regularization method (42), (Fourier) truncation regularization method (35; 32), fundamental solution method (6), Lie-group shooting method (4), meshless method (23), and homotopy analysis method (27). The majority of these works focuses on discussing the approximation accuracy and convergence rates under various different assumptions, while the computational efficiency of each method is inadequately investigated although it does often vary greatly among different approaches. Fast solvers for the resultant discretized regularized linear systems were rarely discussed.

The simple quasi-boundary value method (QBVM) was developed in (9), which improves the earlier quasi-reversibility regularization methods (25; 34). Based on QBVM, a modified quasi-boundary value method (MQBVM) was proposed in (11), which shows a better convergence rate than QBVM under stronger assumptions. Both QBVM and MQBVM regularize the ill-posed backward initial value problem through approximating it by a well-posed boundary value problem that depends on a regularization parameter . For the general non-homogeneous problems, several different extensions based on truncated series were discussed in (35; 40). In our recent work (29), we proposed a virtual optimal control formulation which can unify both QBVM and MQBVM by choosing some special Tikhonov regularization terms. Nevertheless, fast solvers for the nonsymmetric linear systems arising from QBVM and MQBVM were not developed yet in the literature, which are crucial to their applications to large-scale problems.

Meanwhile, with the advent of massively parallel computers, many efficient parallelizable numerical algorithms for solving evolutionary PDEs have been developed in the last few decades. Besides the achieved high parallelism in space, we have seen a lot of recent advances in various parallel-in-time (PinT) algorithms 111We refer to the website for a comprehensive list of references on various PinT algorithms. for solving forward time-dependent PDE problems (14). However, the application of such PinT algorithms to ill-posed backward heat conduction problems were rarely investigated in the literature, except in one short paper (10) about the parareal algorithm for a different parabolic inverse problem and another earlier paper (26) based on numerical (inverse) Laplace transform techniques in time. One obvious difficulty is how to address the underlying regularization treatment in the framework of PinT algorithms, which seems to be highly dependent on the problem structure. Inspired by several recent works (30; 31; 15; 28) on diagonalization-based PinT algorithms, we propose to redesign the existing quasi-boundary value methods in a structured way such that the diagonalization-based PinT direct solver can be directly employed, which can greatly speed up the quasi-boundary value methods without degrading their approximation accuracy.

The paper is organized as follows. In the next section we review two quasi-boundary value methods based on a standard finite difference scheme in space and time In Section 3, two new quasi-boundary value methods are proposed, where the derived block -circulant structured linear systems are then solved by a diagonalization-based PinT direct solver. Convergence analysis with the optimal choice of regularization parameter is discussed in Section 4. Two numerical examples are presented in Section 5 to demonstrate the promising efficiency of our proposed method and some conclusions are given in Section 6.

2 Two quasi-boundary value methods and their finite difference scheme

We now give a brief review of the quasi-boundary value methods, which were widely used in inverse due to its simplicity and effectiveness. We consider a 2D square domain with finite difference discretization in space and time, which can be easily adapted for 1D and 3D regular space domains. With suitable modification, the finite element discretization can also be used within our propose method for general spatial domains. We partition the time interval uniformly into with , and discretize the space domain uniformly into a uniform mesh with and , . For any function

, define the concatenated vector

with being the finite difference approximation of over all interior nodes at time . Let denotes the discretized Laplacian matrix with the five-point center finite difference and the homogeneous Dirichlet boundary conditions being enforced. Let and . Let and be identity matrices.

2.1 The Quasi-Boundary Value Method (QBVM)

Applying the QBVM (9) to BHCP (1), one needs to solve the following quasi-boundary value problem


where the regularization term with a parameter was introduced into the final condition. In noise-free setting the regularized solution converges to the original exact solution uniformly as goes to zero (9). With the backward Euler scheme in time and center finite difference in space, we obtain the full discretization scheme


which can be written into a nonsymmetric sparse linear system



Here the block of is from the regularization term . Notice that has a block Toeplitz or circulant structure except the first block row. Such a block Toeplitz or circulant structure is highly desirable for the development of fast system solvers. The key idea of our proposed PinT-QBVM is to redesign the regularization term such that a block -circulant structure is achieved, which can then be solved by a fast diagonalization-based PinT direct solver.

2.2 The Modified Quasi-Boundary Value Method (MQBVM)

Similarly, applying the MQBVM (11) to our BHCP (1), one needs to solve the quasi-boundary value problem


where a different regularization term was used. After applying the backward Euler scheme in time and center difference in space, we obtain the full discretization scheme


which can be further reformulated into a nonsymmetric sparse linear system




Here has the same entries as except the first block row. Again, the nonzero block of seems to more destructive to the anticipated block Toeplitz/circulant structure. In the next section, we will show this annoying block can be removed by manipulating the PDE itself under a reasonable regularity assumption.

The nonsymmetric sparse linear systems (5) and (9) are of dimension , which can be very costly to solve. The coupling of all time steps further increases the computational burden since we have to solve all time steps simultaneously in one-shot. A general sparse direct solver often has a complexity of order , which is prohibitive for 2D/3D problems with a fine mesh. In such large-scale cases, iterative solvers (e.g. Krylov subspace methods) are often the preferred choice, but effective and efficient preconditioners are required to achieve reasonably fast convergence rates. Moreover, for indefinite nonsymmetric systems, the rigorous convergence analysis of preconditioned iterative methods is a daunting task. In this paper, as the first step to apply PinT algorithms to ill-posed BHCP, we will concentrate on developing efficient diagonalization-based PinT direct solvers, which however can always be used as effective preconditioners within the framework of preconditioned iterative methods, especially for nonlinear problems. Fortunately, regularization is a versatile concept that can be exploited for generating better structured systems that are suitable for fast solvers. This simple yet powerful philosophy gains little attention in literature of BHCP.

3 Two new quasi-boundary value methods and their PinT implementation

Notice the matrices and from QBVM and MQBVM have a block Toeplitz structure except the first row due to the chosen regularization term. In this section, we redesign both QBVM and MQBVM to achieve a block -circulant structure of the system matrices, which consequently leads to a diagonalization-based PinT direct solver. We strive to retain their good regularization effects by keeping the original regularization terms of QBVM and MQBVM. At the same time, additional new terms are introduced to the final condition equation to twist the system structure.

3.1 A PinT quasi-boundary value method (PinT-QBVM)

To obtain the block -circulant structure that suitable for the diagonalization-based PinT algorithm, we propose to add an extra regularization term in QBVM as the following (compare with (2))


where is the same time step size to be used in its time discretization. In certain sense, we can review as an additional mesh-dependent regularization parameter, which would go to zero as the mesh is refined. Hence, we indeed only introduced a small perturbation term to the original regularization term of QBVM. In numerical simulations, we find that adding this extra regularization term indeed also leads to more accurate reconstruction.

Applying the backward Euler scheme in time and center finite difference scheme in space, upon dividing both sides of the final condition equation by we obtain the full scheme


which can be further reformulated into a nonsymmetric sparse linear system


where (with )

Here the matrix has the desired block -circulant structure, which is crucial to our PinT direct solver.

Assuming is differentiable in time at such that , then we can rewrite the final condition in (11) into

which can be interpreted as a weighted combination of both QBVM and MQBVM. This mathematical reformulation itself does not lead to a block -circulant structure, but it explains why PinT-QBVM may work better than QBVM.

3.2 A PinT modified quasi-boundary value method (PinT-MQBVM)

Similarly, we can redesign the MQBVM (11) by adding an additional regularization term (compared with (6))


Clearly, as the mesh step size gets smaller, we need to choose such that the new term is well controlled.

Again, assuming is differentiable at such that , upon dividing both sides by we can rewrite the final condition equation in (15) into

Applying the backward Euler scheme in time and center finite difference scheme in space, we obtain the full scheme


which can be further reformulated into a nonsymmetric sparse linear system


where (with )

Notice is the exactly same as except with a different value, and the first block of and is also different.

3.3 A diagonalization-based PinT direct solver

In contrast to and , both and have the same block -circulant structure, which will be utilized to design a PinT direct solver as explained below (based on ). For concise description, we will use the reshaping operations (17, p. 28): matrix-to-vector vec and vector-to-matrix mat. With the Kronecker product notation, we can write




Let (with and ) be the discrete Fourier matrix. For any given scalar number , define a diagonal matrix . It is well-known (2) that the -circulant matrix admits a explicit diagonalization , where , , and with being the first column of . With this explicit diagonalization , we can factorize into the product form

Hence, let , the solution vector can be computed via the following 3 steps:


where and denotes the -th column of . Here we have used the well-known Kronecker product property for any compatible matrices and . Clearly, the complex-shifted linear systems of size in Step-(b) can be computed in parallel since they are independent of each other. Notice that a different spatial operator would only affects the matrix in Step-(b). Moreover, due to , the matrix multiplication in Step-(a) and Step-(c) can be computed efficiently via FFT in time direction with complexity. In serial computation with sparse direct solver for Step-(b), the total complexity is , which is significantly lower than the complexity of sparse direct solver for the all-at-once system. Further parallel speedup can be achieved in parallel computing, we refer to (16; 3) for similar parallel results. Finally, we highlight that the circulant structure of is not essential to such a PinT direct solver. In fact, any efficient diagonalization of

with a well-conditioned eigenvector matrix

would be sufficient, but the computation of Step-(a) and (c) may need higher complexity unless has a special structure for fast computation.

4 Convergence analysis and the choice of regularization parameter

In this section we present the error estimates for both PinT-QBVM and PinT-MQBVM under suitable assumptions. Let and define a Hilbert function space equipped with the norm . Then admits a set of orthonormal eigenbasis in

, associated to a set of eigenvalues

such that with and . Give any , it has a series expansion with for all . Let be the compact contraction semi-group generated by . The noise-free un-regularized unique solution of the original BHCP (1) has a formal series expression


which is numerically unstable for computing , but very useful in convergence analysis. For a given , it was proved in (9, Lemma 1) that the original problem (1) has a unique (classical) solution if and only if Hence we assume that for a constant , which guarantees a unique solution.

Let ,and denotes the noise-free regularized solution of QBVM, MQBVM, PinT-QBVM, and PinT-MQBVM, respectively. Then we can easily verify the following explicit series representations


where is assumed to be fixed. Clearly with reduces to and with leads to .

We first discuss the stability estimates. Notice for any , for PinT-QBVM we have


which, upon setting , gives the known estimate (9) of QBVM: Similarly, for PinT-MQBVM there holds


Next, we estimate the regularization errors. Based on the series expansions of and , we have


Based on the series expansions of and , we have the error estimate


where the small in denominator will affect the choice of optimal parameter for PinT-MQBVM.

The following theorem shows PinT-QBVM and PinT-MQBVM can achieve the same convergence rate. Let and denotes the noisy regularized solution of each method with be replaced by noisy , respectively. Assume the noisy measurement be in such that holds for some .

Theorem 4.1.

Assume that for some constant and for some . Then for PinT-QBVM with the choice and PinT-MQBVM with the choice there hold


Combining the above estimates (4) and (4) for PinT-QBVM, there holds