Differential imaging is a technique to post-process images captured by ground-based telescopes during an observation campaign, in order to make exoplanets in a distant planetary system directly visible and to remove the so-called quasi-static speckles that dramatically affect detection capabilities.In order to introduce geometric diversity between the exoplanets and the quasi-static speckles, the light is split into spectral channels during the data acquisition process, producing a 4-D data cube with images recorded at many wavelengths and at many times.In this work, we propose to follow an inverse problem approach to model the astronomical data as the contribution of a low-rank component containing the background of quasi-static speckles and a sparse component containing the exoplanets.We then formulate the resulting model as a convex non-smooth optimization model so that an accelerated proximal gradient descent can be used to solve the detection problem.

Optimization over the set of matrices $X$ that satisfy $X^\top B X = I_p$, referred to as the generalized Stiefel manifold, appears in many applications involving sampled covariance matrices such as the canonical correlation analysis (CCA), independent component analysis (ICA), and the generalized eigenvalue problem (GEVP). Solving these problems is typically done by iterative methods that require a fully formed $B$. We propose a cheap stochastic iterative method that solves the optimization problem while having access only to a random estimates of $B$. Our method does not enforce the constraint in every iteration; instead, it produces iterations that converge to critical points on the generalized Stiefel manifold defined in expectation. The method has lower per-iteration cost, requires only matrix multiplications, and has the same convergence rates as its Riemannian optimization counterparts that require the full matrix $B$. Experiments demonstrate its effectiveness in various machine learning applications involving generalized orthogonality constraints, including CCA, ICA, and the GEVP.

Orthogonality constraints naturally appear in many machine learning problems, from principal component analysis to robust neural network training. They are usually solved using Riemannian optimization algorithms, which minimize the objective function while enforcing the constraint. However, enforcing the orthogonality constraint can be the most time-consuming operation in such algorithms. Recently, Ablin and Peyré (2022) proposed the landing algorithm, a method with cheap iterations that does not enforce the orthogonality constraints but is attracted towards the manifold in a smooth manner. This article provides new practical and theoretical developments for the landing algorithm. First, the method is extended to the Stiefel manifold, the set of rectangular orthogonal matrices. We also consider stochastic and variance reduction algorithms when the cost function is an average of many functions. We demonstrate that all these methods have the same rate of convergence as their Riemannian counterparts that exactly enforce the constraint, and converge to the manifold. Finally, our experiments demonstrate the promise of our approach to an array of machine-learning problems that involve orthogonality constraints.

Effective image post-processing algorithms are vital for the successful direct imaging of exoplanets.Existing algorithms use techniques based on a low-rank approximation to separate the rotating planet signal from the quasi-static speckles.In this paper, we present a novel approach that iteratively finds the planet's flux and the low-rank approximation of quasi-static signals, strengthening the existing model based on lowrank approximations.We implement the algorithm with two different norms and test it on data, showing improvement over classical low-rank approaches.Our results highlight the benefits of iterative refinement of low-rank approximation to enhance planet detection.

We propose a direct imaging method for the detection of exo-planets based on a combined low-rank plus structured sparse model. For this task, we develop a dictionary of possible effective circular trajectories a planet can take during the observation time, elements of which can be efficiently computed using rotation and convolution operation. We design a simple alternating iterative hard-thresholding algorithm that jointly promotes a low-rank background and a sparse exoplanet foreground, to solve the non-convex optimisation problem. The experimental comparison on the β-Pictoris exoplanet benchmark dataset shows that our method has the potential to outperform the widely used Annular PCA for specific planet light intensities in terms of the receiver operating characteristic (ROC) curves.

Quality control of meteorological data is an important part of atmospheric analysis and prediction, as missing or erroneous values can degrade the quality of weather and climate information derived from these data.In practice, the presence of missing data in the weather series is quite common and problematic for many uses.We compare the performance of matrix completion methods with the state of the art to solve this missing data problem.The experimental results are carried out using the daily minimum and maximum temperature measurements of the network of weather stations operated by the Royal Meteorological Institute (RMI) of Belgium.

Direct imaging of exoplanets is a challenging task due to the small angular distance and high contrast relative to their host star, and the presence of quasi-static noise. We propose a new statistical method for direct imaging of exoplanets based on a likelihood ratio detection map, which assumes that the noise after the background subtraction step obeys a Laplacian distribution. We compare the method with two detection approaches based on signal-to-noise ratio (SNR) map after performing the background subtraction by the widely used Annular Principal Component Analysis (AnnPCA). The experimental results on the Beta Pictoris data set show the method outperforms SNR maps in terms of achieving the highest true positive rate (TPR) at zero false positive rate (FPR).

We consider model order reduction of linear time-varying systems on a finite time interval using balanced truncation. A standard way to perform MOR is to first numerically integrate the associated pair of differential Lyapunov equations for the two gramians, then compute projection matrices using the square root method, and finally formulate the reduced systems at each time instant of a chosen grid. This approach is well-known for delivering good approximation, but rather costly in computation and storage requirement. Furthermore, if one needs to compute the reduced system for any new time instant that is not included in the chosen grid, the mentioned procedure must be performed again without explicitly making use of the already computed data. For dealing with such a situation, we propose to store the projection matrices corresponding to a simplified sparse time grid and to use them to recover the projection subspaces at any other time instant via curve interpolation on the Grassmann manifold. By doing this, we can avoid the repetition of solving the differential Lyapunov equations which is the most expensive step in the procedure and therefore, as shown in a numerical example, accelerate the online reduction process.

Riemannian optimization is the task of finding an optimum of a real-valued function defined on a Riemannian manifold. Riemannian optimization has been a topic of much interest over the past few years due to many applications including computer vision, signal processing, and numerical linear algebra. The substantial background required to successfully design and apply Riemannian optimization algorithms is a significant impediment for many potential users. Therefore, multiple packages, such as Manopt (in Matlab) and Pymanopt (in Python), have been developed. This article describes ROPTLIB, a C++ library for Riemannian optimization. Unlike prior packages, ROPTLIB simultaneously achieves the following goals: (i) it has user-friendly interfaces in Matlab, Julia, and C++; (ii) users do not need to implement manifold- and algorithm-related objects; (iii) it provides efficient computational time due to its C++ core; (iv) it implements state-of-the-art generic Riemannian optimization algorithms, including quasi-Newton algorithms; and (v) it is based on object-oriented programming, allowing users to rapidly add new algorithms and manifolds.

Merging gene expression datasets is a simple way to increase the number of samples in an analysis. However experimental and data processing conditions, which are proper to each dataset or batch, generally influence the expression values and can hide the biological effect of interest. It is then important to normalize the bigger merged dataset, as failing to adjust for those batch effects may adversely impact statistical inference. Batch effect removal methods are generally based on a location-scale approach, however less widespread methods based on matrix factorization have also been proposed. We investigate on breast cancer data how those batch effect removal methods improve (or possibly degrade) the performance of simple classifiers. Our results indicate that the matrix factorization approach would deserve greater attention, as it gives results at least as good as common location-scale methods, and even significantly better results in specific cases.

We generalize the notion of Bézier surfaces and surface splines to Riemannian manifolds. To this end we put forward and compare three possible alternative definitions of Bézier surfaces. We furthermore investigate how to achieve ${\mathcal{C}}^0$- and ${\mathcal{C}}^1$-continuity of Bézier surface splines. Unlike in Euclidean space and for one-dimensional Bézier splines on manifolds, ${\mathcal{C}}^1$-continuity cannot be ensured by simple conditions on the Bézier control points: it requires an adaptation of the Bézier spline evaluation scheme. Finally, we propose an algorithm to optimize the Bézier control points given a set of points to be interpolated by a Bézier surface spline. We show computational examples on the sphere, the special orthogonal group, and two Riemannian shape spaces.

We present a new framework to fit a path to a given finite set of data points on a Riemannian manifold. The path takes the form of a continuously-differentiable concatenation of Riemannian Bezier segments. The selection of the control points that define the Bezier segments is partly guided by the differentiability requirement and by a minimal mean squared acceleration objective. We illustrate our approach on specific manifolds: the Euclidean plane (for sanity check), the sphere (as a first nonlinear illustration), the special orthogonal group (with rigid body motion applications), and the shape manifold (with 2D shape morphing applications).

It is seemingly paradoxical to the classical definition of the independent component analysis (ICA), that in reality, the true sources are often not strictly uncorrelated. With this in mind, this letter concerns a framework to extract quasi-uncorrelated sources with finite supports by optimizing a range-based contrast function under unit-norm constraints (to handle the inherent scaling indeterminacy of ICA) but without orthogonality constraints. Albeit the appealing contrast properties of the range-based function (e.g., the absence of mixing local optima), the function is not differentiable everywhere. Unfortunately, there is a dearth of literature on derivative-free optimizers that effectively handle such a nonsmooth yet promising contrast function. This is the compelling reason for the design of a nonsmooth optimization algorithm on a manifold of matrices having unit-norm columns with the following objectives: to ascertain convergence to a Clarke stationary point of the contrast function and adhere to the necessary unit-norm constraints more naturally. The proposed nonsmooth optimization algorithm crucially relies on the design and analysis of an extension of the mesh adaptive direct search (MADS) method to handle locally Lipschitz objective functions defined on the sphere. The applicability of the algorithm in the ICA domain is demonstrated with simulations involving natural, face, aerial, and texture images.

We consider the problem of fitting a discrete curve to time-labeled data points on the set Pnof all n-by-n symmetric positive-definite matrices. The quality of a curve is measured by a weighted sum of a term that penalizes its lack of fit to the data and a regularization term that penalizes speed and acceleration. The corresponding objective function depends on the choice of a Riemannian metric on Pn. We consider the Euclidean metric, the Log-Euclidean metric and the affine-invariant metric. For each, we derive a numerical algorithm to minimize the objective function. We compare these in terms of reliability and speed, and we assess the visual appear ance of the solutions on examples for n = 2. Notably, we find that the Log-Euclidean and the affine-invariant metrics tend to yield similar-and sometimes identical-results, while the former allows for much faster and more reliable algorithms than the latter.

We consider large matrices of low rank. We address the problem of recovering such matrices when most of the entries are unknown. Matrix completion finds applications in recommender systems. In this setting, the rows of the matrix may correspond to items and the columns may correspond to users. The known entries are the ratings given by users to some items. The aim is to predict the unobserved ratings. This problem is commonly stated in a constrained optimization framework. We follow an approach that exploits the geometry of the low-rank constraint to recast the problem as an unconstrained optimization problem on the Grassmann manifold. We then apply first- and second-order Riemannian trust-region methods to solve it. The cost of each iteration is linear in the number of known entries. Our methods, RTRMC 1 and 2, outperform state-of-the-art algorithms on a wide range of problem instances.

Subspace tracking methods are widespread in signal and image processing. To reduce the influence of perturbations or outliers on the measurements, some authors have used a stochastic piecewise constant velocity model on the Grassmann manifold. This paper presents an efficient way to simulate such a model using a particular representation of the Grassmann manifold. By doing so, we can reduce the spatial and time complexity of filtering techniques based on this model. We also propose an approximation of this system which can be computed in a finite number of operations and show similar results if the subspace variation is slow.

Newton's method for solving the matrix equation F(X) identical to AX-XX(T) AX = 0 runs up against the fact that its zeros are not isolated. This is due to a symmetry of F by the action of the orthogonal group. We show how differential-geometric techniques can be exploited to remove this symmetry and obtain a "geometric" Newton algorithm that finds the zeros of F. The geometric Newton method does not suffer from the degeneracy issue that stands in the way of the original Newton method.

DNA microarrays provide a huge amount of data and require therefore dimensionality reduction methods to extract meaningful biological information. Independent component analysis (ICA) was proposed by several authors as an interesting means. Unfortunately, experimental data are usually of poor quality because of noise, outliers and lack of samples. Robustness to these hurdles will thus be a key feature for an ICA algorithm. This paper identifies a robust contrast function and proposes a new ICA algorithm.

The quantity of mRNA transcripts in a cell is determined by a complex interplay of cooperative and counteracting biological processes. Independent Component Analysis (ICA) is one of a few number of unsupervised algorithms that have been applied to microarray gene expression data in an attempt to understand phenotype differences in terms of changes in the activation/inhibition patterns of biological pathways. While the ICA model has been shown to outperform other linear representations of the data such as Principal Components Analysis (PCA), a validation using explicit pathway and regulatory element information has not yet been performed. We apply a range of popular ICA algorithms to six of the largest microarray cancer datasets and use pathway-knowledge and regulatory-element databases for validation. We show that ICA outperforms PCA and clustering-based methods in that ICA components map closer to known cancer-related pathways, regulatory modules, and cancer phenotypes. Furthermore, we identify cancer signalling and oncogenic pathways and regulatory modules that play a prominent role in breast cancer and relate the differential activation patterns of these to breast cancer phenotypes. Importantly, we find novel associations linking immune response and epithelial-mesenchymal transition pathways with estrogen receptor status and histological grade, respectively. In addition, we find associations linking the activity levels of biological pathways and transcription factors (NF1 and NFAT) with clinical outcome in breast cancer. ICA provides a framework for a more biologically relevant interpretation of genomewide transcriptomic data. Adopting ICA as the analysis tool of choice will help understand the phenotype-pathway relationship and thus help elucidate the molecular taxonomy of heterogeneous cancers and of other complex genetic diseases.

Several blind source separation algorithms obtain a separating matrix by computing the congruence transformation that "best" diagonalizes a collection of covariance matrices. Recent methods avoid a pre-whitening phase and directly attempt to compute a non-orthogonal diagonalizing congruence. However, since the magnitude of the sources is unknown, there is a fundamental indeterminacy on the norm of the rows of the separating matrix. We show how this indeterminacy can be taken into account by restricting the separating matrix to the oblique manifold. The geometry of this manifold is developed and a trust-region-based algorithm for non-orthogonal joint diagonalization is proposed