Workshop B7 - Numerical Linear Algebra

Organizers: Froilan Dopico (Universidad Carlos III de Madrid, Spain) - Alex Townsend (Massachusetts Institute of Technology, USA) - Volker Mehrmann (Technische Universität Berlin, Germany)

View workshop abstracts PDF



July 13, 14:30 ~ 15:00 - Room B5

Numerical optimization strategies for large-scale (constrained, coupled) matrix/tensor factorization

Lieven De Lathauwer

KU Leuven, Belgium   -

We give an overview of recent developments in numerical optimization-based computation of tensor decompositions. By careful exploitation of tensor product structure in methods such as quasi-Newton and nonlinear least squares, good convergence is combined with fast computation. A modular approach extends the computation to coupled factorizations and structured factors. Compact representations (polyadic, Tucker, ...) of large data sets may be obtained by stochastic optimization, randomization, compressed sensing, ... Careful exploitation of the representation structure allows us to scale the numerical algorithms to large problem size. The discussion is illustrated with application examples.

Vervliet N., Debals O., De Lathauwer L., "Tensorlab 3.0 - Numerical optimization strategies for large-scale constrained and coupled matrix/tensor factorization", in Proc. of the 50th Asilomar Conference on Signals, Systems and Computers, Pacific Grove, CA, Nov. 2016, pp. 1733-1738.

Vervliet N., Debals O., Sorber L., Van Barel M., De Lathauwer L., Tensorlab 3.0, Available online, Mar. 2016. URL:

Joint work with Nico Vervliet (KU Leuven, Belgium), Otto Debals (KU Leuven, Belgium), Laurent Sorber (KU Leuven, Belgium)(now Forespell, Belgium) and Marc Van Barel (KU Leuven, Belgium).

View abstract PDF

July 13, 15:00 ~ 15:30 - Room B5

On the complexity of solving elliptic PDEs using low-rank tensor approximations

Markus Bachmayr

Universität Bonn, Germany   -

In the construction of solvers for high-dimensional problems based on hierarchical tensor representations, one faces a trade-off between ensuring convergence and retaining small ranks of arising intermediate quantities. When approximately solving PDEs, in addition, low-rank approximation errors need to be balanced with discretization errors. We give an overview of available results on the complexity of solvers for hierarchical tensor representations, as well as of some open questions. We also consider what can or cannot be gained by tensor methods applied to different classes of high-dimensional problems.

Joint work with Albert Cohen (UPMC Paris VI), Wolfgang Dahmen (RWTH Aachen) and Reinhold Schneider (TU Berlin).

View abstract PDF

July 13, 15:30 ~ 16:00 - Room B5

Orthogonal tensors and their extremal spectral properties

André Uschmajew

University of Bonn, Germany   -

Let $\| X \|_2$, $\| X \|_F$ and $\| X \|_*$ denote the spectral, Frobenius and nuclear norm of an $n_1 \times \dots \times n_d$ tensor $X$. The sharp upper bounds for $\| X \|_F / \| X \|_2$ and $\| X \|_* / \| X \|_F$ seem to be unknown for many configurations of dimensions $(n_1,\dots,n_d)$ with $d \ge 3$. Assuming $n_1 \le \dots \le n_d$, a trivial upper bound for both ratios is $\sqrt{n_1 \cdots n_d}$. It is sharp for matrices ($d=2$) and achieved if and only if $X$ is a multiple of a matrix with pairwise orthonormal rows. Using a natural definition of orthogonal higher-order tensors, we can generalize this fact in the sense that orthogonality is necessary (up to scaling) and sufficient to achieve the trivial upper bound for both ratios. However, as it turns out, orthogonal tensors do not necessarily exist for every configuration $(n_1,\dots,n_d)$. When $d=3$, the question of existence of real orthogonal tensors admits an equivalent algebraic formulation known as Hurwitz' problem on composition formulas for bilinear forms. A classical result (due to Hurwitz) then implies that real $n \times n \times n$ orthogonal tensors only exist for $n = 1,2,4,8$. Surprisingly, the situation is different in complex spaces: unitary $n \times n \times n$ do not exist for any $n \ge 2$.

Joint work with Zhening Li (University of Portsmouth), Yuji Nakatsukasa (Oxford University) and Tasuku Soma (University of Tokyo).

View abstract PDF

July 13, 16:00 ~ 16:30 - Room B5

Polynomial root-finding using companion matrices

Fernando De Terán

Universidad Carlos III de Madrid, Spain   -

The use of companion matrices for computing the roots of scalar polynomials is a standard approach (it is, for instance, the one followed by the command 'roots' in MATLAB). It consists in computing the roots of a scalar polynomial as the eigenvalues of a companion matrix. In this talk, I will review several numerical and theoretical issues on this topic. I will pay special attention to the backward stability of solving the polynomial root-finding problem using companion matrices. More precisely, the question is the following: Even if the computed eigenvalues of the companion matrix are the eigenvalues of a nearby matrix, does this guarantee that they are the roots of a nearby polynomial? Usually, the companion matrix approach focuses on monic polynomials, since one can always divide by the leading coefficient, if necessary. But, is it enough for the backward stability issue to focus on monic polynomials? I will also pay attention to some other (more theoretical) questions like: How many companion matrices are there and what do they look like?

View abstract PDF

July 13, 17:00 ~ 17:30 - Room B5

Averaging positive definite matrices

Pierre-Antoine Absil

UCLouvain, Belgium   -

This talk concerns the problem of averaging a collection of symmetric positive-definite (SPD) matrices. Generally speaking, averaging methods are required notably to aggregate several noisy measurements of the same object, or to compute the mean of clusters in k-means clustering algorithms, or as a subtask in higher-level tasks such as curve fitting. The problem of averaging SPD matrices arises for example in medical imaging (denoising and segmentation tasks in Diffusion Tensor Imaging), mechanics (elasticity tensor computation), and in video tracking and radar detection tasks. Among several possible definitions for the mean of SPD matrices (including the straightforward arithmetic mean), the "Karcher mean" (specifically the least-squares mean in the sense of the so-called affine-invariant metric) is of widespread interest in the research literature, since it possesses several pleasant properties while being challenging to compute. In this talk, we will review recent advances in iterative methods that converge to the Karcher mean, and in methods that approach it using limited resources.

This work is in collaboration with Xinru Yuan and Kyle Gallivan (Florida State University), Wen Huang (Rice University), and Estelle Massart and Julien Hendrickx (UCLouvain).

View abstract PDF

July 13, 17:30 ~ 18:00 - Room B5

Eigenvalue Optimization with Applications to Semidefinite Programs and Rectangular Eigenvalue Problems

Emre Mengi

Koc University, Turkey   -

The majority of interest into eigenvalue optimization has originated from its intimate connection with convex semidefinite programs. This interest has diminished to a degree after the invention of interior point methods, which made numerical solutions of small convex semidefinite programs possible. Nevertheless to this day the convex semidefinite programs with large matrix variables or with large number of constraints stand as a challenge. On the other side, nonconvex eigenvalue optimization problems have been explored especially after late 1980s, the research in this direction has concentrated on particular problems for instance motivated by control theory and depending on only a few optimization variables.

This talk starts with a brief review of a unified algorithm for nonconvex eigenvalue optimization problems over a few variables. The algorithm is based on global piece-wise quadratic models for the eigenvalue functions. It works very well in practice provided bounds on second derivatives of eigenvalue functions are available analytically or heuristically.

The second part features a greedy subspace framework for eigenvalue optimization problems involving large matrices. At every iteration the subspace framework solves a projected eigenvalue optimization problem onto a small subspace, and expands the subspace with the addition of the eigenvectors at the optimal points of the small problem. The proposed framework converges superlinearly both in theory and in practice.

The third part concerns the adaption of the subspace framework for convex semidefinite programs with large matrix variables. Here we benefit from the eigenvalue optimization characterization of the dual of a semidefinite program. The small problems are solved by interior point methods. One difficulty is to start with a subspace that leads to a feasible semidefinite program, we describe an approach to overcome this difficulty.

Finally, we focus on large generalized rectangular eigenvalue problems whose spectra consist of a few imaginary eigenvalues. Such problems arise for instance from robust stability considerations of dissipative Hamiltonian systems. Extraction of their imaginary eigenvalues can be viewed as a nonconvex eigenvalue optimization problem depending on one variable, on which we employ the subspace framework. The small projected problems are nonconvex and solved by the unified algorithm.

Joint work with Nicat Aliyev (Koc University, Turkey), Fatih Kangal (Koc University, Turkey), Karl Meerbergen (KU Leuven, Belgium), Wim Michiels (KU Leuven, Belgium) and Emre Alper Yildirim (Koc University, Turkey).

View abstract PDF

July 14, 14:30 ~ 15:00 - Room B5

Revisiting the perfect shift strategy

Paul Van Dooren

Catholic University of Louvain, Belgium   -

In this paper we revisit the problem of performing a QR-step on an unreduced Hessenberg matrix when we know an exact eigenvalue of this matrix. Under exact arithmetic, this eigenvalue will appear on diagonal of the transformed Hessenberg matrix and will be decoupled from the remaining part of the Hessenberg matrix, thus resulting in a deflation. But it is well known that in finite precision arithmetic the so-called perfect shift gets blurred and the eigenvalue can not be deflated and/or is perturbed significantly. In this paper, we develop a new strategy for computing such a QR step so that the deflation is always successful.

Joint work with Nicola Mastronardi (CNR Bari, Italy).

View abstract PDF

July 14, 15:00 ~ 15:30 - Room B5

Dynamics via nonlinear pseudospectra

David Bindel

Cornell University, USA   -

Nonlinear eigenvalue problems arise naturally from transform analysis of linear differential and functional differential equations involving effects such as damping, delays, and radiation. But just as linear eigenvalue analysis tells us about asymptotic stability without necessarily providing a full picture of the intermediate transient behavior, so to does nonlinear eigenvalue analysis give an incomplete picture of the dynamics of the associated system. In this talk, we describe how generalizations of pseudospectra and related inclusion regions to the case of nonlinear eigenvalue problems, an show how they can provide information about the non-asymptotic dynamics of systems with radiation and delay.

Joint work with Amanda Hood (Cornell University).

View abstract PDF

July 14, 15:30 ~ 16:00 - Room B5

Vector Spaces of Linearizations for Rectangular Matrix Polynomials



The seminal work [MMMM06] introduced vector spaces of matrix pencils, with the property that almost all the pencils in the vector space are strong linearizations of a given square regular matrix polynomial. This work was subsequently extended to include the case of square singular matrix polynomials in [DTDM09]. Extending these ideas, we construct similar vector spaces of rectangular matrix pencils such that almost every matrix pencil of the space is a strong linearization of a given rectangular matrix polynomial $P$ in a generalized sense. Moreover, the minimal indices of $P$ can be recovered from those of the matrix pencil.

We further show that such pencils can be `trimmed' to form smaller pencils that are unimodular equivalent to $P.$ The resulting pencils are almost always strong linearizations of $P.$ Moreover they are easier to construct and are often smaller than the Fiedler linearizations of $P$ introduced in [DTDM12]. Further, the backward error analysis carried out in [DLPVD16] when applied to these trimmed linearizations shows that under suitable conditions, the computed eigenstructure of the linearizations obtained from some backward stable algorithm yield the exact eigenstructure of a slightly perturbed matrix polynomial.


[DTDM09] F. De Teran, F. M. Dopico, and D. S. Mackey, Linearizations of singular matrix polynomials and the recovery of minimal indices, Electron. J. Linear Algebra, 18(2009), pp. 371-402.

[DTDM12] F. De Teran, F. M. Dopico, and D. S. Mackey, Fiedler companion linearizations for rectangular matrix polynomials, Linear Algebra Appl., 437(2012), pp. 957-991.

[DLPVD16] F. M. Dopico, P. W. Lawrence, J. Perez and P. Van Dooren, Block Kronecker linearizations of matrix polynomials and their backward errors, MIMS Eprint, 2016.34, 2016.

[MMMM06] D. S. Mackey, N. Mackey, C. Mehl and V. Mehrmann, Vector spaces of linearizations for matrix polynomials, SIAM J. Matrix Annal. Appl., 28(4): 971-1004, 2006.

Joint work with Biswajit Das (Indian Institute of Technology Guwahati).

View abstract PDF

July 14, 16:00 ~ 16:30 - Room B5

The computational complexity of duality

Shmuel Friedland

University of Illinois at Chicago, United States   -

We show that for any given norm ball or proper cone, weak membership in its dual ball or dual cone is polynomial-time reducible to weak membership in the given ball or cone. A consequence is that the weak membership or membership problem for a ball or cone is NP-hard if and only if the corresponding problem for the dual ball or cone is NP-hard. In a similar vein, we show that computation of the dual norm of a given norm is polynomial-time reducible to computation of the given norm. This extends to convex functions satisfying a polynomial growth condition: for such a given function, computation of its Fenchel dual/conjugate is polynomial-time reducible to computation of the given function. Hence the computation of a norm or a convex function of polynomial-growth is NP-hard if and only if the computation of its dual norm or Fenchel dual is NP-hard. We discuss implications of these results on the weak membership problem for a symmetric convex body and its polar dual, the polynomial approximability of Mahler volume, and the weak membership problem for the epigraph of a convex function with polynomial growth and that of its Fenchel dual.

Joint work with Lek-Heng Lim (University of Chicago, USA).

View abstract PDF

July 14, 17:00 ~ 17:30 - Room B5

Numerical methods for Lyapunov matrix equations with banded symmetric data

Valeria Simoncini

Dipartimento di Matematica, Universita' di Bologna, Italy   -

We are interested in the numerical solution of the large-scale Lyapunov equation \[ A{\mathbf X}+{\mathbf X}A^T=C, \] where $A,C\in\cal{R}^{n\times n}$ are both large and banded matrices. We suppose that $A$ is symmetric and positive definite and $C$ is symmetric and positive semidefinite. While the case of low rank $C$ has been successfully addressed in the literature, the more general banded setting has not received much attention, in spite of its possible occurrence in applications. In this talk we aim to fill this gap.

It has been recently shown that if $A$ is well conditioned, the entries of the solution matrix $\mathbf X$ decay in absolute value as their indexes move away from the sparsity pattern of $C$. This property can be used in a memory-saving matrix-oriented Conjugate Gradient method to obtain a banded approximate solution.

For $A$ not well conditioned, the entries of $\mathbf X$ do not sufficiently decay to derive a good banded approximation. Nonetheless, we show that it is possible to split $\mathbf X$ as $\mathbf X=\mathbf Z_b + \mathbf Z_r$, where $\mathbf Z_b$ is banded and $\mathbf Z_r$ is numerically low rank. We thus propose a novel strategy that efficiently approximates both $\mathbf Z_b$ and $\mathbf Z_r$ with acceptable memory requirements.

Numerical experiments are reported to illustrate the potential of the discussed methods.

Joint work with Davide Palitta (Universita' di Bologna, Italy).

View abstract PDF

July 14, 17:30 ~ 18:00 - Room B5

On the solution of quadratic matrix equations with infinite-dimensional structured coefficients

Beatrice Meini

University of Pisa, Italy   -

Matrix equations of the kind $A_{-1}+A_0 X+A_1 X^2=X$, where both the matrix coefficients and the unknown are semi-infinite matrices belonging to a Banach algebra, are considered. These equations, where coefficients are nonnegative quasi-Toeplitz matrices, are encountered in certain Quasi-Birth-Death (QBD) stochastic processes, as the tandem Jackson queue or the reflecting random walk in the quarter plane. We provide a numerical framework for approximating the minimal nonnegative solution of these equations which relies on semi-infinite quasi-Toeplitz matrix arithmetic. This arithmetic does not perform any finite size truncation, instead it approximates the infinite matrices through the sum of a banded Toeplitz matrix and a compact correction. In particular, we show that the algorithm of Cyclic Reduction can be effectively applied and can approximate the infinite dimensional solutions with quadratic convergence at a cost which is comparable to that of the finite case. This way, we may compute a finite approximation of the sought solution, as well as of the invariant probability measure of the associated QBD process, within a given accuracy. Numerical experiments, performed on a collection of benchmarks, confirm the theoretical analysis.

Joint work with D.A. Bini (University of Pisa), S. Massei (Scuola Normale Superiore, Pisa) and L. Robol (CNR, Pisa).

View abstract PDF

July 14, 18:00 ~ 18:30 - Room B5

Recent advances in low-rank approximations

Ivan Oseledets

Skolkovo Institute of Science and Technology, Russia   -

In this talk, I will give a brief overview of our recent results in theory, algorithms, and applications of low-rank approximations of matrices and tensors. These include: 1) Jacobi-Davidson method for the optimization over low-rank matrix manifolds 2) Desingularization approach for optimization over low-rank matrix sets 3) Several new applications of tensors in chemistry, physics, machine learning.

View abstract PDF

July 15, 14:30 ~ 15:00 - Room B5

Low-rank updates of matrix functions

Daniel Kressner

EPFL, Switzerland   -

This talk is concerned with the development and analysis of fast algorithms for updating a matrix function $f(A)$ if $A$ undergoes a low-rank change $A+L$. For example, when $A$ is the Laplacian of an undirected graph then removing one edge or vertex of the graph corresponds to a rank-one change of $A$. The evaluation und update of such matrix functions (or parts thereof) is of relevance in the analysis of network community structure of networks and signal graph processing. Our algorithms are based on the tensorization of polynomial or rational Krylov subspaces involving $A$ and $A^T$. The choice of a suitable element from such a tensorized subspace for approximating $f(A+L)-f(A)$ is straightforward in the symmetric case but turns out to be more intricate in the nonsymmetric case. We show that the usual convergence results for Krylov subspace methods for matrix functions can be extended to our algorithms. If time permits, we will also touch upon the closely related problem of applying the Fréchet derivative of a matrix function to a low-rank matrix.

Joint work with Bernhard Beckermann, U de Lille 1, France and Marcel Schweitzer, EPFL, Switzerland.

View abstract PDF

July 15, 15:00 ~ 15:30 - Room B5

Subspace methods for computing the Crawford number and the real pseudospectral abscissa

Bart Vandereycken

University of Geneva, Switzerland, Switzerland   -

Certain eigenvalue or singular value optimization algorithms require repeatedly calculating the spectrum of a smoothly varying matrix. Two examples are the computation of the Crawford number and the real pseudospectral abscissa. In this talk, I will show how the computed eigenvectors and singular vectors can be used to construct subspace methods that approximate the original problems increasingly well. For the Crawford number, the subspace method is equivalent to one in [Kangal et al., A Subspace Method for Large Scale Eigenvalue Optimization, 2017] but our new convergence analysis predicts the numerical order of convergence very well. For the real pseudospectral abscissa, our algorithm is new and we prove its superlinear convergence.

Joint work with Ding Lu (University of Geneva, Switzerland) and Daniel Kressner (EPF Lausanne, Switzerland).

View abstract PDF

July 15, 15:30 ~ 16:00 - Room B5

Separable Nonnegative Matrix Factorization

Nicolas Gillis

University of Mons, Belgium   -

Nonnegative Matrix Factorization (NMF) is a linear dimensionality reduction technique for nonnegative data. NMF has become a very popular technique in data mining and machine learning because it automatically extracts meaningful features through a sparse and part-based representation. NMF consists in approximating a nonnegative data matrix with the product of two nonnegative matrices. In this talk, we first introduce NMF and illustrate its usefulness with some application examples. We then focus on the separability assumption that allows to provably solve the NMF problem (although it is NP-hard in general). We present several recent algorithms, including geometric algorithms and algorithms based on convex optimization. We illustrate the results with some numerical experiments.

Joint work with Robert Luce (EPFL, Switzerland) and Stephen Vavasis (University of Waterloo, Canada).

View abstract PDF

July 15, 16:00 ~ 16:30 - Room B5

Geometry of matrix polynomial spaces

Andrii Dmytryshyn

Umeå University, Sweden   -

We study how small perturbations of matrix polynomials may change their elementary divisors and minimal indices by constructing the closure hierarchy graphs (stratifications) of orbits and bundles of matrix polynomial Fiedler linearizations. We show that the stratification graphs do not depend on the choice of Fiedler linearization which means that all the spaces of the matrix polynomial Fiedler linearizations have the same geometry (topology). We also develop the theory for structure preserving stratification of skew-symmetric matrix polynomials. The results are illustrated by examples using the software tool Stratigraph.

Joint work with Stefan Johansson (Umeå University), Bo Kågström (Umeå University), Paul Van Dooren (Université catholique de Louvain).

View abstract PDF

July 15, 17:00 ~ 17:30 - Room B5

Parameter-Dependent Rank-One Perturbations of Singular Hermitian Pencils

Christian Mehl

TU Berlin, Germany   -

We investigate the effect of perturbations of singular Hermitian pencils $\lambda E+A$ that are (1) of rank one, (2) structure-preserving, (3) generic, (4) parameter-dependent, and (5) regularizing, i.e., generic perturbations that have the form \[ \lambda (E+\tau euu^T)+A+\tau auu^T, \] where the perturbed pencil is regular.

It is known that in this situation the eigenvalues of the perturbed pencil are constant in the parameter $\tau$. One would now expect that for a generic regularizing rank-one perturbation all eigenvalues of the perturbed pencil are simple. While this is indeed the case for pencils without symmetry structure, surprisingly this is different for real symmetric pencils where it happens that generically double eigenvalues will occur. While these eigenvalues will be constant in $\tau$ as expected, this will not be the case for their geometric multiplicities.

In this talk, we will explain why this peculiar behavior occurs and investigate what happens when the geometric multiplicity of eigenvalues changes.

Joint work with Volker Mehrmann (TU Berlin) and Michal Wojtylak (University of Krakow).

View abstract PDF

July 15, 17:30 ~ 18:00 - Room B5

Numerical instability of algebraic-geometric methods to solve polynomial equations

Vanni Noferini

University of Essex, United Kingdom   -

Among the existing algorithms to solve systems of polynomial equations, those based on algebraic geometry have particular relevance. In practice, they are often implemented either in a mixed symbolic-numeric environment, or purely in floating point arithmetic.

We focus on two important classes of methods, those based on resultant matrices and those based on Groebner bases. The common paradigm is to convert first the polynomial root-finding problem into a related eigenvalue problem, and then to find approximate solutions to the latter with standard, and reliable, numerical linear algebra algorithms. One can then recover an approximate solution to the original problem from the computed eigenvalues and eigenvectors. However, even if the conversion step is performed symbolically, the eigenvalue problem is potentially worse conditioned than the original problem. As a result, the polynomial system solver may be numerically unstable, even if the eigenvalue problem is formed without any approximation error and even it is then solved in a numerically stable manner.

This issue has been long known to practitioners. Here, we provide recent theoretical developments that quantitatively describe a worst-case spectacular increase of the condition number. We also comment on the difference between the various variants of these "numerical algebraic geometry" algorithms, discussing whether and how this effect can be mitigated when implementing these methods.

Joint work with Sujit Rao (Cornell University, United States) and Alex Townsend (Cornell University, United States).

View abstract PDF



Fast numerical solver for higher order local platforms based on error correction methods

Sunyoung Bu

Hongik University, South Korea   -

Error correction method (ECM) [1,2] is based on the construction of a local approximation to the solution on each time step. In this study, a higher-order continuous local platform is constructed to develop higher-order semi-explicit one-step ECM due to the excellent convergence order $O(h^{2p+2})$ of ECM, provided the local approximation has a local residual error $O(h^p)$. Unfortunately, the construction generates inevitably a huge size of matrix. To overcome this complexity, a fast numerical matrix solver is also provided. Special choices of parameters for the local platform can lead to the improvement of the well-known explicit fourth and fifth order Runge-Kutta methods. Numerical experiments demonstrate the theoretical results.

[1] P. Kim, X. Piao and S.D. Kim, An error corrected Euler method for solving stiff problems based on Chebyshev collocation, SIAM J. Numer. Anal., 49(2011), 2211--2230.

[2] S.D. Kim, X. Piao, D.H. Kim and P. Kim, Convergence on error correction methods for solving initial value problems, J. Comp. Appl. Math., 236(2012)(17), 4448--4461.

Joint work with Philsu Kim (Kyungpook National University, South Korea).

View abstract PDF

Perron-Frobenius Theorem for Nonnegative Tensors

Antoine Gautier

Saarland University, Germany   -

Recently, the Perron-Frobenius theory has been extended to various class of nonlinear spectral problems involving nonnegative tensors. Relevant examples include the $\ell^{p,q}$-singular vectors of matrices [1], the $\ell^p$-eigenvectors of a square tensor [2], the $\ell^{p,q}$-singular vectors of rectangular tensors [3] and the $\ell^{p_1,\ldots,p_d}$-singular vectors of tensors [4,5]. 

We present the results of [6] where a unified Perron-Frobenius theorem based on the concept of multi-homogeneous mappings is proposed. This general result can be applied to a variety of problems which includes all mentioned spectral problems for nonnegative tensors and allows to identify their common structure. It implies and improves many existing results for each of these nonnegative tensor spectral problems. In particular, when applied to these special cases, our theorem provide new characterizations of the spectral radius as well as novel convergence guarantees of the power method.

More precisely, we discuss a weak Perron-Frobenius theorem for multi-homogeneous maps which implies the existence of a nonnegative eigenvector corresponding to the spectral radius as well as a Collatz-Wielandt characterization and a number of Gelfand type formulas for the spectral radius. Furthermore, we present a strong Perron-Frobenius theorem where conditions for the existence, uniqueness and maximality of a positive eigenvector are given. We also discuss the convergence of the power-method towards the unique positive eigenvector together with a linear convergence rate. When applied to the special case of a multi-homogeneous map associated to a tensor spectral problem, our theorem allows equivalent or wider choices for $p,q$ and $p_1,\ldots,p_d$ than existing results. Moreover, the generality of our formulation allows to unify the definitions of irreducibility, weak irreducibility and weak primitivity for tensors which have been introduced separately by various authors.

[1] D. Boyd, The power method for $\ell^p$ norms, LAA 9 (1974), 95–101.

[2] L. Lim, Singular values and eigenvalues of tensors: a variational approach, CAMSAP’05, 2005, 129–132.

[3] C. Ling and L. Qi, $\ell^{k,s}$-Singular values and spectral radius of rectangular tensors, Front. Math. China 8 (2013), 63–83.

[4] S. Friedland, S. Gaubert, and L. Han, Perron-Frobenius theorem for nonnegative multilinear forms and extensions, LAA 438 (2013), 738–749.

[5] A. Gautier and M. Hein, Tensor norm and maximal singular vectors of nonnegative tensors — A Perron–Frobenius theorem, a Collatz–Wielandt characterization and a generalized power method, LAA 505 (2016), 313-343.

[6] A. Gautier, F. Tudisco and M. Hein, The Perron-Frobenius Theorem for Multi-homogeneous Maps, arXiv:1702.03230 (2017), 1–54. 

Joint work with F. Tudisco (University of Padua, Italy) and M. Hein (Saarland University, Germany).

View abstract PDF

Numerical methods for solving nonlinear matrix equations

Hyun-Min Kim

Pusan National University, Republic of Korea(South Korea)   -

We consider the existence of solutions to several types of nonlinear matrix equation, for example the quadratic matrix equations and the matrix polynomials, etc. The existences of elementwise positive solution, elementwise negative solution and positive definte solution are introduced. Also, we apply the some numerical methods to find such solutions which are Newton's method with and without incorporating exact line searches and relaxation method and functional iterations. Finally we show some numrical results.


\bibitem HHigham N. J., and Kim H.-M., {``Numerical analysis of a quadratic matrix equation''}, {\em IMA J. Numer. Anal.}, Vol. 20, 2000, pp. 499-519.

\bibitem HHigham N. J., and Kim H.-M., {``Solving a quadratic matrix equation by Newton’s method with exact line searches''}, {\em SIAM J. Matrix Anal. Appl.}, Vol. 23, 2001, pp. 303-316.

\bibitem LLatouche G., {``Newton’s iteration for nonlinear equations in Markov chains,''}, {\em IMA J. Numer. Anal.}, Vol. 14, 1994, pp. 583-598.

\bibitem SSeo J.-H., and Kim H.-M., {``Solving matrix polynomials by Newton's method with exact line searches''}, {\em J. Korean Soc. Ind. Appl. Math.}, Vol. 12, 2008, pp. 55-68.

\bibitem SSeo J.-H., and Kim H.-M., {``Convergence of pure and relaxed Newton methods for solving a matrix polynomial equation arising in stochastic models''}, {\em Linear Algebra Appl.}, Vol. 40, 2014, pp. 34-49.

\bibitem MMeng J., and Kim H.-M., {``The positive definite solution to a nonlinear matrix equation''}, {\em Linear and Multilinear Algebra}, Vol. 64, 2016, pp. 653-666.

\bibitem MMeng J., and Kim H.-M., {``The positive definite solution of the nonlinear matrix equation $X^p = A + M(B+X^{-1})^{-1} M^{\ast}$''}, {\em Journal of Computational and Applied Mathematics}, To appear. \end{thebibliography}

Joint work with Jie Meng(Pusan National University).

View abstract PDF

Stable ALS Approximation in the TT-Format for Rank-Adaptive Tensor Completion

Sebastian Kraemer

IGPM at RWTH Aachen University, Germany   -

Low rank tensor completion is a highly ill-posed inverse problem, particularly when the data model is not accurate, and some sort of regularization is required in order to solve it. In this article we focus on the calibration of the data model. For alternating optimization, we observe that existing rank adaption methods do not enable a continuous transition between manifolds of different ranks. We denote this flaw as $\textit{instability (under truncation)}$. As a consequence of this flaw, arbitrarily small changes in the singular values of an iterate can have arbitrarily large influence on the further reconstruction. We therefore introduce a singular value based regularization to the standard alternating least squares (ALS), which is motivated by averaging in micro-steps. We prove its $\textit{stability}$ and derive a natural semi-implicit rank adaption strategy. We further prove that the standard ALS micro-steps are only stable on manifolds of fixed ranks, and only around points that have what we define as $\textit{internal tensor restricted isometry property iTRIP}$. Finally, we provide numerical examples that show improvements of the reconstruction quality up to orders of magnitude in the new Stable ALS Approximation (SALSA) compared to standard ALS.

Joint work with Lars Grasedyck (IGPM at RWTH Aachen University).

View abstract PDF

An upper $J$- Hessenberg reduction of a matrix through symplectic Householder transformations


University Lille Nord de France, France   -

In this work, we introduce a reduction of a matrix to a condensed form, the upper $J$- Hessenberg form, via elementary symplectic Householder transformations, which are rank-one modification of the identity. This reduction is the crucial step for constructing an efficient SR-algorithm. The method is the analog of the reduction of a matrix to Hessenberg form, via Householder transformations, when instead of an Euclidean linear space, one takes a sympletctic one. Features of the reduction are highlighted. Two variants numerically more stables are then derived. Some numerical experiments are given, showing the efficiency of these variants.

Joint work with Haithem Ben Kahla (University Lille Nord de France, France).

View abstract PDF

A new formulation for the nearest stable matrix problem

Punit Sharma

University of Mons, Belgium   -

A square matrix is stable if all its eigenvalues are in the closed left half of the complex plane and those on the imaginary axis are semisimple. In this talk, we consider the problem of computing the nearest stable matrix to an unstable one: given a unstable matrix $A \in \mathbb R^{n,n}$, minimize the Frobenius norm of $\Delta_A \in \mathbb R^{n,n}$ such that $A+\Delta_A$ is a stable matrix. This is the converse of the stability radius problem, where a stable matrix is given and one looks for the smallest perturbation that moves an eigenvalue outside the stability region. This problem occurs for example in system identification where one needs to identify a stable system from observations.

We propose new algorithms to solve this problem based on a reformulation using linear dissipative Hamiltonian systems: we show that a matrix $A$ is stable if and only if it can be written as $A = (J-R)Q$, where $J$ is skew-symmetric, $R$ is positive semidefinite and $Q$ is positive definite. This reformulation results in an equivalent optimization problem with a simple convex feasible set. We propose three strategies to solve the problem in variables $(J,R,Q)$: (i) a block coordinate descent method, (ii) a projected gradient descent method, and (iii) a fast gradient method inspired from smooth convex optimization. These methods require $\mathcal{O}(n^3)$ operations per iteration. We show the effectiveness of the fast gradient method compared to the other approaches and to several state-of-the-art algorithms.

The paper is under review and available from \url{}.

Joint work with Nicolas Gillis (University of Mons, Belgium).

View abstract PDF

The GSVD and CSD: Matrix Trigonometry or Where are the ellipses?

Yuyang Wang

Amazon, United States   -

The SVD ellipse picture for a matrix $A$ is a very familiar visual for the action of $A$ on the unit ball. We are not aware of any ellipse pictures in the literature nor even a notion that a natural ellipse picture exists for the generalized SVD (GSVD) of two matrices $A$ and $B$ or the closely related CSD (CS Decomposition) of an orthogonal matrix. We believe that the lack of a geometric view of the GSVD is part of the reason that the GSVD is not as widely understood or as widely used as it should be. In this paper, we reveal the trigonometry, the ellipse picture, and further geometry of the GSVD. In particular, we show how the GSVD reveals whether $B$ is ``small" relative to $A$, in the same way that $B/A$ is the slope of a vector in 2d. In higher dimensions, $B$ may be small relative to $A$ in some directions, and large in others. We further characterize the connections between the generalized singular values of $(A, B)$, which we call the cotangent form of the gsvd, and the singular values of $AB^\dagger$ (e.g., $A$ times the pseudoinverse of $B$).

Joint work with Alan Edelman (MIT, USA).

View abstract PDF

FoCM 2017, based on a nodethirtythree design.