This paper describes and compares some structure preserving techniques for the solution of linear discrete ill-posed problems with the t-product. A new randomized tensor singular value decomposition (R-tSVD) with a t-product is presented for low tubal rank tensor approximations. Regularization of linear inverse problems by truncated tensor eigenvalue decomposition (T-tEVD), truncated tSVD (T-tSVD), randomized T-tSVD (RT-tSVD), t-product Golub-Kahan bidiagonalization (tGKB) process, and t-product Lanczos (t-Lanczos) process are considered. A solution method that is based on reusing tensor Krylov subspaces generated by the tGKB process is described. The regularization parameter is the number of iterations required by each method. The discrepancy principle is used to determine this parameter. Solution methods that are based on truncated iterations are compared with solution methods that combine Tikhonov regularization with the tGKB and t-Lanczos processes. Computed examples illustrate the performance of these methods when applied to image and gray-scale video restorations. Our new RT-tSVD method is seen to require less CPU time and yields restorations of higher quality than the T-tSVD method.
We provide a global convergence proof of the recently proposed sequential homotopy method with an inexact Krylov--semismooth-Newton method employed as a local solver. The resulting method constitutes an active-set method in function space. After discretization, it allows for efficient application of Krylov-subspace methods. For a certain class of optimal control problems with PDE constraints, in which the control enters the Lagrangian only linearly, we propose and analyze an efficient, parallelizable, symmetric positive definite preconditioner based on a double Schur complement approach. We conclude with numerical results for a badly conditioned and highly nonlinear benchmark optimization problem with elliptic partial differential equations and control bounds. The resulting method is faster than using direct linear algebra for the 2D benchmark and allows for the parallel solution of large 3D problems.
This article considers the extension of two-grid $hp$-version discontinuous Galerkin finite element methods for the numerical approximation of second-order quasilinear elliptic boundary value problems of monotone type to the case when agglomerated polygonal/polyhedral meshes are employed for the coarse mesh approximation. We recall that within the two-grid setting, while it is necessary to solve a nonlinear problem on the coarse approximation space, only a linear problem must be computed on the original fine finite element space. In this article, the coarse space will be constructed by agglomerating elements from the original fine mesh. Here, we extend the existing a priori and a posteriori error analysis for the two-grid $hp$-version discontinuous Galerkin finite element method from 10.1007/s10915-012-9644-1 for coarse meshes consisting of standard element shapes to include arbitrarily agglomerated coarse grids. Moreover, we develop an $hp$-adaptive two-grid algorithm to adaptively design the fine and coarse finite element spaces; we stress that this is undertaken in a fully automatic manner, and hence can be viewed as blackbox solver. Numerical experiments are presented for two- and three-dimensional problems to demonstrate the computational performance of the proposed $hp$-adaptive two-grid method.
We present a practical algorithm to approximate the exponential of skew-Hermitian matrices up to round-off error based on an efficient computation of Chebyshev polynomials of matrices and the corresponding error analysis. It is based on Chebyshev polynomials of degrees 2, 4, 8, 12 and 18 which are computed with only 1, 2, 3, 4 and 5 matrix-matrix products, respectively. For problems of the form $\exp(-iA)$, with $A$ a real and symmetric matrix, an improved version is presented that computes the sine and cosine of $A$ with a reduced computational cost. The theoretical analysis, supported by numerical experiments, indicates that the new methods are more efficient than schemes based on rational Pad\'e approximants and Taylor polynomials for all tolerances and time interval lengths. The new procedure is particularly recommended to be used in conjunction with exponential integrators for the numerical time integration of the Schr\"odinger equation.
Mechanical cloaks are materials engineered to manipulate the elastic response around objects to make them indistinguishable from their homogeneous surroundings. Typically, methods based on material-parameter transformations are used to design optical, thermal and electric cloaks. However, they are not applicable in designing mechanical cloaks, since continuum-mechanics equations are not form-invariant under general coordinate transformations. As a result, existing design methods for mechanical cloaks have so far been limited to a narrow selection of voids with simple shapes. To address this challenge, we present a systematic, data-driven design approach to create mechanical cloaks composed of aperiodic metamaterials using a large pre-computed unit cell database. Our method is flexible to allow the design of cloaks with various boundary conditions, multiple loadings, different shapes and numbers of voids, and different homogeneous surroundings. It enables a concurrent optimization of both topology and properties distribution of the cloak. Compared to conventional fixed-shape solutions, this results in an overall better cloaking performance, and offers unparalleled versatility. Experimental measurements on 3D-printed structures further confirm the validity of the proposed approach. Our research illustrates the benefits of data-driven approaches in quickly responding to new design scenarios and resolving the computational challenge associated with multiscale designs of functional structures. It could be generalized to accommodate other applications that require heterogeneous property distribution, such as soft robots and implants design.
We propose a state redistribution method for high order discontinuous Galerkin methods on curvilinear embedded boundary grids. State redistribution relaxes the overly restrictive CFL condition that results from arbitrarily small cut cells and explicit time stepping. Thus, the scheme can take time steps that are proportional to the size of cells in the background grid. The discontinuous Galerkin scheme is stabilized by postprocessing the numerical solution after each stage or step of an explicit time stepping method. This is done by temporarily merging the small cells into larger, possibly overlapping neighborhoods using a special weighted inner product. Then, the numerical solution on the neighborhoods is returned to the base grid in a conservative fashion. The advantage of this approach is that it uses only basic mesh information that is already available in many cut cell codes and does not require complex geometric manipulations. Finally, we present a number of test problems that demonstrate the encouraging potential of this technique for applications on curvilinear embedded geometries. Numerical experiments reveal that our scheme converges with order $p+1$ in $L_1$ and between $p$ and $p+1$ in $L_\infty$ for problems with smooth solutions. We also demonstrate that state redistribution is capable of capturing shocks.
Many physical and mathematical models involve random fields in their input data. Examples are ordinary differential equations, partial differential equations and integro--differential equations with uncertainties in the coefficient functions described by random fields. They also play a dominant role in problems in machine learning. In this article, we do not assume to have knowledge of the moments or expansion terms of the random fields but we instead have only given discretized samples for them. We thus model some measurement process for this discrete information and then approximate the covariance operator of the original random field. Of course, the true covariance operator is of infinite rank and hence we can not assume to get an accurate approximation from a finite number of spatially discretized observations. On the other hand, smoothness of the true (unknown) covariance function results in effective low rank approximations to the true covariance operator. We derive explicit error estimates that involve the finite rank approximation error of the covariance operator, the Monte-Carlo-type errors for sampling in the stochastic domain and the numerical discretization error in the physical domain. This permits to give sufficient conditions on the three discretization parameters to guarantee that an error below a prescribed accuracy $\varepsilon$ is achieved.
In this paper, a class of arbitrarily high-order linear momentum-preserving and energy-preserving schemes are proposed, respectively, for solving the regularized long-wave equation. For the momentum-preserving scheme, the key idea is based on the extrapolation/prediction-correction technique and the symplectic Runge-Kutta method in time, together with the standard Fourier pseudo-spectral method in space. We show that the scheme is linear, high-order, unconditionally stable and preserves the discrete momentum of the system. For the energy-preserving scheme, it is mainly based on the energy quadratization approach and the analogous linearized strategy used in the construction of the linear momentum-preserving scheme. The proposed scheme is linear, high-order and can preserve a discrete quadratic energy exactly. Numerical results are addressed to demonstrate the accuracy and efficiency of the proposed scheme.
A novel class of high-order linearly implicit energy-preserving integrating factor Runge-Kutta methods are proposed for the nonlinear Schr\"odinger equation. Based on the idea of the scalar auxiliary variable approach, the original equation is first reformulated into an equivalent form which satisfies a quadratic energy. The spatial derivatives of the system are then approximated with the standard Fourier pseudo-spectral method. Subsequently, we apply the extrapolation technique/prediction-correction strategy to the nonlinear terms of the semi-discretized system and a linearized energy-conserving system is obtained. A fully discrete scheme is gained by further using the integrating factor Runge-Kutta method to the resulting system. We show that, under certain circumstances for the coefficients of a Runge-Kutta method, the proposed scheme can produce numerical solutions along which the modified energy is precisely conserved, as is the case with the analytical solution and is extremely efficient in the sense that only linear equations with constant coefficients need to be solved at every time step. Numerical results are addressed to demonstrate the remarkable superiority of the proposed schemes in comparison with other existing structure-preserving schemes.
We introduce a new class of preconditioners to enable flexible GMRES to find a least-squares solution, and potentially the pseudoinverse solution, of large-scale sparse, asymmetric, singular, and potentially inconsistent systems. We develop the preconditioners based on a new observation that generalized inverses (i.e., $\boldsymbol{A}^{g}\in\{\boldsymbol{G}\mid\boldsymbol{A}\boldsymbol{G}\boldsymbol{A}=\boldsymbol{A}\}$) enable the preconditioned Krylov subspaces to converge in a single step. We then compute an approximate generalized inverse (AGI) efficiently using a hybrid incomplete factorization (HIF), which combines multilevel incomplete LU with rank-revealing QR on its final Schur complement. We define the criteria of $\epsilon$-accuracy and stability of AGI to guarantee the convergence of preconditioned GMRES for consistent systems. For inconsistent systems, we fortify HIF with iterative refinement to obtain HIFIR, which allows accurate computations of the null-space vectors. By combining the two techniques, we then obtain a new solver, called PIPIT, for obtaining the pseudoinverse solutions for systems with low-dimensional null spaces. We demonstrate the robustness of HIF and HIFIR and show that they improve both accuracy and efficiency of the prior state of the art by orders of magnitude for systems with up to a million unknowns.
This manuscript presents an efficient solver for the linear system that arises from the Hierarchical Poincar\'e-Steklov (HPS) discretization of three dimensional variable coefficient Helmholtz problems. Previous work on the HPS method has tied it with a direct solver. This work is the first efficient iterative solver for the linear system that results from the HPS discretization. The solution technique utilizes GMRES coupled with an exact block-Jacobi preconditioner. The construction of the block-Jacobi preconditioner involves two nested local solves that are accelerated by local homogenization. The local nature of the discretization and preconditioner naturally yield matrix-free application of the linear system. A distributed memory implementation allows the solution technique to tackle problems approximately $50$ wavelengths in each direction requiring more than a billion unknowns to get approximately 7 digits of accuracy in less than an hour. Additional numerical results illustrate the performance of the solution technique.