Paul Breiding

Paul Breiding, Nick Vannieuwenhoven

The canonical tensor rank approximation problem (TAP) consists of approximating a real-valued tensor by one of low canonical rank, which is a challenging non-linear, non-convex, constrained optimization problem, where the constraint set forms a non-smooth semi-algebraic set. We introduce a Riemannian Gauss-Newton method with trust region for solving small-scale, dense TAPs. The novelty of our approach is threefold. First, we parametrize the constraint set as the Cartesian product of Segre manifolds, hereby formulating the TAP as a Riemannian optimization problem, and we argue why this parametrization is among the theoretically best possible. Second, an original ST-HOSVD-based retraction operator is proposed. Third, we introduce a hot restart mechanism that efficiently detects when the optimization process is tending to an ill-conditioned tensor rank decomposition and which often yields a quick escape path from such spurious decompositions. Numerical experiments show improvements of up to three orders of magnitude in terms of the expected time to compute a successful solution over existing state-of-the-art methods.

Carlos Beltrán, Paul Breiding, Nick Vannieuwenhoven

We prove the existence of an open set of $n_1\times n_2 \times n_3$ tensors of rank $r$ on which a popular and efficient class of algorithms for computing tensor rank decompositions based on a reduction to a linear matrix pencil, typically followed by a generalized eigendecomposition, is arbitrarily numerically forward unstable. Our analysis shows that this problem is caused by the fact that the condition number of the tensor rank decomposition can be much larger for $n_1 \times n_2 \times 2$ tensors than for the $n_1\times n_2 \times n_3$ input tensor. Moreover, we present a lower bound for the limiting distribution of the condition number of random tensor rank decompositions of third-order tensors. The numerical experiments illustrate that for random tensor rank decompositions one should anticipate a loss of precision of a few digits.

Paul Breiding, Nick Vannieuwenhoven

The join set of a finite collection of smooth embedded submanifolds of a mutual vector space is defined as their Minkowski sum. Join decompositions generalize some ubiquitous decompositions in multilinear algebra, namely tensor rank, Waring, partially symmetric rank and block term decompositions. This paper examines the numerical sensitivity of join decompositions to perturbations; specifically, we consider the condition number for general join decompositions. It is characterized as a distance to a set of ill-posed points in a supplementary product of Grassmannians. We prove that this condition number can be computed efficiently as the smallest singular value of an auxiliary matrix. For some special join sets, we characterized the behavior of sequences in the join set converging to the latter's boundary points. Finally, we specialize our discussion to the tensor rank and Waring decompositions and provide several numerical experiments confirming the key results.

Paul Breiding, Nick Vannieuwenhoven

We obtain estimates of the multiplicative constants appearing in local convergence results of the Riemannian Gauss-Newton method for least squares problems on manifolds and relate them to the geometric condition number of [P. Bürgisser and F. Cucker, Condition: The Geometry of Numerical Algorithms, 2013].

Paul Breiding, Nick Vannieuwenhoven

We compute the expected value of powers of the geometric condition number of random tensor rank decompositions. It is shown in particular that the expected value of the condition number of $n_1\times n_2 \times 2$ tensors with a random rank-$r$ decomposition, given by factor matrices with independent and identically distributed standard normal entries, is infinite. This entails that it is expected and probable that such a rank-$r$ decomposition is sensitive to perturbations of the tensor. Moreover, it provides concrete further evidence that tensor decomposition can be a challenging problem, also from the numerical point of view. On the other hand, we provide strong theoretical and empirical evidence that tensors of size $n_1~\times~n_2~\times~n_3$ with all $n_1,n_2,n_3 \ge 3$ have a finite average condition number. This suggests there exists a gap in the expected sensitivity of tensors between those of format $n_1\times n_2 \times 2$ and other order-3 tensors. For establishing these results, we show that a natural weighted distance from a tensor rank decomposition to the locus of ill-posed decompositions with an infinite geometric condition number is bounded from below by the inverse of this condition number. That is, we prove one inequality towards a so-called condition number theorem for the tensor rank decomposition.

Paul Breiding

Let $A \in (\mathbb{C}^{n})^{\otimes p}$ be a complex tensor of order $p$. The pair $(v,η)\in\mathbb{C}^n\times \mathbb{C}$ is called an h-eigenpair of $A$, if $v\neq0$ and it satisfies $Av^{p-1}=η^{p-2} v$, where $Av^{p-1}$ is the contraction of $A$ by $v$ in all but the first modes. We describe a randomized algorithm to compute approximations of h-eigenpairs of complex tensors. Assuming random input, the average number of arithmetic operations it performs is polynomially bounded in the input size.

Paul Breiding, Samantha Fairchild, Pierpaola Santarsiero et al.

The essential variety is an algebraic subvariety of dimension $5$ in real projective space $\mathbb R\mathrm P^{8}$ which encodes the relative pose of two calibrated pinhole cameras. The $5$-point algorithm in computer vision computes the real points in the intersection of the essential variety with a linear space of codimension $5$. The degree of the essential variety is $10$, so this intersection consists of 10 complex points in general. We compute the expected number of real intersection points when the linear space is random. We focus on two probability distributions for linear spaces. The first distribution is invariant under the action of the orthogonal group $\mathrm{O}(9)$ acting on linear spaces in $\mathbb R\mathrm P^{8}$. In this case, the expected number of real intersection points is equal to $4$. The second distribution is motivated from computer vision and is defined by choosing 5 point correspondences in the image planes $\mathbb R\mathrm P^2\times \mathbb R\mathrm P^2$ uniformly at random. A Monte Carlo computation suggests that with high probability the expected value lies in the interval $(3.95 - 0.05,\ 3.95 + 0.05)$.

Giovanni Luca Marchetti, Erin Connelly, Paul Breiding et al.

We study the critical points over an algebraic variety of an optimization problem defined by a quadratic objective that is degenerate. This scenario arises in machine learning when the dataset size is small with respect to the model, and is typically referred to as overparametrization. Our main result relates the degenerate optimization problem to a nondegenerate one via a projection. In the highly-degenerate regime, we find that a central role is played by the ramification locus of the projection. Additionally, we provide tools for counting the number of critical points over projective varieties, and discuss specific cases arising from deep learning. Our work bridges tools from algebraic geometry with ideas from machine learning, and it extends the line of literature around the Euclidean distance degree to the degenerate setting.

Paul Breiding, Taylor Brysiewicz, David K. Johnson

We introduce a low-memory framework for certifying numerical solutions to polynomial systems which uses solution iterators and spatial partitioning trees to reduce memory requirements. We provide a prototypical algorithm, analyze its complexity, and demonstrate the memory reduction on a large example.

Paul Breiding, John Cobb, Aviva K. Englander et al.

Many hypersurfaces in algebraic geometry, such as discriminants, arise as the projection of another variety. The real complement of such a hypersurface partitions its ambient space into open regions. In this paper, we propose a new method for computing these regions. Existing methods for computing regions require the explicit equation of the hypersurface as input. However, computing this equation by elimination can be computationally demanding or even infeasible. Our approach instead derives from univariate interpolation by computing the intersection of the hypersurface with a line. Such an intersection can be done using so-called pseudo-witness sets without computing a defining equation for the hypersurface - we perform elimination without actually eliminating. We implement our approach in a forthcoming Julia package and demonstrate, on several examples, that the resulting algorithm accurately recovers all regions of the real complement of a hypersurface.