Schur Complement Reduction Research Articles

An ADMM-Newton-CNN numerical approach to a TV model for identifying discontinuous diffusion coefficients in elliptic equations: convex case with gradient observations

Identifying the discontinuous diffusion coefficient in an elliptic equation with observation data of the gradient of the solution is an important nonlinear and ill-posed inverse problem. Models with total variational (TV) regularization have been widely studied for this problem, while the theoretically required nonsmoothness property of the TV regularization and the hidden convexity of the models are usually sacrificed when numerical schemes are considered in the literature. In this paper, we show that the favorable nonsmoothness and convexity properties can be entirely kept if the well-known alternating direction method of multipliers (ADMM) is applied to the TV-regularized models, hence it is meaningful to consider designing numerical schemes based on the ADMM. Moreover, we show that one of the ADMM subproblems can be well solved by the active-set Newton method along with the Schur complement reduction method, and the other one can be efficiently solved by the deep convolutional neural network (CNN). The resulting ADMM-Newton-CNN approach is demonstrated to be easily implementable and very efficient even for higher-dimensional spaces with fine mesh discretization.

Efficient solution techniques for two-phase flow in heterogeneous porous media using exact Jacobians

Two efficient and scalable numerical solution methods will be compared using exact Jacobians to solve the fully coupled Newton systems arising during fully implicit discretization of the equations for two-phase flow in porous media. These methods use algebraic multigrid (AMG) to solve the linear systems in every Newton step. The algebraic multigrid methods rely on (i) a Schur Complement Reduction (SCR-AMG) and (ii) a Constrained Pressure Residual method (CPR-AMG) to decouple elliptic and hyperbolic contributions. Both methods employ automatic differentiation (AD) to calculate exact Jacobians and are compared with finite difference (FD) approximations of the Jacobian. The superiority of AD is shown by several numerical test cases from the field of CO2 geo-sequestration comprising two- and three-dimensional examples. A weak scaling test on JUQUEEN, a BlueGene/Q supercomputer, demonstrates the efficiency and scalability of both methods. To achieve maximum comparability and reproducibility, the Portable Extensible Toolkit for Scientific Computation (PETSc) is used as framework for the comparison of all solvers.

A robust and efficient iterative method for hyper-elastodynamics with nested block preconditioning

We develop a robust and efficient iterative method for hyper-elastodynamics based on a novel continuum formulation recently developed in [1]. The numerical scheme is constructed based on the variational multiscale formulation and the generalized-α method. Within the nonlinear solution procedure, a block factorization is performed for the consistent tangent matrix to decouple the kinematics from the balance laws. Within the linear solution procedure, another block factorization is performed to decouple the mass balance equation from the linear momentum balance equations. A nested block preconditioning technique is proposed to combine the Schur complement reduction approach with the fully coupled approach. This preconditioning technique, together with the Krylov subspace method, constitutes a novel iterative method for solving hyper-elastodynamics. We demonstrate the efficacy of the proposed preconditioning technique by comparing with the SIMPLE preconditioner and the one-level domain decomposition preconditioner. Two representative examples are studied: the compression of an isotropic hyperelastic cube and the tensile test of a fully-incompressible anisotropic hyperelastic arterial wall model. The robustness with respect to material properties and the parallel performance of the preconditioner are examined.

English

Two-by-two block matrices of special form with square matrix blocks arise in important applications, such as in optimal control of partial differential equations and in high order time integration methods. Two solution methods involving very efficient preconditioned matrices, one based on a Schur complement reduction of the given system and one based on a transformation matrix with a perturbation of one of the given matrix blocks are presented. The first method involves an additional inner solution with the pivot matrix block but gives a very tight condition number bound when applied for a time integration method. The second method does not involve this matrix block but only inner solutions with a linear combination of the pivot block and the off-diagonal matrix blocks. Both the methods give small condition number bounds that hold uniformly in all parameters involved in the problem, i.e. are fully robust. The paper presents shorter proofs, extended and new results compared to earlier publications.

Robust operator estimates and the application to substructuring methods for first-order systems

We discuss a family of discontinuous Petrov–Galerkin (DPG) schemes for quite general partial differential operators. The starting point of our analysis is the DPG method introduced by [Demkowicz et al. , SIAM J. Numer. Anal. 49 (2011) 1788–1809; Zitelli et al. , J. Comput. Phys. 230 (2011) 2406–2432]. This discretization results in a sparse positive definite linear algebraic system which can be obtained from a saddle point problem by an element-wise Schur complement reduction applied to the test space. Here, we show that the abstract framework of saddle point problems and domain decomposition techniques provide stability and a priori estimates. To obtain efficient numerical algorithms, we use a second Schur complement reduction applied to the trial space. This restricts the degrees of freedom to the skeleton. We construct a preconditioner for the skeleton problem, and the efficiency of the discretization and the solution method is demonstrated by numerical examples.

Projected Krylov Methods for Solving Non-Symmetric Two-by-Two Block Linear Systems Arising from Fictitious Domain Formulations

The paper deals with the solution of large non-symmetric two-by-two block linear systems with a singular leading submatrix. Our algorithm consists of two levels. The outer level combines the Schur complement reduction with the orthogonal projectors that leads to the linear equation on subspaces. To solve this equation, we use a Krylov-type method representing the inner level of the algorithm. We propose a general technique how to get from the standard Krylov methods their projected variants generating iterations on subspaces. Then we derive the projected GMRES. The efficiency of our approach is illustrated by examples arising from the combination of the fictitious domain and FETI method.

Application of the iterative solution method with Schur complement reduction to mixed finite elements based in a tetrahedral discretization

The iterative solution method for mixed finite element methods is applied to a 3-D domain partitioned with tetrahedral elements. For the particular discretization technique of first partitioning the domain with hexahedral cells, and then subsequently partitioning cells with five tetrahedral elements, a Schur complement decomposition is devised wherein the actual number of equations solved is reduced by 80%. Although this Schur complement reduction requires a fair amount of computational overhead, its application within the iterative solution method can reduce overall solution time by about 44%, depending on closure criterion and other factors.

Minimum Rank of Graphs with Loops

A loop graph $\mf G$ is a finite undirected graph that allows loops but does not allow multiple edges. The set $\sym(\lG)$ of real symmetric matrices associated with a loop graph $\lG$ of order $n$ is the set of symmetric matrices $A=[a_{ij}]\in\Rnn$ such that $a_{ij}\ne 0$ if and only if $ij\in E(\lG)$. The minimum (maximum) rank of a loop graph is the minimum (maximum) of the ranks of the matrices in $\sym(\lG)$. We characterize loop graphs having minimum rank at most two (by forbidden induced subgraphs and graph complements) and loop graphs having minimum rank equal to the order of the graph. A Schur complement reduction technique is used to determine the minimum ranks of cycles with various loop configurations; we also determine the minimum ranks of complete graphs and paths with various configurations of loops. Unlike simple graphs, loop graphs can have maximum rank less than the order of the graph. We present some results on maximum rank and which ranks between minimum and maximum can be realized. Interesting open questions remain.

Robust preconditioners for the high-contrast Stokes equation

We consider the Stokes equation with high-contrast viscosity coefficients. We construct a preconditioner that is robust with respect to contrast size and mesh size simultaneously based on the preconditioner proposed by Aksoylu et al. (2008). We examine the performance of our preconditioner against multigrid and provide a comparative study reflecting the effect of the underlying discretization and the aspect ratio of the mesh. We address the rigorous justification of the solver methods, p-Uzawa and p-Minres, used in Aksoylu and Unlu (2013), and compare the results with additional solver method, Schur complement reduction (SCR). We observe that our preconditioner is only contrast size robust under the p-SCR solver. The inexact p-Uzawa solver remains to be the best choice for the most effective performance of our preconditioner as we observe contrast size and mesh size robustness simultaneously in this case. As the contrast size grows asymptotically, we prove and numerically demonstrate that the inexact p-Uzawa solver converges to the exact one. Finally, we show that our preconditioner is contrast size and mesh size robust under p-Minres when the Schur complement solve is accurate enough.

On large-scale generalized inverses in solving two-by-two block linear systems

The goal is to analyze a role of generalized inverses in the solution of two-by-two block linear systems by means of an algorithm combining the Schur complement reduction with the null-space method. The efficient procedures for computing actions of generalized inversions to rectangular large-scale matrices are presented with a special attention to the Moore–Penrose inverse. A regularization technique for matrices with the known null-space basis is introduced that enables to get generalized inverses without necessity to recognize zero pivots. The theoretical results are confirmed by numerical examples.

Preconditioned iterative methods for Stokes flow problems arising in computational geodynamics

The numerical solution of a discrete Stokes flow problem is a common component of many computational geodynamic models. We compare using the Schur complement reduction (SCR) and a fully coupled (FC) approach to solve geodynamic problems which are characterized by incompressible fluids containing large variations in viscosity. The scalability and robustness of these methods is assessed by examining their sensitivity to the discretization parameter h and the viscosity contrast Δ η , respectively. We demonstrate that the scaled BFBt preconditioner [Elman, H.C., Howle, V.E., Shadid, J., Shuttleworth, R., Tuminaro, R., 2006. Block preconditioners based on approximate commutators. SIAM J. Sci. Comput. 27 (5), 1651–1668] can be an effective preconditioner for S = G T K − 1 G , even when the fluid viscosity varies in space, provided an appropriate scaling is used. The performance of this preconditioner is strongly linked to the scaling employed. Through a number of numerical experiments, we demonstrate that our new scaling strategy ensures that the scaled BFBt preconditioner remains effective in scenarios where the viscosity variations are either locally smooth or discontinuous. Of the solvers considered, applying FGMRES to the pressure Schur complement in combination with BFBt and our choice of scaling matrices was the most scalable and robust approach.

Maximum Attainable Accuracy of Inexact Saddle Point Solvers

In this paper we study numerical behavior of several iterative Krylov subspace solvers applied to the solution of large-scale saddle point problems. Two main representatives of segregated solution approach are analyzed: the Schur complement reduction method based on the elimination of primary unknowns and the null-space projection method, which relies on a basis for the subspace described by the constraints. We show that the choice of the back-substitution formula may considerably influence the maximum attainable accuracy of approximate solutions computed in finite precision arithmetic.

Projected Schur complement method for solving non‐symmetric systems arising from a smooth fictitious domain approach

AbstractThis paper deals with a fast method for solving large‐scale algebraic saddle‐point systems arising from fictitious domain formulations of elliptic boundary value problems. A new variant of the fictitious domain approach is analyzed. Boundary conditions are enforced by control variables introduced on an auxiliary boundary located outside the original domain. This approach has a significantly higher convergence rate; however, the algebraic systems resulting from finite element discretizations are typically non‐symmetric. The presented method is based on the Schur complement reduction. If the stiffness matrix is singular, the reduced system can be formulated again as another saddle‐point problem. Its modification by orthogonal projectors leads to an equation that can be efficiently solved using a projected Krylov subspace method for non‐symmetric operators. For this purpose, the projected variant of the BiCGSTAB algorithm is derived from the non‐projected one. The behavior of the method is illustrated by examples, in which the BiCGSTAB iterations are accelerated by a multigrid strategy. Copyright © 2007 John Wiley & Sons, Ltd.

Limiting accuracy of segregated solution methods for nonsymmetric saddle point problems

Nonsymmetric saddle point problems arise in a wide variety of applications in computational science and engineering. The aim of this paper is to discuss the numerical behavior of several nonsymmetric iterative methods applied for solving the saddle point systems via the Schur complement reduction or the null-space projection approach. Krylov subspace methods often produce the iterates which fluctuate rather strongly. Here we address the question whether large intermediate approximate solutions reduce the final accuracy of these two-level (inner–outer) iteration algorithms. We extend our previous analysis obtained for symmetric saddle point problems and distinguish between three mathematically equivalent back-substitution schemes which lead to a different numerical behavior when applied in finite precision arithmetic. Theoretical results are then illustrated on a simple model example.

Schur complement reduction in the mixed-hybrid approximation of Darcy's law: rounding error analysis

Mixed-hybrid finite element approximation of the potential fluid flow problem leads to the solution of a large symmetric indefinite system for the velocity and potential head vector components. Such discretization gives rise to a very accurate approximation of the continuity equation in every element, and for low-order discretizations, the structural properties of the discrete matrix blocks allow cheap block elimination of the positive-definite diagonal block and subsequent reduction to the Schur complement system for the pressure and Lagrangian vector components. This system is then frequently solved by the iterative conjugate gradient-type method. Whereas this approach is well known, considerably less attention has been paid to the numerical stability aspects of such transformation. It was shown in [5] that block LU factorization can be unstable even when the system matrix is symmetric positive definite. In this paper we examine this type of conditional stability for a particular application in the underground water flow modelling. We show that the actual error of the computed approximate solution depends not only on the user-defined tolerance in the conjugate gradient process but also on the spectral properties of the corresponding matrix blocks eliminated during the Schur complement reduction. It is often observed that although the backward error of the approximate solution in the iterative part is reduced to the level of machine accuracy, the total residual norm after the back-substitution process remains at certain accuracy level. We give a bound for this maximal attainable accuracy and illustrate our theoretical results on a model example.

An asymptotically optimal Schur complement reduction for the Stokes equation

In this paper we develop an efficient Schur complement method for solving the 2D Stokes equation. As a basic algorithm, we apply a decomposition approach with respect to the trace of the pressure. The alternative stream function-vorticity reduction is also discussed. The original problem is reduced to solving the equivalent boundary (interface) equation with symmetric and positive definite operator in the appropriate trace space. We apply a mixed finite element approximation to the interface operator by $P_1$ iso $P_2/P_1$ triangular elements and prove the optimal error estimates in the presence of stabilizing bubble functions. The norm equivalences for the corresponding discrete operators are established. Then we propose an asymptotically optimal compression technique for the related stiffness matrix (in the absence of bubble functions) providing a sparse factorized approximation to the Schur complement. In this case, the algorithm is shown to have an optimal complexity of the order $O(N \log^q N)$ , q = 2 or q = 3, depending on the geometry, where N is the number of degrees of freedom on the interface. In the presence of bubble functions, our method has the complexity $O(N^2 \log N)$ arithmetical operations. The Schur complement interface equation is resolved by the PCG iterations with an optimal preconditioner.

Fast computations with the harmonic Poincaré–Steklov operators on nested refined meshes

In this paper we develop asymptotically optimal algorithms for fast computations with the discrete harmonic Poincare–Steklov operators (Dirichlet–Neumann mapping) for interior and exterior problems in the presence of a nested mesh refinement. Our approach is based on the multilevel interface solver applied to the Schur complement reduction onto the nested refined interface associated with a nonmatching decomposition of a polygon by rectangular substructures. This paper extends methods from Khoromskij and Prossdorf (1995), where the finite element approximations of interior problems on quasi‐uniform grids have been considered. For both interior and exterior problems, the matrix–vector multiplication with the compressed Schur complement matrix on the interface is shown to have a complexity of the order O(Nr log3Nu), where Nr = O((1 + pr) Nu) is the number of degrees of freedom on the polygonal boundary under consideration, Nu is the boundary dimension of a finest quasi‐uniform level and pr ⩾ 0 defines the refinement depth. The corresponding memory needs are estimated by O(Nr logqNu), where q = 2 or q = 3 in the case of interior and exterior problems, respectively.