Efficient algorithms for computing rank-revealing factorizations on a GPU

Standard rank-revealing factorizations such as the singular value decomposition (SVD) and column pivoted QR factorization are challenging to implement efficiently on a GPU. A major difficulty in this regard is the inability of standard algorithms to cast most operations in terms of the Level-3 BLAS. This article presents two alternative algorithms for computing a rank-revealing factorization of the form $\mathsf{A}=\mathsf{U}\mathsf{T}{\mathsf{V}}^{\ast }$, where $\mathsf{U}$ and $\mathsf{V}$ are orthogonal and $\mathsf{T}$ is trapezoidal (or triangular if $\mathsf{A}$ is square). Both algorithms use randomized projection techniques to cast most of the flops in terms of matrix-matrix multiplication, which is exceptionally efficient on the GPU. Numerical experiments illustrate that these algorithms achieve significant acceleration over finely tuned GPU implementations of the SVD while providing low rank approximation errors close to that of the SVD.     

Backward error analysis of specified eigenpairs for sparse matrix polynomials

This article studies the unstructured and structured backward error analysis of specified eigenpairs for matrix polynomials. The structures we discuss include $T$-symmetric, $T$-skew-symmetric, Hermitian, skew Hermitian, $T$-even, $T$-odd, $H$-even, $H$-odd, $T$-palindromic, $T$-anti-palindromic, $H$-palindromic, and $H$-anti-palindromic matrix polynomials. Minimally structured perturbations are constructed with respect to Frobenius norm such that specified eigenpairs become exact eigenpairs of an appropriately perturbed matrix polynomial that also preserves sparsity. Further, we have used our results to solve various quadratic inverse eigenvalue problems that arise from real-life applications.     

Constructing the field of values of decomposable and general matrices using the ZNN based path following method

This paper describes and develops a fast and accurate path following algorithm that computes the field of values boundary curve $\partial F(A)$ for every conceivable complex or real square matrix $A$. It relies on the matrix flow decomposition algorithm that finds a proper block-diagonal flow representation for the associated hermitean matrix flow ${ℱ}_A(t)=\cos (t)H+\sin (t)K$ under unitary similarity if that is possible. Here ${ℱ}_A(t)$ is the 1-parameter-varying linear combination of the real and skew part matrices $H=(A+A^{\ast })/2$ and $K=(A-A^{\ast })/(2i)$ of $A$. For indecomposable matrix flows, ${ℱ}_A(t)$ has just one block and the ZNN based field of values algorithm works with ${ℱ}_A(t)$ directly. For decomposing flows ${ℱ}_A(t)$, the algorithm decomposes the given matrix $A$ unitarily into block-diagonal form $U^{\ast }AU=\text{diag}(A_j)$ with $j>1$ diagonal blocks $A_j$ whose individual sizes add up to the size of $A$. It then computes the field of values boundaries separately for each diagonal block $A_j$ using the path following ZNN eigenvalue method. The convex hull of all sub-fields of values boundary points $\partial F(A_j)$ finally determines the field of values boundary curve correctly for decomposing matrices $A$. The algorithm removes standard restrictions for path following FoV methods that generally cannot deal with decomposing matrices $A$ due to possible eigencurve crossings of ${ℱ}_A(t)$. Tests and numerical comparisons are included. Our ZNN based method is coded for sequential and parallel computations and both versions run very accurately and fast when compared with Johnson's Francis QR eigenvalue and Bendixson rectangle based method and compute global eigenanalyses of ${ℱ}_A(t_k)$ for large discrete sets of angles $t_k\in [0,2\pi ]$ more slowly.     

Fast solution of three-dimensional elliptic equations with randomly generated jumping coefficients by using tensor-structured preconditioners

In this paper, we propose and analyze the numerical algorithms for fast solution of periodic elliptic problems in random media in ${ℝ}^d$, $d=2,3$. Both the two-dimensional (2D) and three-dimensional (3D) elliptic problems are considered for the jumping equation coefficients built as a checkerboard type configuration of bumps randomly distributed on a large $L\times L$, or $L\times L\times L$ lattice, respectively. The finite element method discretization procedure on a 3D $n\times n\times n$ uniform tensor grid is described in detail, and the Kronecker tensor product approach is proposed for fast generation of the stiffness matrix. We introduce tensor techniques for the construction of the low Kronecker rank spectrally equivalent preconditioner in a periodic setting to be used in the framework of the preconditioned conjugate gradient iteration. The discrete 3D periodic Laplacian pseudo-inverse is first diagonalized in the Fourier basis, and then the diagonal matrix is reshaped into a fully populated third-order tensor of size $n\times n\times n$. The latter is approximated by a low-rank canonical tensor by using the multigrid Tucker-to-canonical tensor transform. As an example, we apply the presented solver in numerical analysis of stochastic homogenization method where the 3D elliptic equation should be solved many hundred times, and where for every random sampling of the equation coefficient one has to construct the new stiffness matrix and the right-hand side. The computational characteristics of the presented solver in terms of a lattice parameter $L$ and the grid-size, $n^d$, in both 2D and 3D cases are illustrated in numerical tests. Our solver can be used in various applications where the elliptic problem should be solved for a number of different coefficients for example, in many-particle dynamics, protein docking problems or stochastic modeling.     

Asymptotics for the eigenvalues of Toeplitz matrices with a symbol having a power singularity

The present work is devoted to the construction of an asymptotic expansion for the eigenvalues of a Toeplitz matrix $T_n(a)$ as $n$ goes to infinity, with a continuous and real-valued symbol $a$ having a power singularity of degree $\gamma$ with , at one point. The resulting matrix is dense and its entries decrease slowly to zero when moving away from the main diagonal, we apply the so called simple-loop (SL) method for constructing and justifying a uniform asymptotic expansion for all the eigenvalues. Note however, that the considered symbol does not fully satisfy the conditions imposed in previous works, but only in a small neighborhood of the singularity point. In the present work: (i) We construct and justify the asymptotic formulas of the SL method for the eigenvalues ${\lambda }_j(T_n(a))$ with $j⩾\epsilon n$, where the eigenvalues are arranged in nondecreasing order and $\epsilon$ is a sufficiently small fixed number. (ii) We show, with the help of numerical calculations, that the obtained formulas give good approximations in the case . (iii) We numerically show that the main term of the asymptotics for eigenvalues with , formally obtained from the formulas of the SL method, coincides with the main term of the asymptotics constructed and justified in the classical works of Widom and Parter.     

High relative accuracy with some special matrices related to Γ and β functions

For some families of totally positive matrices using $\mathrm{\Gamma }$ and $\beta$ functions, we provide their bidiagonal factorization. Moreover, when these functions are defined over integers, we prove that the bidiagonal factorization can be computed with high relative accuracy and so we can compute with high relative accuracy their eigenvalues, singular values, inverses and the solutions of some associated linear systems. We provide numerical examples illustrating this high relative accuracy.     

Tensor rank bounds and explicit QTT representations for the inverses of circulant matrices

In this paper, we are concerned with the inversion of circulant matrices and their quantized tensor-train (QTT) structure. In particular, we show that the inverse of a complex circulant matrix $A$, generated by the first column of the form ${(a_0,\dots ,a_{m-1},0,\dots ,0,a_{-n},\dots ,a_{-1})}^{\top }$ admits a QTT representation with the QTT ranks bounded by $(m+n)$. Under certain assumptions on the entries of $A$, we also derive an explicit QTT representation of $A^{-1}$. The latter can be used, for instance, to overcome stability issues arising when numerically solving differential equations with periodic boundary conditions in the QTT format.     

Lower and upper bounds of condition number for Vandermonde-wise matrices and method of fundamental solutions using pseudo radial-lines

Consider the method of fundamental solutions (MFS) for 2D Laplace's equation in a bounded simply connected domain $S$. In the standard MFS, the source nodes are located on a closed contour outside the domain boundary $\mathrm{\Gamma }(=\partial S)$, which is called pseudo-boundary. For circular, elliptic, and general closed pseudo-boundaries, analysis and computation have been studied extensively. New locations of source nodes are proposed along two pseudo radial-lines outside $\mathrm{\Gamma }$. Numerical results are very encouraging and promising. Since the success of the MFS mainly depends on stability, our efforts are focused on deriving the lower and upper bounds of condition number (Cond). The study finds stability properties of new Vandermonde-wise matrices on nodes $x_i\in [a,b]$ with . The Vandermonde-wise matrix is called in this article if it can be decomposed into the standard Vandermonde matrix. New lower and upper bounds of Cond are first derived for the standard Vandermonde matrix, and then for new algorithms of the MFS using two pseudo radial-lines. Both lower and upper bounds of Cond are intriguing in the stability study for the MFS. Numerical experiments are carried out to verify the stability analysis made. Since the fundamental solutions (as $\{\ln |\overline{{PQ}_i}|\}$) are the basis functions of the MFS, new Vandermonde-wise matrices are found. Since the nodes $x_i\in [a,b]$ with may come from approximations and interpolations by the Laurent polynomials with singular part, the conclusions in this article are important not only to the MFS but also to matrix analysis.     

Data-driven linear complexity low-rank approximation of general kernel matrices: A geometric approach

A general, rectangular kernel matrix may be defined as $K_{ij}=\kappa (x_i,y_j)$ where $\kappa (x,y)$ is a kernel function and where $X={\{x_i\}}_{i=1}^m$ and $Y={\{y_i\}}_{i=1}^n$ are two sets of points. In this paper, we seek a low-rank approximation to a kernel matrix where the sets of points $X$ and $Y$ are large and are arbitrarily distributed, such as away from each other, “intermingled”, identical, and so forth. Such rectangular kernel matrices may arise, for example, in Gaussian process regression where $X$ corresponds to the training data and $Y$ corresponds to the test data. In this case, the points are often high-dimensional. Since the point sets are large, we must exploit the fact that the matrix arises from a kernel function, and avoid forming the matrix, and thus ruling out most algebraic techniques. In particular, we seek methods that can scale linearly or nearly linearly with respect to the size of data for a fixed approximation rank. The main idea in this paper is to geometrically select appropriate subsets of points to construct a low rank approximation. An analysis in this paper guides how this selection should be performed.     

Recovering orthogonal tensors under arbitrarily strong,but locally correlated,noise

We consider the problem of recovering an orthogonally decomposable tensor with a subset of elements distorted by noise with arbitrarily large magnitude. We focus on the particular case where each mode in the decomposition is corrupted by noise vectors with components that are correlated locally, that is, with nearby components. We show that this deterministic tensor completion problem has the unusual property that it can be solved in polynomial time if the rank of the tensor is sufficiently large. This is the polar opposite of the low-rank assumptions of typical low-rank tensor and matrix completion settings. We show that our problem can be solved through a system of coupled Sylvester-like equations and show how to accelerate their solution by an alternating solver. This enables recovery even with a substantial number of missing entries, for instance for $n$-dimensional tensors of rank $n$ with up to $40\%$ missing entries.     

Nonlinear approximation of functions based on nonnegative least squares solver

In computational practice, most attention is paid to rational approximations of functions and approximations by the sum of exponents. We consider a wide enough class of nonlinear approximations characterized by a set of two required parameters. The approximating function is linear in the first parameter; these parameters are assumed to be positive. The individual terms of the approximating function represent a fixed function that depends nonlinearly on the second parameter. A numerical approximation minimizes the residual functional by approximating function values at individual points. The second parameter's value is set on a more extensive set of points of the interval of permissible values. The proposed approach's key feature consists in determining the first parameter on each separate iteration of the classical nonnegative least squares method. The computational algorithm is used to rational approximate the function . The second example concerns the approximation of the stretching exponential function at $x\ge 0$ by the sum of exponents.     

A parameterized extended shift-splitting preconditioner for nonsymmetric saddle point problems

In this article, a parameterized extended shift-splitting (PESS) method and its induced preconditioner are given for solving nonsingular and nonsymmetric saddle point problems with nonsymmetric positive definite (1,1) part. The convergence analysis of the $PESS$ iteration method is discussed. The distribution of eigenvalues of the preconditioned matrix is provided. A number of experiments are given to verify the efficiency of the $PESS$ method for solving nonsymmetric saddle-point problems.     

Multi-view side information-incorporated tensor completion

Tensor completion originates in numerous applications where data utilized are of high dimensions and gathered from multiple sources or views. Existing methods merely incorporate the structure information, ignoring the fact that ubiquitous side information may be beneficial to estimate the missing entries from a partially observed tensor. Inspired by this, we formulate a sparse and low-rank tensor completion model named SLRMV. The ${\ell }_0$-norm instead of its relaxation is used in the objective function to constrain the sparseness of noise. The CP decomposition is used to decompose the high-quality tensor, based on which the combination of Schatten $p$-norm on each latent factor matrix is employed to characterize the low-rank tensor structure with high computation efficiency. Diverse similarity matrices for the same factor matrix are regarded as multi-view side information for guiding the tensor completion task. Although SLRMV is a nonconvex and discontinuous problem, the optimality analysis in terms of Karush-Kuhn-Tucker (KKT) conditions is accordingly proposed, based on which a hard-thresholding based alternating direction method of multipliers (HT-ADMM) is designed. Extensive experiments remarkably demonstrate the efficiency of SLRMV in tensor completion.     

The rank of sparse random matrices

We determine the asymptotic normalized rank of a random matrix $\mathit{A}$ over an arbitrary field with prescribed numbers of nonzero entries in each row and column. As an application we obtain a formula for the rate of low-density parity check codes. This formula vindicates a conjecture of Lelarge (2013). The proofs are based on coupling arguments and a novel random perturbation, applicable to any matrix, that diminishes the number of short linear relations.     

Vanishing of higher order Alexander-type invariants of plane curves

The higher order degrees are Alexander-type invariants of complements to an affine plane curve. In this paper, we characterize the vanishing of such invariants for a curve C given as a transversal union of plane curves $C^{\prime }$ and $C^{\prime \prime }$ in terms of the finiteness and the vanishing properties of the invariants of $C^{\prime }$ and $C^{\prime \prime }$, and whether or not they are irreducible. As a consequence, we prove that the multivariable Alexander polynomial ${\mathrm{\Delta }}_C^{multi}$ is a power of $(t-1)$, and we characterize when ${\mathrm{\Delta }}_C^{multi}=1$ in terms of the defining equations of $C^{\prime }$ and $C^{\prime \prime }$. Our results impose obstructions on the class of groups that can be realized as fundamental groups of complements of a transversal union of curves.     

Elliptic Weingarten hypersurfaces of Riemannian products

Let $M^n$ be either a simply connected space form or a rank-one symmetric space of the noncompact type. We consider Weingarten hypersurfaces of $M\times \mathbb{R}$, which are those whose principal curvatures $k_1,\cdots ,k_n$ and angle function $𝛩$ satisfy a relation $W(k_1,\cdots ,k_n,{𝛩}^2)=0$, being W a differentiable function which is symmetric with respect to $k_1,\cdots ,k_n$. When $\partial W/\partial k_i>0$ on the positive cone of ${\mathbb{R}}^n$, a strictly convex Weingarten hypersurface determined by W is said to be elliptic. We show that, for a certain class of Weingarten functions W, there exist rotational strictly convex Weingarten hypersurfaces of $M\times \mathbb{R}$ which are either topological spheres or entire graphs over M. We establish a Jellett–Liebmann-type theorem by showing that a compact, connected and elliptic Weingarten hypersurface of either ${\mathbb{S}}^n\times \mathbb{R}$ or ${\mathbb{H}}^n\times \mathbb{R}$ is a rotational embedded sphere. Other uniqueness results for complete elliptic Weingarten hypersurfaces of these ambient spaces are obtained. We also obtain existence results for constant scalar curvature hypersurfaces of ${\mathbb{S}}^n\times \mathbb{R}$ and ${\mathbb{H}}^n\times \mathbb{R}$ which are either rotational or invariant by translations (parabolic or hyperbolic). We apply our methods to give new proofs of the main results by Manfio and Tojeiro on the classification of constant sectional curvature hypersurfaces of ${\mathbb{S}}^n\times \mathbb{R}$ and ${\mathbb{H}}^n\times \mathbb{R}$.     

On the volume functional of compact manifolds with harmonic Weyl tensor

The main aim of this article is to give the complete classification of critical metrics of the volume functional on a compact manifold M with boundary $\partial M$ and under the harmonic Weyl tensor condition. In particular, we prove that a critical metric with a harmonic Weyl tensor on a simply connected compact manifold with the boundary isometric to a standard sphere ${\mathbb{S}}^{n-1}$ must be isometric to a geodesic ball in a simply connected space form ${\mathbb{R}}^n$, ${\mathbb{H}}^n$, and ${\mathbb{S}}^n$. To this end, we first conclude the classification of such critical metrics under the Bach-flat assumption and then we prove that both geometric conditions are equivalent in this situation.