Xiaozhe Hu

Xuan Zhao, Xiaozhe Hu, Wei Cai et al.

A robust and fast solver for the fractional differential equation (FDEs) involving the Riesz fractional derivative is developed using an adaptive finite element method on non-uniform meshes. It is based on the utilization of hierarchical matrices ($\mathcal{H}$-Matrices) for the representation of the stiffness matrix resulting from the finite element discretization of the FDEs. We employ a geometric multigrid method for the solution of the algebraic system of equations. We combine it with an adaptive algorithm based on a posteriori error estimation to deal with general-type singularities arising in the solution of the FDEs. Through various test examples we demonstrate the efficiency of the method and the high-accuracy of the numerical solution even in the presence of singularities. The proposed technique has been verified effectively through fundamental examples including Riesz, Left/Right Riemann-Liouville fractional derivative and, furthermore, it can be readily extended to more general fractional differential equations with different boundary conditions and low-order terms. To the best of our knowledge, there are currently no other methods for FDEs that resolve singularities accurately at linear complexity as the one we propose here.

Yicong Ma, Kaibo Hu, Xiaozhe Hu et al.

In this paper, we develop two classes of robust preconditioners for the structure-preserving discretization of the incompressible magnetohydrodynamics (MHD) system. By studying the well-posedness of the discrete system, we design block preconditioners for them and carry out rigorous analysis on their performance. We prove that such preconditioners are robust with respect to most physical and discretization parameters. In our proof, we improve the existing estimates of the block triangular preconditioners for saddle point problems by removing the scaling parameters, which are usually difficult to choose in practice. This new technique is not only applicable to the MHD system, but also to other problems. Moreover, we prove that Krylov iterative methods with our preconditioners preserve the divergence-free condition exactly, which complements the structure-preserving discretization. Another feature is that we can directly generalize this technique to other discretizations of the MHD system. We also present preliminary numerical results to support the theoretical results and demonstrate the robustness of the proposed preconditioners.

Xiaozhe Hu, Shuhong Wu, Xiao-Hui Wu et al.

We develop a simple algorithmic framework to solve large-scale symmetric positive definite linear systems. At its core, the framework relies on two components: (1) a norm-convergent iterative method (i.e. smoother) and (2) a preconditioner. The resulting preconditioner, which we refer to as a combined preconditioner, is much more robust and efficient than the iterative method and preconditioner when used in Krylov subspace methods. We prove that the combined preconditioner is positive definite and show estimates on the condition number of the preconditioned system. We combine an algebraic multigrid method and an incomplete factorization preconditioner to test the proposed framework on problems in petroleum reservoir simulation. Our numerical experiments demonstrate noticeable speed-up when we compare our combined method with the standalone algebraic multigrid method or the incomplete factorization preconditioner.

Xiaozhe Hu, Panayot S. Vassilevski, Jinchao Xu

The main purpose of this paper is to provide a comprehensive convergence analysis of nonlinear AMLI-cycle multigrid method for symmetric positive definite problems. Based on classical assumptions for approximation and smoothing properties, we show that the nonlinear AMLI-cycle MG method is uniformly convergent. Furthermore, under only the assumption that the smoother is convergent, we show that the nonlinear AMLI-cycle method is always better (or not worse) than the respective V-cycle MG method. Finally, numerical experiments are presented to illustrate the theoretical results.

Xiaozhe Hu, Carmen Rodrigo, Francisco J. Gaspar et al.

A stable finite element scheme that avoids pressure oscillations for a three-field Biot's model in poroelasticity is considered. The involved variables are the displacements, fluid flux (Darcy velocity), and the pore pressure, and they are discretized by using the lowest possible approximation order: Crouzeix-Raviart finite elements for the displacements, lowest order Raviart-Thomas-Nedelec elements for the Darcy velocity, and piecewise constant approximation for the pressure. Mass lumping technique is introduced for the Raviart-Thomas-Nedelec elements in order to eliminate the Darcy velocity and, therefore, reduce the computational cost. We show convergence of the discrete scheme which is implicit in time and use these types of elements in space with and without mass lumping. Finally, numerical experiments illustrate the convergence of the method and show its effectiveness to avoid spurious pressure oscillations when mass lumping for the Raviart-Thomas-Nedelec elements is used.

Paola F. Antonietti, Paul Houston, Xiaozhe Hu et al.

In this paper we analyze the convergence properties of two-level and W-cycle multigrid solvers for the numerical solution of the linear system of equations arising from hp-version symmetric interior penalty discontinuous Galerkin discretizations of second-order elliptic partial differential equations on polygonal/polyhedral meshes. We prove that the two-level method converges uniformly with respect to the granularity of the grid and the polynomial approximation degree p, provided that the number of smoothing steps, which depends on p, is chosen sufficiently large. An analogous result is obtained for the W-cycle multigrid algorithm, which is proved to be uniformly convergent with respect to the mesh size, the polynomial approximation degree, and the number of levels, provided the number of smoothing steps is chosen sufficiently large. Numerical experiments are presented which underpin the theoretical predictions; moreover, the proposed theoretical assumptions are not fully satisfied.

Nathaniel Trask, Martin Maxey, Xiaozhe Hu

A stable numerical solution of the steady Stokes problem requires compatibility between the choice of velocity and pressure approximation that has traditionally proven problematic for meshless methods. In this work, we present a discretization that couples a staggered scheme for pressure approximation with a divergence-free velocity reconstruction to obtain an adaptive, high-order, finite difference-like discretization that can be efficiently solved with conventional algebraic multigrid techniques. We use analytic benchmarks to demonstrate equal-order convergence for both velocity and pressure when solving problems with curvilinear geometries. In order to study problems in dense suspensions, we couple the solution for the flow to the equations of motion for freely suspended particles in an implicit monolithic scheme. The combination of high-order accuracy with fully-implicit schemes allows the accurate resolution of stiff lubrication forces directly from the solution of the Stokes problem without the need to introduce sub-grid lubrication models.

James Brannick, Fei Cao, Karsten Kahl et al.

In this paper, we consider a classical form of optimal algebraic multigrid (AMG) interpolation that directly minimizes the two-grid convergence rate and compare it with the so-called ideal form that minimizes a certain weak approximation property of the coarse space. We study compatible relaxation type estimates for the quality of the coarse grid and derive a new sharp measure using optimal interpolation that provides a guaranteed lower bound on the convergence rate of the resulting two-grid method for a given grid. In addition, we design a generalized bootstrap algebraic multigrid setup algorithm that computes a sparse approximation to the optimal interpolation matrix. We demonstrate numerically that the BAMG method with sparse interpolation matrix (and spanning multiple levels) outperforms the two-grid method with the standard ideal interpolation (a dense matrix) for various scalar diffusion problems with highly varying diffusion coefficient.

Xiaozhe Hu, Young-Ju Lee, Jinchao Xu et al.

In this paper, we consider the adaptive Eulerian--Lagrangian method (ELM) for linear convection-diffusion problems. Unlike the classical a posteriori error estimations, we estimate the temporal error along the characteristics and derive a new a posteriori error bound for ELM semi-discretization. With the help of this proposed error bound, we are able to show the optimal convergence rate of ELM for solutions with minimal regularity. Furthermore, by combining this error bound with a standard residual-type estimator for the spatial error, we obtain a posteriori error estimators for a fully discrete scheme. We present numerical tests to demonstrate the efficiency and robustness of our adaptive algorithm.

James Brannick, Yao Chen, Xiaozhe Hu et al.

We design and implement a parallel algebraic multigrid method for isotropic graph Laplacian problems on multicore Graphical Processing Units (GPUs). The proposed AMG method is based on the aggregation framework. The setup phase of the algorithm uses a parallel maximal independent set algorithm in forming aggregates and the resulting coarse level hierarchy is then used in a K-cycle iteration solve phase with a $\ell^1$-Jacobi smoother. Numerical tests of a parallel implementation of the method for graphics processors are presented to demonstrate its effectiveness.

Xiaozhe Hu, Junyuan Lin, Ludmil T. Zikatanov

We propose a path cover adaptive algebraic multigrid (PC-$α$AMG) method for solving linear systems of weighted graph Laplacians and can also be applied to discretized second order elliptic partial differential equations. The PC-$α$AMG is based on unsmoothed aggregation AMG (UA-AMG). To preserve the structure of smooth error down to the coarse levels, we approximate the level sets of the smooth error by first forming vertex-disjoint path cover with paths following the level sets. The aggregations are then formed by matching along the paths in the path cover. In such manner, we are able to build a multilevel structure at a low computational cost. The proposed PC-$α$AMG provides a mechanism to efficiently re-build the multilevel hierarchy during the iterations and leads to a fast nonlinear multilevel algorithm. Traditionally, UA-AMG requires more sophisticated cycling techniques, such as AMLI-cycle or K-cycle, but as our numerical results show, the PC-$α$AMG proposed here leads to nearly optimal standard V-cycle algorithm for solving linear systems with weighted graph Laplacians. Numerical experiments for some real world graph problems also demonstrate PC-$α$AMG's effectiveness and robustness, especially for ill-conditioned graphs.

Carmen Rodrigo, Francisco Gaspar, Xiaozhe Hu et al.

In this work we derive equivalence relations between mimetic finite difference schemes on simplicial grids and modified Nédélec-Raviart-Thomas finite element methods for model problems in $\mathbf{H}(\operatorname{\mathbf{curl}})$ and $H(\operatorname{div})$. This provides a simple and transparent way to analyze such mimetic finite difference discretizations using the well-known results from finite element theory. The finite element framework that we develop is also crucial for the design of efficient multigrid methods for mimetic finite difference discretizations, since it allows us to use canonical inter-grid transfer operators arising from the finite element framework. We provide special Local Fourier Analysis and numerical results to demonstrate the efficiency of such multigrid methods.

James H. Adler, Xiaozhe Hu, Ludmil T. Zikatanov

In this paper, we design robust and efficient linear solvers for the numerical approximation of solutions to Maxwell's equations with dissipative boundary conditions. We consider a structure-preserving finite-element approximation with standard Nedelec--Raviart--Thomas elements in space and a Crank--Nicolson scheme in time to approximate the electric and magnetic fields. We focus on two types of block preconditioners. The first type is based on the well-posedness results of the discrete problem. The second uses an exact block factorization of the linear system, for which the structure-preserving discretization yields sparse Schur complements. We prove robustness and optimality of these block preconditioners, and provide supporting numerical tests.

Bin Zheng, Luoping Chen, Xiaozhe Hu et al.

In this paper, we study fast iterative solvers for the solution of fourth order parabolic equations discretized by mixed finite element methods. We propose to use consistent mass matrix in the discretization and use lumped mass matrix to construct efficient preconditioners. We provide eigenvalue analysis for the preconditioned system and estimate the convergence rate of the preconditioned GMRes method. Furthermore, we show that these preconditioners only need to be solved inexactly by optimal multigrid algorithms. Our numerical examples indicate that the proposed preconditioners are very efficient and robust with respect to both discretization parameters and diffusion coefficients. We also investigate the performance of multigrid algorithms with either collective smoothers or distributive smoothers when solving the preconditioner systems.

Lu Wang, Xiaozhe Hu, Jonathan Cohen et al.

In this paper, we develop a new parallel auxiliary grid algebraic multigrid (AMG) method to leverage the power of graphic processing units (GPUs). In the construction of the hierarchical coarse grid, we use a simple and fixed coarsening procedure based on a region quadtree generated from an auxiliary grid. This allows us to explicitly control the sparsity patterns and operator complexities of the AMG solver. This feature provides (nearly) optimal load balancing and predictable communication patterns, which makes our new algorithm suitable for parallel computing, especially on GPU. We also design a parallel smoother based on the special coloring of the quadtree to accelerate the convergence rate and improve the parallel performance of this solver. Based on the CUDA toolkit [40], we implemented our new parallel auxiliary grid AMG method on GPU and the numerical results of this implementation demonstrate the efficiency of our new method. The results achieve an average speedup of over 4 on quasi-uniform grids and 2 on shape regular grids when compared to the AMG implementation in CUSP.

Ana Budiša, Xiaozhe Hu

In this paper, we are interested in an efficient numerical method for the mixed-dimensional approach to modeling single-phase flow in fractured porous media. The model introduces fractures and their intersections as lower-dimensional structures, and the mortar variable is used for flow coupling between the matrix and fractures. We consider a stable mixed finite element discretization of the problem, which results in a parameter-dependent linear system. For this, we develop block preconditioners based on the well-posedness of the discretization choice. The preconditioned iterative method demonstrates robustness with regards to discretization and physical parameters. The analytical results are verified on several examples of fracture network configurations, and notable results in reduction of number of iterations and computational time are obtained.

Xiaozhe Hu, John C. Urschel, Ludmil T. Zikatanov

Graph disaggregation is a technique used to address the high cost of computation for power law graphs on parallel processors. The few high-degree vertices are broken into multiple small-degree vertices, in order to allow for more efficient computation in parallel. In particular, we consider computations involving the graph Laplacian, which has significant applications, including diffusion mapping and graph partitioning, among others. We prove results regarding the spectral approximation of the Laplacian of the original graph by the Laplacian of the disaggregated graph. In addition, we construct an alternate disaggregation operator whose eigenvalues interlace those of the original Laplacian. Using this alternate operator, we construct a uniform preconditioner for the original graph Laplacian.

Qilong Zhai, Xiaozhe Hu, Ran Zhang

This paper proposes and analyzes a new weak Galerkin method for the eigenvalue problem by using the shifted-inverse power technique. A high order lower bound can be obtained at a relatively low cost via the proposed method. The error estimates for both eigenvalue and eigenfunction are provided and asymptotic lower bounds are shown as well under some conditions. Numerical examples are presented to validate the theoretical analysis.

James H. Adler, Samuel Hocking, Xiaozhe Hu et al.

Machine learning techniques have recently been of great interest for solving differential equations. Training these models is classically a data-fitting task, but knowledge of the expression of the differential equation can be used to supplement the training objective, leading to the development of physics-informed scientific machine learning. In this article, we focus on one class of models called nonlinear vector autoregression (NVAR) to solve ordinary differential equations (ODEs). Motivated by connections to numerical integration and physics-informed neural networks, we explicitly derive the physics-informed NVAR (piNVAR) which enforces the right-hand side of the underlying differential equation regardless of NVAR construction. Because NVAR and piNVAR completely share their learned parameters, we propose an augmented procedure to jointly train the two models. Then, using both data-driven and ODE-driven metrics, we evaluate the ability of the piNVAR model to predict solutions to various ODE systems, such as the undamped spring, a Lotka-Volterra predator-prey nonlinear model, and the chaotic Lorenz system.

Brooks Kinch, Xiaozhe Hu, Yilong Huang et al.

For autoregressive modeling of chaotic dynamical systems over long time horizons, the stability of both training and inference is a major challenge in building scientific foundation models. We present a hybrid technique in which an autoregressive transformer is embedded within a novel shooting-based mixed finite element scheme, exposing topological structure that enables provable stability. For forward problems, we prove preservation of discrete energies, while for training we prove uniform bounds on gradients, provably avoiding the exploding gradient problem. Combined with a vision transformer, this yields latent tokens admitting structure-preserving dynamics. We outperform modern foundation models with a $65\times$ reduction in model parameters and long-horizon forecasting of chaotic systems. A "mini-foundation" model of a fusion component shows that 12 simulations suffice to train a real-time surrogate, achieving a $9{,}000\times$ speedup over particle-in-cell simulation.

Li-Ping Liu, Ruiyuan Gu, Xiaozhe Hu

Polynomial functions have plenty of useful analytical properties, but they are rarely used as learning models because their function class is considered to be restricted. This work shows that when trained properly polynomial functions can be strong learning models. Particularly this work constructs polynomial feedforward neural networks using the product activation, a new activation function constructed from multiplications. The new neural network is a polynomial function and provides accurate control of its polynomial order. It can be trained by standard training techniques such as batch normalization and dropout. This new feedforward network covers several previous polynomial models as special cases. Compared with common feedforward neural networks, the polynomial feedforward network has closed-form calculations of a few interesting quantities, which are very useful in Bayesian learning. In a series of regression and classification tasks in the empirical study, the proposed model outperforms previous polynomial models.

Lenore Cowen, Kapil Devkota, Xiaozhe Hu et al.

Data-dependent metrics are powerful tools for learning the underlying structure of high-dimensional data. This article develops and analyzes a data-dependent metric known as diffusion state distance (DSD), which compares points using a data-driven diffusion process. Unlike related diffusion methods, DSDs incorporate information across time scales, which allows for the intrinsic data structure to be inferred in a parameter-free manner. This article develops a theory for DSD based on the multitemporal emergence of mesoscopic equilibria in the underlying diffusion process. New algorithms for denoising and dimension reduction with DSD are also proposed and analyzed. These approaches are based on a weighted spectral decomposition of the underlying diffusion process, and experiments on synthetic datasets and real biological networks illustrate the efficacy of the proposed algorithms in terms of both speed and accuracy. Throughout, comparisons with related methods are made, in order to illustrate the distinct advantages of DSD for datasets exhibiting multiscale structure.

Anuththari Gamage, Brian Rappaport, Shuchin Aeron et al.

Real-world data sets often provide multiple types of information about the same set of entities. This data is well represented by multi-view graphs, which consist of several distinct sets of edges over the same nodes. These can be used to analyze how entities interact from different viewpoints. Combining multiple views improves the quality of inferences drawn from the underlying data, which has increased interest in developing efficient multi-view graph embedding methods. We propose an algorithm, C-RSP, that generates a common (C) embedding of a multi-view graph using Randomized Shortest Paths (RSP). This algorithm generates a dissimilarity measure between nodes by minimizing the expected cost of a random walk between any two nodes across all views of a multi-view graph, in doing so encoding both the local and global structure of the graph. We test C-RSP on both real and synthetic data and show that it outperforms benchmark algorithms at embedding and clustering tasks while remaining computationally efficient.

Xiaozhe Hu, Jinchao Xu, Ludmil Zikatanov

In this paper, we consider the iterative method of subspace corrections with random ordering. We prove identities for the expected convergence rate, which can provide sharp estimates for the error reduction per iteration. We also study the fault-tolerant feature of the randomized successive subspace correction method by simply rejecting all the corrections when error occurs and show that the results iterative method converges with probability one. Moreover, we also provide sharp estimates on the expected convergence rate for the fault-tolerant, randomized, subspace correction method.

Carmen Rodrigo, Xiaozhe Hu, Peter Ohm et al.

In this work, we consider the popular P1-RT0-P0 discretization of the three-field formulation of Biot's consolidation problem. Since this finite-element formulation does not satisfy an inf-sup condition uniformly with respect to the physical parameters, several issues arise in numerical simulations. For example, when the permeability is small with respect to the mesh size, volumetric locking may occur. Thus, we propose a stabilization technique that enriches the piecewise linear finite-element space of the displacement with the span of edge/face bubble functions. We show that for Biot's model this does give rise to discretizations that are uniformly stable with respect to the physical parameters. We also propose a perturbation of the bilinear form, which allows for local elimination of the bubble functions and provides a uniformly stable scheme with the same number of degrees of freedom as the classical P1-RT0-P0 approach. We prove optimal stability and error estimates for this discretization. Finally, we show that this scheme can also be successfully applied to Stokes' equations, yielding a discrete problem with optimal approximation properties and with minimum number of degrees of freedom (equivalent to a P1-P0 discretization). Numerical tests confirm the theory for both poroelastic and Stokes' test problems.

James H. Adler, Francisco J. Gaspar, Xiaozhe Hu et al.

In this paper, we design robust and efficient block preconditioners for the two-field formulation of Biot's consolidation model, where stabilized finite-element discretizations are used. The proposed block preconditioners are based on the well-posedness of the discrete linear systems. Block diagonal (norm-equivalent) and block triangular preconditioners are developed, and we prove that these methods are robust with respect to both physical and discretization parameters. Numerical results are presented to support the theoretical results.

Carmen Rodrigo, Francisco Gaspar, Xiaozhe Hu et al.

We consider finite element discretizations of the Biot's consolidation model in poroelasticity with MINI and stabilized P1-P1 elements. We analyze the convergence of the fully discrete model based on spatial discretization with these types of finite elements and implicit Euler method in time. We also address the issue related to the presence of non-physical oscillations in the pressure approximation for low permeabilities and/or small time steps. We show that even in 1D a Stokes-stable finite element pair fails to provide a monotone discretization for the pressure in such regimes. We then introduce a stabilization term which removes the oscillations. We present numerical results confirming the monotone behavior of the stabilized schemes.

James Brannick, Xiaozhe Hu, Carmen Rodrigo et al.

We focus on the study of multigrid methods with aggressive coarsening and polynomial smoothers for the solution of the linear systems corresponding to finite difference/element discretizations of the Laplace equation. Using local Fourier analysis we determine automatically the optimal values for the parameters involved in defining the polynomial smoothers and achieve fast convergence of cycles with aggressive coarsening. We also present numerical tests supporting the theoretical results and the heuristic ideas. The methods we introduce are highly parallelizable and efficient multigrid algorithms on structured and semi-structured grids in two and three spatial dimensions.

John C. Urschel, Xiaozhe Hu, Jinchao Xu et al.

In this paper, we develop a cascadic multigrid algorithm for fast computation of the Fiedler vector of a graph Laplacian, namely, the eigenvector corresponding to the second smallest eigenvalue. This vector has been found to have applications in fields such as graph partitioning and graph drawing. The algorithm is a purely algebraic approach based on a heavy edge coarsening scheme and pointwise smoothing for refinement. To gain theoretical insight, we also consider the related cascadic multigrid method in the geometric setting for elliptic eigenvalue problems and show its uniform convergence under certain assumptions. Numerical tests are presented for computing the Fiedler vector of several practical graphs, and numerical results show the efficiency and optimality of our proposed cascadic multigrid algorithm.