Showing papers in "SIAM Journal on Matrix Analysis and Applications in 2014"

Nearly Linear Time Algorithms for Preconditioning and Solving Symmetric, Diagonally Dominant Linear Systems

Abstract: We present a randomized algorithm that on input a symmetric, weakly diagonally dominant $n$-by-$n$ matrix $A$ with $m$ nonzero entries and an $n$-vector $b$ produces an ${\tilde{x}} $ such that $\|{\tilde{x}} - A^{\dagger} {b} \|_{A} \leq \epsilon \|A^{\dagger} {b}\|_{A}$ in expected time $O (m \log^{c}n \log (1/\epsilon))$ for some constant $c$. By applying this algorithm inside the inverse power method, we compute approximate Fiedler vectors in a similar amount of time. The algorithm applies subgraph preconditioners in a recursive fashion. These preconditioners improve upon the subgraph preconditioners first introduced by Vaidya in 1990. For any symmetric, weakly diagonally dominant matrix $A$ with nonpositive off-diagonal entries and $k \geq 1$, we construct in time $O (m \log^{c} n)$ a preconditioner $B$ of $A$ with at most $2 (n - 1) + O ((m/k) \log^{39} n)$ nonzero off-diagonal entries such that the finite generalized condition number $\kappa_{f} (A,B)$ is at most $k$, for some other constant $c$. I...

Robust Low-Rank Tensor Recovery: Models and Algorithms

Abstract: Robust tensor recovery plays an instrumental role in robustifying tensor decompositions for multilinear data analysis against outliers, gross corruptions, and missing values and has a diverse array of applications. In this paper, we study the problem of robust low-rank tensor recovery in a convex optimization framework, drawing upon recent advances in robust principal component analysis and tensor completion. We propose tailored optimization algorithms with global convergence guarantees for solving both the constrained and the Lagrangian formulations of the problem. These algorithms are based on the highly efficient alternating direction augmented Lagrangian and accelerated proximal gradient methods. We also propose a nonconvex model that can often improve the recovery results from the convex models. We investigate the empirical recoverability properties of the convex and nonconvex formulations and compare the computational performance of the algorithms on simulated data. We demonstrate through a number o...

$M$-Tensors and Some Applications

Abstract: We introduce M -tensors. This concept extends the concept of M -matrices. We denote Z-tensors as the tensors with nonpositive off-diagonal entries. We show that M -tensors must be Z- tensors and the maximal diagonal entry must be nonnegative. The diagonal elements of a symmetric M -tensor must be nonnegative. A symmetric M -tensor is copositive. Based on the spectral theory of nonnegative tensors, we show that the minimal value of the real parts of all eigenvalues of an M - tensor is its smallest H + -eigenvalue and also is its smallest H-eigenvalue. We show that a Z-tensor is an M -tensor if and only if all its H + -eigenvalues are nonnegative. Some further spectral properties of M -tensors are given. We also introduce strong M -tensors, and some corresponding conclusions are given. In particular, we show that all H-eigenvalues of strong M -tensors are positive. We apply this property to study the positive definiteness of a class of multivariate forms associated with Z-tensors. We also propose an algorithm for testing the positive definiteness of such a multivariate form.

FEAST As A Subspace Iteration Eigensolver Accelerated By Approximate Spectral Projection

Abstract: The calculation of a segment of eigenvalues and their corresponding eigenvectors of a Hermitian matrix or matrix pencil has many applications. A new density-matrix-based algorithm has been proposed recently and a software package FEAST has been developed. The density-matrix approach allows FEAST's implementation to exploit a key strength of modern computer architectures, namely, multiple levels of parallelism. Consequently, the software package has been well received, especially in the electronic structure community. Nevertheless, theoretical analysis of FEAST has lagged. For instance, the FEAST algorithm has not been proven to converge. This paper offers a detailed numerical analysis of FEAST. In particular, we show that the FEAST algorithm can be understood as an accelerated subspace iteration algorithm in conjunction with the Rayleigh--Ritz procedure. The novelty of FEAST lies in its accelerator, which is a rational matrix function that approximates the spectral projector onto the eigenspace in questio...

All Real Eigenvalues of Symmetric Tensors

Abstract: This paper studies how to compute all real eigenvalues, associated to real eigenvectors, of a symmetric tensor. As is well known, the largest or smallest eigenvalue can be found by solving a polynomial optimization problem, while the other middle ones cannot. We propose a new approach for computing all real eigenvalues sequentially, from the largest to the smallest. It uses Jacobian semidefinite relaxations in polynomial optimization. We show that each eigenvalue can be computed by solving a finite hierarchy of semidefinite relaxations. Numerical experiments are presented to show how to do this.

Canonical Polyadic Decomposition of Third-Order Tensors: Reduction to Generalized Eigenvalue Decomposition

Abstract: Canonical polyadic decomposition (CPD) of a third-order tensor is decomposition in a minimal number of rank-1 tensors. We call an algorithm algebraic if it is guaranteed to find the decomposition when it is exact and if it relies only on standard linear algebra (essentially sets of linear equations and matrix factorizations). The known algebraic algorithms for the computation of the CPD are limited to cases where at least one of the factor matrices has full column rank. In this paper we present an algebraic algorithm for the computation of the CPD in cases where none of the factor matrices has full column rank. In particular, we show that if the famous Kruskal condition holds, then the CPD can be found algebraically.

Semidefinite Relaxations for Best Rank-1 Tensor Approximations

Abstract: This paper studies the problem of finding best rank-1 approximations for both symmetric and nonsymmetric tensors For symmetric tensors, this is equivalent to optimizing homogeneous polynomials over unit spheres; for nonsymmetric tensors, this is equivalent to optimizing multiquadratic forms over multispheres We propose semidefinite relaxations, based on sum of squares representations, to solve these polynomial optimization problems Their special properties and structures are studied In applications, the resulting semidefinite programs are often large scale The recent Newton-CG augmented Lagrangian method by Zhao, Sun, and Toh [SIAM J Optim, 20 (2010), pp 1737--1765] is suitable for solving these semidefinite relaxations Extensive numerical experiments are presented to show that this approach is efficient in getting best rank-1 approximations

Efficient and stable Arnoldi restarts for matrix functions based on quadrature

Abstract: When using the Arnoldi method for approximating $f(A){\mathbf b}$, the action of a matrix function on a vector, the maximum number of iterations that can be performed is often limited by the storage requirements of the full Arnoldi basis. As a remedy, different restarting algorithms have been proposed in the literature, none of which has been universally applicable, efficient, and stable at the same time. We utilize an integral representation for the error of the iterates in the Arnoldi method which then allows us to develop an efficient quadrature-based restarting algorithm suitable for a large class of functions, including the so-called Stieltjes functions and the exponential function. Our method is applicable for functions of Hermitian and non-Hermitian matrices, requires no a priori spectral information, and runs with essentially constant computational work per restart cycle. We comment on the relation of this new restarting approach to other existing algorithms and illustrate its efficiency and numer...

An algorithm for generic and low-rank specific identifiability of complex tensors ∗

Abstract: We propose a new sufficient condition for verifying whether general rank-$r$ complex tensors of arbitrary order admit a unique decomposition as a linear combination of rank-$1$ tensors. A practical algorithm is proposed for verifying this condition, with which it was established that in all spaces of dimension less than 15000, with a few known exceptions, listed in the paper, generic identifiability holds for ranks up to one less than the generic rank of the space. This is the largest possible rank value for which generic identifiability can hold, except for spaces with a perfect shape. The algorithm can also verify the identifiability of a given specific rank-$r$ decomposition, provided that it can be shown to correspond to a nonsingular point of the $r$th order secant variety. For sufficiently small rank, which nevertheless improves upon the known bounds for specific identifiability, some local equations of this variety are known, allowing us to verify this property. As a particular example of our appro...

Superfast and Stable Structured Solvers for Toeplitz Least Squares via Randomized Sampling

Abstract: We present some superfast ($O((m+n)\log^{2}(m+n))$ complexity) and stable structured direct solvers for $m\times n$ Toeplitz least squares problems. Based on the displacement equation, a Toeplitz matrix $T$ is first transformed into a Cauchy-like matrix $\mathcal{C}$, which can be shown to have small off-diagonal numerical ranks when the diagonal blocks are rectangular. We generalize standard hierarchically semiseparable (HSS) matrix representations to rectangular ones, and construct a rectangular HSS approximation to $\mathcal{C}$ in nearly linear complexity with randomized sampling and fast multiplications of $\mathcal{C}$ with vectors. A new URV HSS factorization and a URV HSS solution are designed for the least squares solution. We also present two structured normal equation methods. Systematic error and stability analysis for our HSS methods is given, which is also useful for studying other HSS and rank structured methods. We derive the growth factors and the backward error bounds in the HSS factoriz...

Nonnegative Tensor Factorization, Completely Positive Tensors, and a Hierarchical Elimination Algorithm

Abstract: Nonnegative tensor factorization has applications in statistics, computer vision, exploratory multiway data analysis, and blind source separation. A symmetric nonnegative tensor, which has an exact symmetric nonnegative factorization, is called a completely positive tensor. This concept extends the concept of completely positive matrices. A classical result in the theory of completely positive matrices is that a symmetric, diagonally dominated nonnegative matrix is a completely positive matrix. In this paper, we introduce strongly symmetric tensors and show that a symmetric tensor has a symmetric binary decomposition if and only if it is strongly symmetric. Then we show that a strongly symmetric, hierarchically dominated nonnegative tensor is a completely positive tensor, and present a hierarchical elimination algorithm for checking this. Numerical examples are given to illustrate this. Some other properties of completely positive tensors are discussed. In particular, we show that the completely positive ...

A New Analysis of Block Preconditioners for Saddle Point Problems

Abstract: We consider symmetric saddle point matrices. We analyze block preconditioners based on the knowledge of a good approximation for both the top left block and the Schur complement resulting from its elimination. We obtain bounds on the eigenvalues of the preconditioned matrix that depend only of the quality of these approximations, as measured by the related condition numbers. Our analysis applies to indefinite block diagonal preconditioners, block triangular preconditioners, inexact Uzawa preconditioners, block approximate factorization preconditioners, and a further enhancement of these preconditioners based on symmetric block Gauss--Seidel-type iterations. The analysis is unified and allows the comparison of these different approaches. In particular, it reveals that block triangular and inexact Uzawa preconditioners lead to identical eigenvalue distributions. These theoretical results are illustrated on the discrete Stokes problem. It turns out that the provided bounds allow one to localize accurately bo...

An Adaptive Shifted Power Method for Computing Generalized Tensor Eigenpairs

Abstract: Several tensor eigenpair definitions have been put forth in the past decade, but these can all be unified under generalized tensor eigenpair framework, introduced by Chang, Pearson, and Zhang [J. Math. Anal. Appl., 350 (2009), pp. 416--422]. Given mth-order, n-dimensional real-valued symmetric tensors ${\mathscr{A}}$ and $\boldsymbol{\mathscr{B}}$, the goal is to find $\lambda \in \mathbb{R}$ and $\mathbf{x} \in \mathbb{R}^{n}, \mathbf{x} eq 0$ such that ${\mathscr{A}}\mathbf{x}^{m-1} = \lambda {\mathscr{B}}\mathbf{x}^{m-1}$. Different choices for ${\mathscr{B}}$ yield different versions of the tensor eigenvalue problem. We present our generalized eigenproblem adaptive power (GEAP) method for solving the problem, which is an extension of the shifted symmetric higher-order power method (SS-HOPM) for finding Z-eigenpairs. A major drawback of SS-HOPM is that its performance depended on choosing an appropriate shift, but our GEAP method also includes an adaptive method for choosing the shift automatically.

A Symmetry Preserving Algorithm for Matrix Scaling

Abstract: We present an iterative algorithm which asymptotically scales the $\infty$-norm of each row and each column of a matrix to one. This scaling algorithm preserves symmetry of the original matrix and shows fast linear convergence with an asymptotic rate of 1/2. We discuss extensions of the algorithm to the 1-norm, and by inference to other norms. For the 1-norm case, we show again that convergence is linear, with the rate dependent on the spectrum of the scaled matrix. We demonstrate experimentally that the scaling algorithm improves the conditioning of the matrix and that it helps direct solvers by reducing the need for pivoting. In particular, for symmetric matrices the theoretical and experimental results highlight the potential of the proposed algorithm over existing alternatives.

Numerical Optimization of Eigenvalues of Hermitian Matrix Functions

Abstract: This work concerns the global minimization of a prescribed eigenvalue or a weighted sum of prescribed eigenvalues of a Hermitian matrix-valued function depending on its parameters analytically in a box. We describe how the analytical properties of eigenvalue functions can be put into use to derive piecewise quadratic functions that underestimate the eigenvalue functions. These piecewise quadratic underestimators lead us to a global minimization algorithm, originally due to Breiman and Cutler. We prove the global convergence of the algorithm and show that it can be effectively used for the minimization of extreme eigenvalues, e.g., the largest eigenvalue or the sum of the largest specified number of eigenvalues. This is particularly facilitated by the analytical formulas for the first derivatives of eigenvalues, as well as analytical lower bounds on the second derivatives that can be deduced for extreme eigenvalue functions. The applications that we have in mind also include the ${\rm H}_\infty$-norm of a ...

Exact solutions in structured low-rank approximation ∗

Abstract: Structured low-rank approximation is the problem of minimizing a weighted Frobenius distance to a given matrix among all matrices of fixed rank in a linear space of matrices. We study the critical points of this optimization problem using algebraic geometry. A particular focus lies on Hankel matrices, Sylvester matrices, and generic linear spaces.

On the Convergence of the Self-Consistent Field Iteration in Kohn--Sham Density Functional Theory

Abstract: It is well known that the self-consistent field (SCF) iteration for solving the Kohn--Sham (KS) equation often fails to converge, yet there is no clear explanation. In this paper, we investigate the SCF iteration from the perspective of minimizing the corresponding KS total energy functional. By analyzing the second-order Taylor expansion of the KS total energy functional and estimating the relationship between the Hamiltonian and the part of the Hessian which is not used in the SCF iteration, we are able to prove global convergence from an arbitrary initial point and local linear convergence from an initial point sufficiently close to the solution of the KS equation under the assumptions that the gap between the occupied states and unoccupied states is sufficiently large and the second-order derivatives of the exchange correlation functional are uniformly bounded from above. Although these conditions are very stringent and are almost never satisfied in reality, our analysis is interesting in the sense th...

A Residual Replacement Strategy for Improving the Maximum Attainable Accuracy of $s$-Step Krylov Subspace Methods

Abstract: Krylov subspace methods are a popular class of iterative methods for solving linear systems with large, sparse matrices. On modern computer architectures, both sequential and parallel performance of classical Krylov methods is limited by costly data movement, or communication, required to update the approximate solution in each iteration. These motivated communication-avoiding Krylov methods, based on $s$-step formulations, reduce data movement by a factor of $O(s)$ by reordering the computations in classical Krylov methods to exploit locality. Studies on the finite precision behavior of communication-avoiding Krylov methods in the literature have thus far been empirical in nature; in this work, we provide the first quantitative analysis of the maximum attainable accuracy of communication-avoiding Krylov subspace methods in finite precision. Following the analysis for classical Krylov methods, we derive a bound on the deviation of the true and updated residuals in communication-avoiding conjugate gradient...

Convergence of restarted Krylov subspace methods for Stieltjes functions of matrices

Abstract: To approximate $f(A){b}$---the action of a matrix function on a vector---by a Krylov subspace method, restarts may become mandatory due to storage requirements for the Arnoldi basis or due to the growing computational complexity of evaluating $f$ on a Hessenberg matrix of growing size. A number of restarting methods have been proposed in the literature in recent years and there has been substantial algorithmic advancement concerning their stability and computational efficiency. However, the question under which circumstances convergence of these methods can be guaranteed has remained largely unanswered. In this paper we consider the class of Stieltjes functions and a related class, which contain important functions like the (inverse) square root and the matrix logarithm. For these classes of functions we present new theoretical results which prove convergence for Hermitian positive definite matrices $A$ and arbitrary restart lengths. We also propose a modification of the Arnoldi approximation which guaran...

On the Global Convergence of the Alternating Least Squares Method for Rank-One Approximation to Generic Tensors

Abstract: Tensor decomposition has important applications in various disciplines, but it remains an extremely challenging task even to this date. A slightly more manageable endeavor has been to find a low rank approximation in place of the decomposition. Even for this less stringent undertaking, it is an established fact that tensors beyond matrices can fail to have best low rank approximations, with the notable exception that the best rank-one approximation always exists for tensors of any order. Toward the latter, the most popular approach is the notion of alternating least squares whose specific numerical scheme appears in the form as a variant of the power method. Though the limiting behavior of the objective values is well understood, a proof of global convergence for the iterates themselves has been elusive. This paper partially addresses the missing piece by showing that for almost all tensors, the iterates generated by the alternating least squares method for the rank-one approximation converge globally. Th...

Adaptive Tangential Interpolation in Rational Krylov Subspaces for MIMO Dynamical Systems

Abstract: Model reduction approaches have been shown to be powerful techniques in the numerical simulation of very large dynamical systems. The presence of multiple inputs and outputs (MIMO systems) makes the reduction process even more challenging. We consider projection-based approaches where the reduction of complexity is achieved by direct projection of the problem onto a rational Krylov subspace of significantly smaller dimension. We present an effective way to treat multiple inputs by dynamically choosing the next direction vectors to expand the space. We apply the new strategy to the approximation of the transfer matrix function and to the solution of the Lyapunov matrix equation. Numerical results confirm that the new approach is competitive with respect to state-of-the-art methods both in terms of CPU time and memory requirements.

Subspace Methods For Computing The Pseudospectral Abscissa And The Stability Radius

Abstract: The pseudospectral abscissa and the stability radius are well-established tools for quantifying the stability of a matrix under unstructured perturbations. Based on first-order eigenvalue expansions, Guglielmi and Overton [SIAM J. Matrix Anal. Appl., 32 (2011), pp. 1166--1192] recently proposed a linearly converging iterative method for computing the pseudospectral abscissa. In this paper, we propose to combine this method and its variants with subspace acceleration. Each extraction step computes the pseudospectral abscissa of a small rectangular matrix pencil, which is comparably cheap and guarantees monotonicity. We observe local quadratic convergence and prove local superlinear convergence of the resulting subspace methods. Moreover, these methods extend naturally to computing the stability radius. A number of numerical experiments demonstrate the robustness and efficiency of the subspace methods.

Geometric Measure of Entanglement and U-Eigenvalues of Tensors

Abstract: We study tensor analysis problems motivated by the geometric measure of quantum entanglement. We define the concept of the unitary eigenvalue (U-eigenvalue) of a complex tensor, the unitary symmetric eigenvalue (US-eigenvalue) of a symmetric complex tensor, and the best complex rank-one approximation. We obtain an upper bound on the number of distinct US-eigenvalues of symmetric tensors and count all US-eigenpairs with nonzero eigenvalues of symmetric tensors. We convert the geometric measure of the entanglement problem to an algebraic equation system problem. A numerical example shows that a symmetric real tensor may have a best complex rank-one approximation that is better than its best real rank-one approximation, which implies that the absolute-value largest Z-eigenvalue is not always the geometric measure of entanglement.

A Logarithmic Minimization Property of the Unitary Polar Factor in the Spectral and Frobenius Norms

Abstract: The unitary polar factor $Q=U_p$ in the polar decomposition of $Z=U_p \, H$ is the minimizer over unitary matrices $Q$ for both $\|{\rm Log}(Q^* Z)\|^2$ and its Hermitian part $\|{{\rm sym}{_{_*}}\!}({\rm Log}(Q^* Z))\|^2$ over both $\mathbb{R}$ and $\mathbb{C}$ for any given invertible matrix $Z\in\mathbb{C}^{n\times n}$ and any matrix logarithm Log, not necessarily the principal logarithm log. We prove this for the spectral matrix norm for any $n$ and for the Frobenius matrix norm for $n\leq 3$. The result shows that the unitary polar factor is the nearest orthogonal matrix to $Z$ not only in the normwise sense but also in a geodesic distance. The derivation is based on Bhatia's generalization of Bernstein's trace inequality for the matrix exponential and a new sum of squared logarithms inequality. Our result generalizes the fact for scalars that for any complex logarithm and for all $z\in\mathbb{C}{\setminus}\{0\} \min_{\vartheta\in(-\pi,\pi]}{|{\rm Log}_{\mathbb{C}}(e^{-i\vartheta} z)|}^2={|{\rm log}|...

Orthogonal Invariance and Identifiability

Abstract: Matrix variables are ubiquitous in modern optimization, in part because variational properties of useful matrix functions often expedite standard optimization algorithms. Convexity is one important such property: permutation-invariant convex functions of the eigenvalues of a symmetric matrix are convex, leading to the wide applicability of semidefinite programming algorithms. We prove the analogous result for the property of “identifiability,” a notion central to many active-set-type optimization algorithms.

Block GMRES Method with Inexact Breakdowns and Deflated Restarting

Abstract: We consider the solution of large linear systems with multiple right-hand sides using a block GMRES approach. We introduce a new algorithm that effectively handles the situation of almost rank deficient block generated by the block Arnoldi procedure and that enables the recycling of spectral information at restart. The first feature is inherited from an algorithm introduced by Robbe and Sadkane [Linear Algebra Appl., 419 (2006), pp. 265--285], while the second one is obtained by extending the deflated restarting strategy proposed by Morgan [Appl. Numer. Math., 54 (2005), pp. 222--236]. Through numerical experiments, we show that the new algorithm combines efficiently the attractive numerical features of its two parents and outperforms them.

Deflated and Augmented Krylov Subspace Methods: A Framework for Deflated BiCG and Related Solvers

Abstract: We present an extension of the framework of Gaul et al. [SIAM J. Matrix Anal. Appl., 34 (2013), pp. 495--518] for deflated and augmented Krylov subspace methods satisfying a Galerkin condition to more general Petrov--Galerkin conditions. The main goal is to apply the framework to the biconjugate gradient method (BiCG) and some of its generalizations, including BiCGStab and IDR$(s)$. For such applications the assumptions of Gaul et al. were too restrictive. Our abstract approach does not depend on particular recurrences and thus simplifies the derivation of theoretical results. It easily leads to a variety of realizations by specific algorithms. We do not go into algorithmic details, but we show that for every method there are two different approaches for extending it by augmentation and deflation: one that explicitly takes care of the augmentation space in every step, and one that applies the appropriate Krylov subspace method to a projected problem but requires a correction step at the end. The latter ap...

On Generic Nonexistence of the Schmidt--Eckart--Young Decomposition for Complex Tensors

Abstract: The Schmidt--Eckart--Young theorem for matrices states that the optimal rank-$r$ approximation of a matrix is obtained by retaining the first $r$ terms from the singular value decomposition of that matrix. This paper considers a generalization of this optimal truncation property to the rank decomposition (Candecomp/Parafac) of tensors and establishes a necessary orthogonality condition. We prove that this condition is not satisfied at least by an open set of positive Lebesgue measure in complex tensor spaces. It is proved, moreover, that for complex tensors of small rank this condition can be satisfied only by a set of tensors of Lebesgue measure zero. Finally, we demonstrate that generic tensors in cubic tensor spaces are not optimally truncatable.

Birkhoff--von Neumann Theorem for Multistochastic Tensors

Abstract: In this paper, we study the Birkhoff--von Neumann theorem for a class of multistochastic tensors. In particular, we give a necessary and sufficient condition such that a multistochastic tensor is a convex combination of finitely many permutation tensors. It is well-known that extreme points in the set of doubly stochastic matrices are just permutation matrices. However, we find that extreme points in the set of multistochastic tensors are not just permutation tensors. We provide the other types of tensors contained in the set of extreme points.

Fast Enclosure for All Eigenvalues and Invariant Subspaces in Generalized Eigenvalue Problems

Abstract: Two fast algorithms for enclosing all eigenvalues and invariant subspaces in generalized eigenvalue problems are proposed. In these algorithms, individual eigenvectors and invariant subspaces are enclosed when eigenvalues are well separated and closely clustered, respectively. The first algorithm involves only cubic complexity and automatically determines eigenvalue clusters. The second algorithm is applicable even for defective eigenvalues. Numerical results show the properties of the proposed algorithms.