On the perturbation of rank‐one symmetric tensors

The problem of symmetric rank‐one approximation of symmetric tensors is important in independent components analysis, also known as blind source separation, as well as polynomial optimization. We derive several perturbative results that are relevant to the well‐posedness of recovering rank‐one structure from approximately‐rank‐one symmetric tensors. We also specialize the analysis of the shifted symmetric higher‐order power method, an algorithm for computing symmetric tensor eigenvectors, to approximately‐rank‐one symmetric tensors. Copyright © 2012 John Wiley & Sons, Ltd.     

Shifted power method for computing tensor H‐eigenpairs

In this paper, we propose a shifted symmetric higher‐order power method for computing the H‐eigenpairs of a real symmetric even‐order tensor. The local convergence of the method is proved. In addition, by utilizing the fixed‐point analysis, we can characterize exactly which H‐eigenpairs can be found and which cannot be found by the method. Numerical examples are presented to illustrate the performance of the method. Copyright © 2016 John Wiley & Sons, Ltd.     

A semidefinite program approach for computing the maximum eigenvalue of a class of structured tensors and its applications in hypergraphs and copositivity test

Finding the maximum eigenvalue of a symmetric tensor is an important topic in tensor computation and numerical multilinear algebra. In this paper, we introduce a new class of structured tensors called W‐tensors, which not only extends the well‐studied nonnegative tensors by allowing negative entries but also covers several important tensors arising naturally from spectral hypergraph theory. We then show that finding the maximum H‐eigenvalue of an even‐order symmetric W‐tensor is equivalent to solving a structured semidefinite program and hence can be validated in polynomial time. This yields a highly efficient semidefinite program algorithm for computing the maximum H‐eigenvalue of W‐tensors and is based on a new structured sums‐of‐squares decomposition result for a nonnegative polynomial induced by W‐tensors. Numerical experiments illustrate that the proposed algorithm can successfully find the maximum H‐eigenvalue of W‐tensors with dimension up to 10,000, subject to machine precision. As applications, we provide a polynomial time algorithm for computing the maximum H‐eigenvalues of large‐size Laplacian tensors of hyperstars and hypertrees, where the algorithm can be up to 13 times faster than the state‐of‐the‐art numerical method introduced by Ng, Qi, and Zhou in 2009. Finally, we also show that the proposed algorithm can be used to test the copositivity of a multivariate form associated with symmetric extended Z‐tensors, whose order may be even or odd.     

Fast tensor method for summation of long‐range potentials on 3D lattices with defects

In this paper, we present a method for fast summation of long‐range potentials on 3D lattices with multiple defects and having non‐rectangular geometries, based on rank‐structured tensor representations. This is a significant generalization of our recent technique for the grid‐based summation of electrostatic potentials on the rectangular L × L × L lattices by using the canonical tensor decompositions and yielding the O(L) computational complexity instead of O(L³) by traditional approaches. The resulting lattice sum is calculated as a Tucker or canonical representation whose directional vectors are assembled by the 1D summation of the generating vectors for the shifted reference tensor, once precomputed on large N × N × N representation grid in a 3D bounding box. The tensor numerical treatment of defects is performed in an algebraic way by simple summation of tensors in the canonical or Tucker formats. To diminish the considerable increase in the tensor rank of the resulting potential sum, the ?‐rank reduction procedure is applied based on the generalized reduced higher‐order SVD scheme. For the reduced higher‐order SVD approximation to a sum of canonical/Tucker tensors, we prove the stable error bounds in the relative norm in terms of discarded singular values of the side matrices. The required storage scales linearly in the 1D grid‐size, O(N), while the numerical cost is estimated by O(NL). The approach applies to a general class of kernel functions including those for the Newton, Slater, Yukawa, Lennard‐Jones, and dipole‐dipole interactions. Numerical tests confirm the efficiency of the presented tensor summation method; we demonstrate that a sum of millions of Newton kernels on a 3D lattice with defects/impurities can be computed in seconds in Matlab implementation. The tensor approach is advantageous in further functional calculus with the lattice potential sums represented on a 3D grid, like integration or differentiation, using tensor arithmetics of 1D complexity. Copyright © 2015 John Wiley & Sons, Ltd.     

Application of Richardson extrapolation to the numerical solution of partial differential equations

Richardson extrapolation is a methodology for improving the order of accuracy of numerical solutions that involve the use of a discretization size h. By combining the results from numerical solutions using a sequence of related discretization sizes, the leading order error terms can be methodically removed, resulting in higher order accurate results. Richardson extrapolation is commonly used within the numerical approximation of partial differential equations to improve certain predictive quantities such as the drag or lift of an airfoil, once these quantities are calculated on a sequence of meshes, but it is not widely used to determine the numerical solution of partial differential equations. Within this article, Richardson extrapolation is applied directly to the solution algorithm used within existing numerical solvers of partial differential equations to increase the order of accuracy of the numerical result without referring to the details of the methodology or its implementation within the numerical code. Only the order of accuracy of the existing solver and certain interpolations required to pass information between the mesh levels are needed to improve the order of accuracy and the overall solution accuracy. Using the proposed methodology, Richardson extrapolation is used to increase the order of accuracy of numerical solutions of the linear heat and wave equations and of the nonlinear St. Venant equations in one‐dimension. © 2008 Wiley Periodicals, Inc. Numer Methods Partial Differential Eq, 2009     

A steady barycentric Lagrange interpolation method for the 2D higher‐order time‐fractional telegraph equation with nonlocal boundary condition with error analysis

In this paper, we consider the numerical solution of the time‐fractional telegraph equation with a nonlocal boundary condition. A novel barycentric Lagrange interpolation collocation method is developed to solve this equation. Two difficulties have been sorted: the singularity of the integration and the higher accuracy. At the same, we put forward a steady barycentric Lagrange interpolation technique to overcome the new “Runge” phenomenon in computation. Error estimates of the barycentric Lagrange interpolation and the time‐fractional telegraph system for the present method are presented in Sobolev spaces. High convergence rates of the proposed method are obtained and are consisted with the numerical values. Especially in the time dimension, we get the error bound, for h‐refinement and for n_t‐density in the L₂ norms. The numerical results obtained show that the proposed numerical algorithm is accurate and computationally efficient for solving time‐fractional telegraph equation. Experiments demonstrate the high convergence rates of the proposed method are consisted with the theoretical values.     

Two‐level Fourier analysis of multigrid for higher‐order finite‐element discretizations of the Laplacian

In this paper, we employ local Fourier analysis (LFA) to analyze the convergence properties of multigrid methods for higher‐order finite‐element approximations to the Laplacian problem. We find that the classical LFA smoothing factor, where the coarse‐grid correction is assumed to be an ideal operator that annihilates the low‐frequency error components and leaves the high‐frequency components unchanged, fails to accurately predict the observed multigrid performance and, consequently, cannot be a reliable analysis tool to give good performance estimates of the two‐grid convergence factor. While two‐grid LFA still offers a reliable prediction, it leads to more complex symbols that are cumbersome to use to optimize parameters of the relaxation scheme, as is often needed for complex problems. For the purposes of this analytical optimization as well as to have simple predictive analysis, we propose a modification that is “between” two‐grid LFA and smoothing analysis, which yields reasonable predictions to help choose correct damping parameters for relaxation. This exploration may help us better understand multigrid performance for higher‐order finite element discretizations, including for Q₂‐Q₁ (Taylor‐Hood) elements for the Stokes equations. Finally, we present two‐grid and multigrid experiments, where the corrected parameter choice is shown to yield significant improvements in the resulting two‐grid and multigrid convergence factors.     

Further results on the Drazin inverse of even‐order tensors

The notion of the Drazin inverse of an even‐order tensors with the Einstein product was introduced, very recently [J. Ji and Y. Wei. Comput. Math. Appl., 75(9), (2018), pp. 3402‐3413]. In this article, we further elaborate this theory by establishing a few characterizations of the Drazin inverse and $\mathrm{??}$‐weighted Drazin inverse of tensors. In addition to these, we compute the Drazin inverse of tensors using different types of generalized inverses and full rank decomposition of tensors. We also address the solution of multilinear systems by using the Drazin inverse and iterative (higher order Gauss‐Seidel) method of tensors. Besides these, the convergence analysis of the iterative technique is also investigated within the framework of the Einstein product.     

A high order finite volume method for one dimensional nonlocal reactive flows of parabolic type

In this paper, we consider a high order finite volume approximation of one‐dimensional nonlocal reactive flows of parabolic type. The method is obtained by discretizing in space by arbitrary order vertex‐centered finite volumes, followed by a modified Simpson quadrature scheme for the time stepping. Compared to the existed finite volume methods, this new finite volume scheme could achieve the desired accuracy with less data storage by employing higher‐order trial spaces. The finite volume approximations are proved to possess optimal order convergence rates in the H¹‐norm and L²‐norm, which are also confirmed by numerical tests. Copyright © 2015 John Wiley & Sons, Ltd.     

An approach to the Gummel map by vector extrapolation methods

The numerical approximation of nonlinear partial differential equations requires the computation of large nonlinear systems, that are typically solved by iterative schemes. At each step of the iterative process, a large and sparse linear system has to be solved, and the amount of time elapsed per step grows with the dimensions of the problem. As a consequence, the convergence rate may become very slow, requiring massive cpu-time to compute the solution. In all such cases, it is important to improve the rate of convergence of the iterative scheme. This can be achieved, for instance, by vector extrapolation methods. In this work, we apply some vector extrapolation methods to the electronic device simulation to improve the rate of convergence of the family of Gummel decoupling algorithms. Furthermore, a different approach to the topological ε-algorithm is proposed and preliminary results are presented.     

The modified natural power method for principal component computation

A modified version of the natural power method (NP) for fast estimation and tracking of the principal eigenvectors of a vector sequence is Presented. It is an extension of the natural power method because it is a solution to obtain the principal eigenvectors and not only for tracking of the principal subspace. As compared with some power-based methods such as Oja method, the projection approximation subspace tracking (PAST) method, and the novel information criterion (NIC) method, the modified natural power method (MNP) has the fastest convergence rate and can be easily implemented with only O(np) flops of computation at each iteration, where n is the dimension of the vector sequence and p is the dimension of the principal subspace or the number of the principal eigenvectors. Furthermore, it is guaranteed to be globally and exponentially convergent in contrast with some non-power-based methods such as MALASE and OPERA. Selected from Journal of Fudan University (Natural Science), 2004, 43(3): 275–284     

On the successive supersymmetric rank‐1 decomposition of higher‐order supersymmetric tensors

In this paper, a successive supersymmetric rank‐1 decomposition of a real higher‐order supersymmetric tensor is considered. To obtain such a decomposition, we design a greedy method based on iteratively computing the best supersymmetric rank‐1 approximation of the residual tensors. We further show that a supersymmetric canonical decomposition could be obtained when the method is applied to an orthogonally diagonalizable supersymmetric tensor, and in particular, when the order is 2, this method generates the eigenvalue decomposition for symmetric matrices. Details of the algorithm designed and the numerical results are reported in this paper. Copyright © 2007 John Wiley & Sons, Ltd.     

Convergence analysis for second‐order accurate schemes for the periodic nonlocal Allen‐Cahn and Cahn‐Hilliard equations

In this paper, we provide a detailed convergence analysis for fully discrete second‐order (in both time and space) numerical schemes for nonlocal Allen‐Cahn and nonlocal Cahn‐Hilliard equations. The unconditional unique solvability and energy stability ensures ? ⁴ stability. The convergence analysis for the nonlocal Allen‐Cahn equation follows the standard procedure of consistency and stability estimate for the numerical error function. For the nonlocal Cahn‐Hilliard equation, because of the complicated form of the nonlinear term, a careful expansion of its discrete gradient is undertaken, and an H ^?1 inner‐product estimate of this nonlinear numerical error is derived to establish convergence. In addition, an a priori bound of the numerical solution at the discrete level is needed in the error estimate. Such a bound can be obtained by performing a higher order consistency analysis by using asymptotic expansions for the numerical solution. Following the technique originally proposed by Strang (eg, 1964), instead of the standard comparison between the exact and numerical solutions, an error estimate between the numerical solution and the constructed approximate solution yields an O (s ³+h ⁴) convergence in norm, in which s and h denote the time step and spatial mesh sizes, respectively. This in turn leads to the necessary bound under a standard constraint s ≤C h . Here, we also prove convergence of the scheme in the maximum norm under the same constraint.     

A higher order difference method for singularly perturbed parabolic partial differential equations

This paper studies a higher order numerical method for the singularly perturbed parabolic convection-diffusion problems where the diffusion term is multiplied by a small perturbation parameter. In general, the solutions of these type of problems have a boundary layer. Here, we generate a spatial adaptive mesh based on the equidistribution of a positive monitor function. Implicit Euler method is used to discretize the time variable and an upwind scheme is considered in space direction. A higher order convergent solution with respect to space and time is obtained using the postprocessing based extrapolation approach. It is observed that the convergence is independent of perturbation parameter. This technique enhances the order of accuracy from first order uniform convergence to second order uniform convergence in space as well as in time. Comparative study with the existed meshes show the highly effective behavior of the present method.     

Low‐rank approximation of tensors via sparse optimization

The goal of this paper is to find a low‐rank approximation for a given nth tensor. Specifically, we give a computable strategy on calculating the rank of a given tensor, based on approximating the solution to an NP‐hard problem. In this paper, we formulate a sparse optimization problem via an l₁‐regularization to find a low‐rank approximation of tensors. To solve this sparse optimization problem, we propose a rescaling algorithm of the proximal alternating minimization and study the theoretical convergence of this algorithm. Furthermore, we discuss the probabilistic consistency of the sparsity result and suggest a way to choose the regularization parameter for practical computation. In the simulation experiments, the performance of our algorithm supports that our method provides an efficient estimate on the number of rank‐one tensor components in a given tensor. Moreover, this algorithm is also applied to surveillance videos for low‐rank approximation.     

Spectral Algorithms for Tensor Completion

In the tensor completion problem, one seeks to estimate a low‐rank tensor based on a random sample of revealed entries. In terms of the required sample size, earlier work revealed a large gap between estimation with unbounded computational resources (using, for instance, tensor nuclear norm minimization) and polynomial‐time algorithms. Among the latter, the best statistical guarantees have been proved, for third‐order tensors, using the sixth level of the sum‐of‐squares (sos ) semidefinite programming hierarchy. However, the sos approach does not scale well to large problem instances. By contrast, spectral methods—based on unfolding or matricizing the tensor—are attractive for their low complexity, but have been believed to require a much larger sample size. This paper presents two main contributions. First, we propose a new method, based on unfolding, which outperforms naive ones for symmetric k^th‐order tensors of rank r. For this result we make a study of singular space estimation for partially revealed matrices of large aspect ratio, which may be of independent interest. For third‐order tensors, our algorithm matches the sos method in terms of sample size (requiring about rd^3/2 revealed entries), subject to a worse rank condition (r ≪ d^3/4 rather than r ≪ d^3/2). We complement this result with a different spectral algorithm for third‐order tensors in the overcomplete (r ≥ d) regime. Under a random model, this second approach succeeds in estimating tensors of rank d ≤ r ≪ d^3/2 from about rd^3/2 revealed entries. © 2018 Wiley Periodicals, Inc.     

Stability and convergence of two‐grid Crank‐Nicolson extrapolation scheme for the time‐dependent natural convection equations

In this paper, we consider the Crank‐Nicolson extrapolation scheme for the 2D/3D unsteady natural convection problem. Our numerical scheme includes the implicit Crank‐Nicolson scheme for linear terms and the recursive linear method for nonlinear terms. Standard Galerkin finite element method is used to approximate the spatial discretization. Stability and optimal error estimates are provided for the numerical solutions. Furthermore, a fully discrete two‐grid Crank‐Nicolson extrapolation scheme is developed, the corresponding stability and convergence results are derived for the approximate solutions. Comparison from aspects of the theoretical results and computational efficiency, the two‐grid Crank‐Nicolson extrapolation scheme has the same order as the one grid method for velocity and temperature in H¹‐norm and for pressure in L²‐norm. However, the two‐grid scheme involves much less work than one grid method. Finally, some numerical examples are provided to verify the established theoretical results and illustrate the performances of the developed numerical schemes.     

Iterative minimization of the Rayleigh quotient by block steepest descent iterations

The topic of this paper is the convergence analysis of subspace gradient iterations for the simultaneous computation of a few of the smallest eigenvalues plus eigenvectors of a symmetric and positive definite matrix pair (A,M). The methods are based on subspace iterations for A ^{? 1}M and use the Rayleigh‐Ritz procedure for convergence acceleration. New sharp convergence estimates are proved by generalizing estimates, which have been presented for vectorial steepest descent iterations (see SIAM J. Matrix Anal. Appl., 32(2):443‐456, 2011). Copyright © 2013 John Wiley & Sons, Ltd.     

Convergence of a transition probability tensor of a higher‐order Markov chain to the stationary probability vector

In this paper, first we introduce a new tensor product for a transition probability tensor originating from a higher‐order Markov chain. Subsequently, some properties of the new tensor product are explained, and its relationship with the stationary probability vector is studied. Also, similarity between results obtained by this new product and the first‐order case is shown. Furthermore, we prove the convergence of a transition probability tensor to the stationary probability vector. Finally, we show how to achieve a stationary probability vector with some numerical examples and make some comparison between the proposed method and another existing method for obtaining stationary probability vectors. Copyright © 2016 John Wiley & Sons, Ltd.     

A unitary joint diagonalization algorithm for nonsymmetric higher‐order tensors based on Givens‐like rotations

Based on Givens‐like rotations, we present a unitary joint diagonalization algorithm for a set of nonsymmetric higher‐order tensors. Each unitary rotation matrix only depends on one unknown parameter which can be analytically obtained in an independent way following a reasonable assumption and a complex derivative technique. It can serve for the canonical polyadic decomposition of a higher‐order tensor with orthogonal factors. Furthermore, based on cross‐high‐order cumulants of observed signals, we show that the proposed algorithm can be applied to solve the joint blind source separation problem. The simulation results reveal that the proposed algorithm has a competitive performance compared with those of several existing related methods.