Motivated by a wide range of real-world problems whose solutions exhibit boundary and interior layers, the numerical analysis of discretizations of singularly perturbed differential equations is an established sub-discipline within the study of the numerical approximation of solutions to differential equations. Consequently, much is known about how to accurately and stably discretize such equations on \textit{a priori} adapted meshes, in order to properly resolve the layer structure present in their continuum solutions. However, despite being a key step in the numerical simulation process, much less is known about the efficient and accurate solution of the linear systems of equations corresponding to these discretizations. In this paper, we discuss problems associated with the application of direct solvers to these discretizations, and we propose a preconditioning strategy that is tuned to the matrix structure induced by using layer-adapted meshes for convection-diffusion equations, proving a strong condition-number bound on the preconditioned system in one spatial dimension, and a weaker bound in two spatial dimensions. Numerical results confirm the efficiency of the resulting preconditioners in one and two dimensions, with time-to-solution of less than one second for representative problems on $1024\times 1024$ meshes and up to $40\times$ speedup over standard sparse direct solvers.
Orthogonal polynomials of several variables have a vector-valued three-term recurrence relation, much like the corresponding one-dimensional relation. This relation requires only knowledge of certain recurrence matrices, and allows simple and stable evaluation of multivariate orthogonal polynomials. In the univariate case, various algorithms can evaluate the recurrence coefficients given the ability to compute polynomial moments, but such a procedure is absent in multiple dimensions. We present a new Multivariate Stieltjes (MS) algorithm that fills this gap in the multivariate case, allowing computation of recurrence matrices assuming moments are available. The algorithm is essentially explicit in two and three dimensions, but requires the numerical solution to a non-convex problem in more than three dimensions. Compared to direct Gram-Schmidt-type orthogonalization, we demonstrate on several examples in up to three dimensions that the MS algorithm is far more stable, and allows accurate computation of orthogonal bases in the multivariate setting, in contrast to direct orthogonalization approaches.
Sparse superposition codes were originally proposed as a capacity-achieving communication scheme over the gaussian channel, whose coding matrices were made of i.i.d. gaussian entries.We extend this coding scheme to more generic ensembles of rotational invariant coding matrices with arbitrary spectrum, which include the gaussian ensemble as a special case. We further introduce and analyse a decoder based on vector approximate message-passing (VAMP).Our main findings, based on both a standard replica symmetric potential theory and state evolution analysis, are the superiority of certain structured ensembles of coding matrices (such as partial row-orthogonal) when compared to i.i.d. matrices, as well as a spectrum-independent upper bound on VAMP's threshold. Most importantly, we derive a simple "spectral criterion " for the scheme to be at the same time capacity-achieving while having the best possible algorithmic threshold, in the "large section size" asymptotic limit. Our results therefore provide practical design principles for the coding matrices in this promising communication scheme.
Imposing orthogonality on the layers of neural networks is known to facilitate the learning by limiting the exploding/vanishing of the gradient; decorrelate the features; improve the robustness. This paper studies theoretical properties of orthogonal convolutional layers. We establish necessary and sufficient conditions on the layer architecture guaranteeing the existence of an orthogonal convolutional transform. The conditions prove that orthogonal convolutional transforms exist for almost all architectures used in practice for 'circular' padding.We also exhibit limitations with 'valid' boundary condition and 'same' boundary condition with zero padding. Recently, a regularization term imposing the orthogonality of convolutional layers has been proposed, and impressive empirical results have been obtained in different applications (Wang et al. 2020).The second motivation of the present paper is to specify the theory behind this.We make the link between this regularization term and orthogonality measures. In doing so, we show that this regularization strategy is stable with respect to numerical and optimization errors and that, in the presence of small errors and when the size of the signal/image is large, the convolutional layers remain close to isometric.The theoretical results are confirmed with experiments, the landscape of the regularization term is studied and the regularization strategy is validated on real datasets. Altogether, the study guarantees that the regularization with L_{orth} (Wang et al. 2020) is an efficient, flexible and stable numerical strategy to learn orthogonal convolutional layers.
Given an array of distinct integers $A[1\ldots n]$, the Range Minimum Query (RMQ) problem requires us to construct a data structure from $A$, supporting the RMQ query: given an interval $[a,b]\subseteq[1,n]$, return the index of the minimum element in subarray $A[a\ldots b]$, i.e. return $\text{argmin}_{i\in[a,b]}A[i]$. The fundamental problem has a long history. The textbook solution which uses $O(n)$ words of space and $O(1)$ time by Gabow, Bentley, Tarjan (STOC 1984) and Harel, Tarjan (SICOMP 1984) dates back to 1980s. The state-of-the-art solution is presented by Fischer, Heun (SICOMP 2011) and Navarro, Sadakane (TALG 2014). The solution uses $2n-1.5\log n+n/\left(\frac{\log n}{t}\right)^t+\tilde{O}(n^{3/4})$ bits of space and $O(t)$ query time, where the additive $\tilde{O}(n^{3/4})$ is a pre-computed lookup table used in the RAM model, assuming the word-size is $\Theta(\log n)$ bits. On the other hand, the only known lower bound is proved by Liu and Yu (STOC 2020). They show that any data structure which solves RMQ in $t$ query time must use $2n-1.5\log n+n/(\log n)^{O(t^2\log^2t)}$ bits of space, assuming the word-size is $\Theta(\log n)$ bits. In this paper, we prove nearly tight lower bound for this problem. We show that, for any data structure which solves RMQ in $t$ query time, $2n-1.5\log n+n/(\log n)^{O(t\log^2t)}$ bits of space is necessary in the cell-probe model with word-size $\Theta(\log n)$ bits. We emphasize that, in terms of time complexity, our lower bound is tight up to a polylogarithmic factor.
In our previous work [SIAM J. Sci. Comput. 43(3) (2021) B784-B810], an accurate hyper-singular boundary integral equation method for dynamic poroelasticity in two dimensions has been developed. This work is devoted to studying the more complex and difficult three-dimensional problems with Neumann boundary condition and both the direct and indirect methods are adopted to construct combined boundary integral equations. The strongly-singular and hyper-singular integral operators are reformulated into compositions of weakly-singular integral operators and tangential-derivative operators, which allow us to prove the jump relations associated with the poroelastic layer potentials and boundary integral operators in a simple manner. Relying on both the investigated spectral properties of the strongly-singular operators, which indicate that the corresponding eigenvalues accumulate at three points whose values are only dependent on two Lam\'e constants, and the spectral properties of the Calder\'on relations of the poroelasticity, we propose low-GMRES-iteration regularized integral equations. Numerical examples are presented to demonstrate the accuracy and efficiency of the proposed methodology by means of a Chebyshev-based rectangular-polar solver.
In this work, we consider fracture propagation in nearly incompressible and (fully) incompressible materials using a phase-field formulation. We use a mixed form of the elasticity equation to overcome volume locking effects and develop a robust, nonlinear and linear solver scheme and preconditioner for the resulting system. The coupled variational inequality system, which is solved monolithically, consists of three unknowns: displacements, pressure, and phase-field. Nonlinearities due to coupling, constitutive laws, and crack irreversibility are solved using a combined Newton algorithm for the nonlinearities in the partial differential equation and employing a primal-dual active set strategy for the crack irreverrsibility constraint. The linear system in each Newton step is solved iteratively with a flexible generalized minimal residual method (GMRES). The key contribution of this work is the development of a problem-specific preconditioner that leverages the saddle-point structure of the displacement and pressure variable. Four numerical examples in pure solids and pressure-driven fractures are conducted on uniformly and locally refined meshes to investigate the robustness of the solver concerning the Poisson ratio as well as the discretization and regularization parameters.
A singularly perturbed parabolic problem of convection-diffusion type with a discontinuous initial condition is examined. An analytic function is identified which matches the discontinuity in the initial condition and also satisfies the homogenous parabolic differential equation associated with the problem. The difference between this analytical function and the solution of the parabolic problem is approximated numerically, using an upwind finite difference operator combined with an appropriate layer-adapted mesh. The numerical method is shown to be parameter-uniform. Numerical results are presented to illustrate the theoretical error bounds established in the paper.
A singularly perturbed parabolic problem of convection-diffusion type with a discontinuous initial condition is examined. A particular complimentary error function is identified which matches the discontinuity in the initial condition. The difference between this analytical function and the solution of the parabolic problem is approximated numerically. A co-ordinate transformation is used so that a layer-adapted mesh can be aligned to the interior layer present in the solution. Numerical analysis is presented for the associated numerical method, which establishes that the numerical method is a parameter-uniform numerical method. Numerical results are presented to illustrate the pointwise error bounds established in the paper.
Consider solving large sparse range symmetric singular linear systems $A{\bf x} = {\bf b}$ which arise, for instance, in the discretization of convection diffusion equations with periodic boundary conditions, and partial differential equations for electromagnetic fields using the edge-based finite element method. In theory, the Generalized Minimal Residual (GMRES) method converges to the least squares solution for inconsistent systems if the coefficient matrix $A$ is range symmetric, i.e. ${\rm R}(A)={ {\rm R}(A^{\rm T} })$, where ${\rm R}(A)$ is the range space of $A$. However, in practice, GMRES may not converge due to numerical instability. In order to improve the convergence, we propose using the pseudo-inverse for the solution of the severely ill-conditioned Hessenberg systems in GMRES. Numerical experiments on semi-definite inconsistent systems indicate that the method is efficient and robust. Finally, we further improve the convergence of the method, by reorthogonalizing the Modified Gram-Schmidt procedure.
Methods that align distributions by minimizing an adversarial distance between them have recently achieved impressive results. However, these approaches are difficult to optimize with gradient descent and they often do not converge well without careful hyperparameter tuning and proper initialization. We investigate whether turning the adversarial min-max problem into an optimization problem by replacing the maximization part with its dual improves the quality of the resulting alignment and explore its connections to Maximum Mean Discrepancy. Our empirical results suggest that using the dual formulation for the restricted family of linear discriminators results in a more stable convergence to a desirable solution when compared with the performance of a primal min-max GAN-like objective and an MMD objective under the same restrictions. We test our hypothesis on the problem of aligning two synthetic point clouds on a plane and on a real-image domain adaptation problem on digits. In both cases, the dual formulation yields an iterative procedure that gives more stable and monotonic improvement over time.