Numerical methods for tensor decompositions. From least squares to beta-divergence, from batch to updating.

Over the years, many algorithms have been designed to decompose a matrix into a product of other matrices. These matrix decompositions can be used to compress data with a minimal loss of information or for extracting meaningful components. More recently, tensor decompositions such as the canonical
polyadic decomposition (CPD) and the low multilinear rank approximation (LMLRA) have been designed as higher-order generalizations of these matrix decompositions leading to better results in certain applications. While matrix methods offer a natural way to process data involving two variables, they
often fail to recognize higher-order patterns. Tensors, in contrast, can have any number of modes and can preserve the higher-order structure of the data, which makes them highly useful to analyze multi-variable data. Second, tensor decompositions are often unique in cases where their matrix counterparts are
not. As fewer (artificial) restrictions such as orthogonality or triangle structure have to be imposed on the terms to enforce their uniqueness, this can lead to a decomposition in more meaningful terms.

Current tensor decomposition methods commonly make two assumptions about the data supplied by the user. First, they assume that all data is available at once and second, they consider noise on the data to be additive and Gaussian distributed. In many applications, these assumptions are justified and one can decompose the tensor using standard batch methods that minimize the least-squares distance between the tensor and the low-rank model. However, as tensor methods are applied to more and more real-world problems, the number of cases where these assumptions are clearly violated increases as well. Any real-time application, for instance, violates the assumption of having all data available at the start. In such applications, such as process control or patient monitoring, new data arrive at certain time intervals. This data should be included in the tensor model without recomputing the full decomposition from scratch, as such approaches are not time- nor resource-efficient. Similarly, it is not hard to find applications where the assumption of additive Gaussian noise is not satisfied. Audio data, for instance, is generally modelled using a non-Gaussian noise distribution while low-light imaging assumes Poisson distributed pixel intensities. When the real and imaginary parts of a signal are collected independently, as is typically the case in MRI imaging, the noise distribution of the signal magnitude is Rician instead of Gaussian. For these applications, fitting a low-rank model to the tensor with the least-squares distance as cost function will lead to suboptimal solutions while a suitable choice of an alternative cost function can provide a model that corresponds better to the underlying components.

In this thesis, we provide methods that can still efficiently compute a suitable tensor decomposition when either of the aforementioned assumptions is not satisfied. First, we propose updating methods for both the CPD and the LMLRA. These methods start from an existing tensor decomposition and efficiently update this decomposition when new data arrives. By exploiting the multilinear structure of the tensor models, these methods are both efficient and accurate in tracking low-rank representations of the data. When desired, one can even use a weighted least-squares approach to weight certain tensor entries more heavily than others, e.g., when the uncertainty on the collected tensor entries varies between sensors. The performance of the CPD method is demonstrated for the monitoring of brain hemodynamics coupling in neonates. A second contribution is that we extend the existing least-squares CPD algorithms to other cost functions. This makes it possible to use a dedicated cost function for the application at hand and to obtain better solutions than when standard least-squares approaches are used. We focus in particular on the class of beta-divergence cost functions, of which the least-squares distance and the Kullback—Leibler and Itakura—Saito divergences are special cases, but we also supply a dedicated method for the Rician cost function. Additionally, any twice-differentiable entry-wise cost function is supported. We demonstrate the improved effectiveness of the method over least-squares approaches by blindly separating a series of fMRI signals. Finally, we propose an improved algebraic algorithm to compute the CPD of a tensor in the least-squares sense. By recursively splitting the generalized eigenspace into successively smaller subspaces instead of determining all generalized eigenvalues at once, the decomposition is more accurate, yet can still be obtained efficiently. The algebraic solution can be used as an initialization for more costly optimization algorithms and significantly reduce the number of iterations and, hence, also the runtime of these methods.