This paper considers regularized block multiconvex optimization, where the feasible set and objective function are generally nonconvex but convex in each block of variables. It also accepts nonconvex blocks and requires these blocks to be updated by proximal minimization. We review some interesting applications and propose a generalized block coordinate descent method. Under certain conditions, we show that any limit point satisfies the Nash equilibrium conditions. Furthermore, we establish global convergence and estimate the asymptotic convergence rate of the method by assuming a property based on the Kurdyka-Lojasiewicz inequality. The proposed algorithms are tested on nonnegative matrix and tensor factorization, as well as matrix and tensor recovery from incomplete observations. The tests include synthetic data and hyperspectral data, as well as image sets from the CBCL and ORL databases. Compared to the existing state-of-the-art algorithms, the proposed algorithms demonstrate superior performance in both speed and solution quality. The MATLAB code of nonnegative matrix/tensor decomposition and completion, along with a few demos, are accessible from the authors' homepages.

No abstract uploaded!

A recently developed nonlocal vector calculus is exploited to provide a variational analysis for a general class of nonlocal diffusion problems described by a linear integral equation on bounded domains in Rn. The nonlocal vector calculus also enables striking analogies to be drawn between the nonlocal model and classical models for diffusion, including a notion of nonlocal flux. The ubiquity of the nonlocal operator in applications is illustrated by a number of examples ranging from continuum mechanics to graph theory. In particular, it is shown that fractional Laplacian and fractional derivative models for anomalous diffusion are special cases of the nonlocal model for diffusion that we consider. The numerous applications elucidate different interpretations of the operator and the associated governing equations. For example, a probabilistic perspective explains that the nonlocal spatial operator appearing in our model corresponds to the infinitesimal generator for a symmetric jump process. Sufficient conditions on the kernel of the nonlocal operator and the notion of volume constraints are shown to lead to a well-posed problem. Volume constraints are a proxy for boundary conditions that may not be defined for a given kernel. In particular, we demonstrate for a general class of kernels that the nonlocal operator is a mapping between a volume constrained subspace of a fractional Sobolev subspace and its dual. We also demonstrate for other particular kernels that the inverse of the operator does not smooth but does correspond to diffusion. The impact of our results is that both a continuum analysis and a numerical method for the modeling of anomalous diffusion on bounded domains in Rn are provided. The analytical framework allows us to consider finite-dimensional approximations using discontinuous and continuous Galerkin methods, both of which are conforming for the nonlocal diffusion equation we consider; error and condition number estimates are derived.

Finite difference WENO schemes have established themselves as very worthy performers for entire classes of applications that involve hyperbolic conservation laws. In this paper we report on two major advances that make finite difference WENO schemes more efficient.

No abstract uploaded!

In this paper, we survey our recent work on designing high order positivity-preserving well-balanced finite difference and finite volume WENO (weighted essentially non-oscillatory) schemes, and discontinuous Galerkin finite element schemes for solving the shallow water equations with a non-flat bottom topography. These schemes are genuinely high order accurate in smooth regions for general solutions, are essentially non-oscillatory for general solutions with discontinuities, and at the same time they preserve exactly the water at rest or the more general moving water steady state solutions. A simple positivity-preserving limiter, valid under suitable CFL condition, has been introduced in one dimension and reformulated to two dimensions with triangular meshes, and we prove that the resulting schemes guarantee the positivity of the water depth.

For solving time-dependent convection-dominated partial differential equations (PDEs), which arise frequently in computational physics, high order numerical methods, including finite difference, finite volume, finite element and spectral methods, have been undergoing rapid developments over the past decades. In this article we give a brief survey of two selected classes of high order methods, namely the weighted essentially non-oscillatory (WENO) finite difference and finite volume schemes and discontinuous Galerkin (DG) finite element methods, emphasizing several of their recent developments: bound-preserving limiters for DG, finite volume and finite difference schemes, which address issues in robustness and accuracy; WENO limiters for DG methods, which address issues in non-oscillatory performance when there are strong shocks, and inverse Lax-Wendroff type boundary treatments for finite difference schemes, which address issues in solving complex geometry problems using Cartesian meshes.

Calculation of band structures of three dimensional photonic crystals amounts to solving large-scale Maxwell eigenvalue problems, which are notoriously challenging due to high multiplicity of zero eigenvalues. In this paper, we try to address this problem in such a broad context that band structures of three dimensional isotropic photonic crystals in all $14$ Bravais lattices can be efficiently computed in a unified framework. In this work, we uncover the delicate machinery behind several key results of our framework and on the basis of this new understanding we drastically simplify the derivations, proofs and arguments. Particular effort is made on reformulating the Bloch condition for all $14$ Bravais lattices in the redefined orthogonal coordinate system, and establishing eigen-decomposition of discrete partial derivative operators by identifying the hierarchical structure of the underlying normal (block) companion matrix, and reducing the eigen-decomposition of the double-curl operator to a simple factorization of a $3$-by-$3$ complex skew-symmetric matrix. With the validity of the novel nullspace free method in the broad context, we perform some calculations on one benchmark system to demonstrate the accuracy and efficiency of our algorithm to solve Maxwell eigenvalue problems.

Nonconvex optimization arises in many areas of computational science and engineering. However, most nonconvex optimization algorithms are only known to have local convergence or subsequence convergence properties. In this paper, we propose an algorithm for nonconvex optimization and establish its global convergence (of the whole sequence) to a critical point. In addition, we give its asymptotic convergence rate and numerically demonstrate its efficiency. In our algorithm, the variables of the underlying problem are either treated as one block or multiple disjoint blocks. It is assumed that each non-differentiable component of the objective function, or each constraint, applies only to one block of variables. The differentiable components of the objective function, however, can involve multiple blocks of variables together. Our algorithm updates one block of variables at a time by minimizing a certain prox-linear surrogate, along with an extrapolation to accelerate its convergence. The order of update can be either deterministically cyclic or randomly shuffled for each cycle. In fact, our convergence analysis only needs that each block be updated at least once in every fixed number of iterations. We show its global convergence (of the whole sequence) to a critical point under fairly loose conditions including, in particular, the Kurdyka–Łojasiewicz condition, which is satisfied by a broad class of nonconvex/nonsmooth applications. These results, of course, remain valid when the underlying problem is convex.We apply our convergence results to the coordinate descent iteration for non-convex regularized linear regression, as well as amodified rank-one residue iteration for nonnegative matrix factorization. We show that both applications have global convergence. Numerically,we tested our algorithm on nonnegative matrix and tensor factorization problems, where random shuffling clearly improves the chance to avoid low-quality local solutions.

AbstractWe introduce a sub-cell shock capturing method for scalar conservation laws built upon the Jin-Xin relaxation framework. Here, sub-cell shock capturing is achieved thanks to an original defect measure correction technique. The proposed correction exactly restores entropy shock solutions of the exact Riemann problem and moreover, it produces monotone and entropy satisfying approximate self-similar solutions. These solutions are then sampled using Glimm's random choice method to advance in time. The resulting scheme combines the simplicity of the Jin-Xin relaxation method with the resolution of the Glimm's scheme in the sharp capture of discontinuities. The strong benefit of using defect measure corrections over usual sub-cell shock capturing methods is that the scheme can be easily made entropy satisfying with respect to infinitely many entropy pairs. Consequently, under a classical CFL condition, the method is proved to converge to the unique entropy weak solution of the Cauchy problem for general non-linear flux functions.

No abstract uploaded!

In this paper, we establish sublinear and linear convergence of fixed point iterations generated by averaged operators in a Hilbert space. Our results are achieved under a bounded Holder regularity assumption which generalizes the well-known notion of bounded linear regularity. As an application of our results, we provide a convergence rate analysis for many important iterative methods in solving broad mathematical problems such as convex feasibility problems and variational inequality problems. These include Krasnoselskii–Mann iterations, the cyclic projection algorithm, forward-backward splitting and the Douglas–Rachford feasibility algorithm along with some variants. In the important case in which the underlying sets are convex sets described by convex polynomials in a finite dimensional space, we show that the Holder regularity properties are automatically satisfied, from which sublinear convergence follows.

In this paper, we propose and analyze a mixed $H^1$-conforming finite element method for solving Maxwell's equations in terms of electric field and Lagrange multiplier, where the multiplier is introduced accounting for the divergence constraint. We mainly focus on the case that the physical domain is non-convex and its boundary includes reentrant corners or edges, which may lead the solution of Maxwell's equations to be a non-$H^1$ very weak function and thus causes many numerical difficulties. The proposed method is formulated in the stabilized form by adding an additional mesh-dependent stabilization term to the mixed variational formulation. A general framework of stability and error analysis is established. Specifically, a pair of $H^1$-conforming finite elements, namely the $CP_2$-$P_1$ elements, for electric field and multiplier is studied and its stability and error bounds are also derived. Numerical experiments for source problems as well as eigenvalue problems on the $L$-shaped and cracked domains are presented to illustrate the high performance of the proposed mixed $H^1$-conforming finite element method.

In this work a simple method to enforce the positivity-preserving property for general high-order conservative schemes is proposed. The method detects critical numerical fluxes which may lead to negative density and pressure, and then imposes a simple flux limiter combining the high-order numerical flux with the first-order Lax-Friedrichs flux to satisfy a sufficient condition for preserving positivity. Though an extra time-step size condition is required to maintain the formal order of accuracy, it is less restrictive than those in previous works. A number of numerical examples suggest that this method, when applied on an essentially non-oscillatory scheme, can be used to prevent positivity failure when the flow involves vacuum or near vacuum and very strong discontinuities.

No abstract uploaded!

We analyze discontinuous Galerkin methods using upwind-biased numerical fluxes for time-dependent linear conservation laws. In one dimension, optimal a priori error estimates of order $k+1$ are obtained for the semidiscrete scheme when piecewise polynomials of degree at most $k$ $(k \ge 0)$ are used. Our analysis is valid for arbitrary nonuniform regular meshes and for both periodic boundary conditions and for initial-boundary value problems. We extend the analysis to the multidimensional case on Cartesian meshes when piecewise tensor product polynomials are used, and to the fully discrete scheme with explicit Runge--Kutta time discretization. Numerical experiments are shown to demonstrate the theoretical results.

No abstract uploaded!

In this paper, we design a new class of central compact schemes based on the cell-centered compact schemes of Lele {[}S.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992) 16-42{]}. These schemes equate a weighted sum of the nodal derivatives of a smooth function to a weighted sum of the function on both the grid points (cell boundaries) and the cell-centers. In our approach, instead of using a compact interpolation to compute the values on cell-centers, the physical values on these half grid points are stored as independent variables and updated using the same scheme as the physical values on the grid points. This approach increases the memory requirement but not the computational costs. Through systematic Fourier analysis and numerical tests, we observe that the schemes have excellent properties of high order, high resolution and low dissipation. It is an ideal class of schemes for the simulation of multi-scale problems such as aeroacoustics and turbulence.

No abstract uploaded!

Finding the stationary states of a free energy functional is an important problem in phase field crystal (PFC) models. Many efforts have been devoted to designing numerical schemes with energy dissipation and mass conservation properties. However, most existing approaches are time-consuming due to the requirement of small effective step sizes. In this paper, we discretize the energy functional and propose efficient numerical algorithms for solving the constrained nonconvex minimization problem. A class of gradient-based approaches, which are the so-called adaptive accelerated Bregman proximal gradient (AA-BPG) methods, is proposed, and the convergence property is established without the global Lipschitz constant requirements. A practical Newton method is also designed to further accelerate the local convergence with convergence guarantee. One key feature of our algorithms is that the energy dissipation and mass conservation properties hold during the iteration process. Moreover, we develop a hybrid acceleration framework to accelerate the AA-BPG methods and most of the existing approaches through coupling with the practical Newton method. Extensive numerical experiments, including two three-dimensional periodic crystals in the Landau--Brazovskii (LB) model and a two-dimensional quasicrystal in the Lifshitz--Petrich (LP) model, demonstrate that our approaches have adaptive step sizes which lead to a significant acceleration over many existing methods when computing complex structures.

Finding a fixed point to a nonexpansive operator, i.e., x = Tx, abstracts many problems in numerical linear algebra, optimization, and other areas of data sciences. To solve xed- point problems, we propose ARock, an algorithmic framework in which multiple agents (machines, processors, or cores) update x in an asynchronous parallel fashion. Asynchrony is crucial to parallel computing since it reduces synchronization wait, relaxes communication bottleneck, and thus speeds up computing significantly. At each step of ARock, an agent updates a randomly selected coordinate xi based on possibly out-of-date information on x. The agents share x through either global memory or communication. If writing xi is atomic, the agents can read and write x without memory locks. We prove that if the nonexpansive operator T has a fixed point, then with probability one, ARock generates a sequence that converges to a fixed point of T. Our conditions on T and step sizes are weaker than comparable work. Linear convergence is obtained under suitable assumptions. We propose special cases of ARock for linear systems, convex optimization, machine learning, as well as distributed and decentralized consensus problems. Numerical experiments of solving sparse logistic regression problems are presented.

Whether the three-dimensional incompressible Euler equations can develop a singularity in finite time from smooth initial data is one of the most challenging problems in mathematical fluid dynamics. This work attempts to provide an affirmative answer to this long-standing open question from a numerical point of view by presenting a class of potentially singular solutions to the Euler equations computed in axisymmetric geometries. The solutions satisfy a periodic boundary condition along the axial direction and a no-flow boundary condition on the solid wall. The equations are discretized in space using a hybrid 6th-order Galerkin and 6th-order finite difference method on specially designed adaptive (moving) meshes that are dynamically adjusted to the evolving solutions. With a maximum effective resolution of over $(3 × 10^{12})^2$ near the point of the singularity, we are able to advance the solution up to $\tau_2 = 0.003505$ and predict a singularity time of $t_s ≈ 0.0035056$, while achieving a pointwise relative error of $O(10^{−4})$ in the vorticity vector ω and observing a $(3 × 10^8)$-fold increase in the maximum vorticity $\|ω\|_{\infty}$. The numerical data are checked against all major blowup/non-blowup criteria, including Beale–Kato–Majda, Constantin–Fefferman–Majda, and Deng–Hou–Yu, to confirm the validity of the singularity. A local analysis near the point of the singularity also suggests the existence of a self-similar blowup in the meridian plane.

No abstract uploaded!

We develop an alternative formulation of conservative finite difference weighted essentially non-oscillatory (WENO) schemes to solve conservation laws. In this formulation, the WENO interpolation of the solution and its derivatives are used to directly construct the numerical flux, instead of the usual practice of reconstructing the flux functions. Even though this formulation is more expensive than the standard formulation, it does have several advantages. The first advantage is that arbitrary monotone fluxes can be used in this framework, while the traditional practice of reconstructing flux functions can be applied only to smooth flux splitting. The second advantage, which is fully explored in this paper, is that a narrower effective stencil is used compared with previous high order finite difference WENO schemes based on the reconstruction of flux functions, with a Lax-Wendroff time discretization. We will describe the scheme formulation and present numerical tests for one- and two-dimensional scalar and system conservation laws demonstrating the designed high order accuracy and non-oscillatory performance of the schemes constructed in this paper.

H(div) conforming and discontinuous Galerkin (DG) methods are designed for incompressible Euler's equation in two and three dimension. Error estimates are proved for both the semi-discrete method and fully-discrete method using backward Euler time stepping. Numerical examples exhibiting the performance of the methods are given.