The special orthogonal group 𝕊𝕆_n is a Lie group whose geometry and local structure are encoded by the exponential map in its Lie algebra 𝐒𝐤𝐞𝐰_n, the set of skew-symmetric matrices. The associated multi-valued inverse problem – the matrix logarithm – in 𝕊𝕆_n exhibits a highly nontrivial local diffeomorphism structure, which differs from the matrix logarithm for invertible matrices. This work characterizes the local diffeomorphism structure of the exponential in the set of skew-symmetric matrices where its derivative is invertible. We show that this set with an invertible derivative can be organized into diffeomorphic regions, using a canonical alignment of Schur decompositions. In particular, the region that contains the principal logarithm has a special multiplicity structure: each matrix in 𝕊𝕆_n admits at most two skew-symmetric preimages in this region. Based on this geometric framework, we introduce the diffeomorphic logarithm of special orthogonal matrices together with an efficient and stable algorithm. Moreover, it is applied to the Karcher mean problem in 𝕊𝕆_n, demonstrating continuous behavior of the mean under perturbations of the data, which is not captured by the principal logarithm.

We propose a new method for computing the eigenvalue decomposition of a dense real normal matrix $A$ through the decomposition of its skew-symmetric part. The method relies on algorithms that are known to be efficiently implemented, such as the bidiagonal singular value decomposition and the symmetric eigenvalue decomposition. The advantages of this method stand for normal matrices with few real eigenvalues, such as random orthogonal matrices. We provide a stability and a complexity analysis of the method. The numerical performance is compared with existing algorithms. In most cases, the method has the same operation count as the Hessenberg factorization of a dense matrix. Finally, we provide experiments for the application of computing a Riemannian barycenter on the special orthogonal group.

The matrix exponential restricted to skew-symmetric matrices has numerous applications, notably in view of its interpretation as the Lie group exponential and Riemannian exponential for the special orthogonal group. We characterize the invertibility of the derivative of the skew-restricted exponential, thereby providing a simple expression of the tangent conjugate locus of the orthogonal group. In view of the skew restriction, this characterization differs from the classic result on the invertibility of the derivative of the exponential of real matrices. Based on this characterization, for every skew-symmetric matrix A outside the (zero-measure) tangent conjugate locus, we explicitly construct the domain and image of a smooth inverse – which we term nearby logarithm – of the skew-restricted exponential around A. This nearby logarithm reduces to the classic principal logarithm of special orthogonal matrices when A=0. The symbolic formulae for the differentiation and its inverse are derived and implemented efficiently. The extensive numerical experiments show that the proposed formulae are up to 3.9-times and 3.6-times faster than the current state-of-the-art robust formulae for the differentiation and its inversion, respectively.

This paper considers the problem of minimizing a differentiable function with locally Lipschitz continuous gradient on the algebraic variety of real matrices of upper-bounded rank. This problem is known to enable the formulation of various machine learning or signal processing tasks such as dimensionality reduction, collaborative filtering, and signal recovery. Several definitions of stationarity exist for this nonconvex problem. Among them, Bouligand stationarity is the strongest necessary condition for local optimality. This paper proposes two first-order methods that generate a sequence in the variety whose accumulation points are Bouligand stationary. The first method combines the well-known projected projected-gradient descent map with a rank reduction mechanism. The second method is a hybrid of projected gradient descent and projected projected-gradient descent. Both methods stand out in the field of low-rank optimization methods when considering their convergence properties, their streamlined design, their typical computational cost per iteration, and their empirically observed numerical performance. The theoretical framework used to analyze the proposed methods is of independent interest.

This paper considers the optimization problem in the form of $\min_{X \in \mathcal{F}_v} f(x) + \lambda \|X\|_1,$ where $f$ is smooth, $\mathcal{F}_v = \{X \in \mathbb{R}^{n \times q} : X^T X = I_q, v \in \mathrm{span}(X)\}$, and $v$ is a given positive vector. The clustering models including but not limited to the models used by $k$-means, community detection, and normalized cut can be reformulated as such optimization problems. It is proven that the domain $\mathcal{F}_v$ forms a compact embedded submanifold of $\mathbb{R}^{n \times q}$ and optimization-related tools including a family of computationally efficient retractions and an orthonormal basis of any normal space of $\mathcal{F}_v$ are derived. An inexact accelerated Riemannian proximal gradient method that allows adaptive step size is proposed and its global convergence is established. Numerical experiments on community detection in networks and normalized cut for image segmentation are used to demonstrate the performance of the proposed method.

Clustering is a widely used technique with a long and rich history in a variety of areas. However, most existing algorithms do not scale well to large datasets, or are missing theoretical guarantees of convergence. This paper introduces a provably robust clustering algorithm based on loss minimization that performs well on Gaussian mixture models with outliers. It provides theoretical guarantees that the algorithm obtains high accuracy with high probability under certain assumptions. Moreover, it can also be used as an initialization strategy for $k$-means clustering. Experiments on real-world large-scale datasets demonstrate the effectiveness of the algorithm when clustering a large number of clusters, and a $k$-means algorithm initialized by the algorithm outperforms many of the classic clustering methods in both speed and accuracy, while scaling well to large datasets such as ImageNet.

We consider the problem of minimizing a differentiable function with locally Lipschitz continuous gradient on a stratified set and present a first-order algorithm designed to find a stationary point of that problem. Our assumptions on the stratified set are satisfied notably by the determinantal variety (i.e., matrices of bounded rank), its intersection with the cone of positive-semidefinite matrices, and the set of nonnegative sparse vectors. The iteration map of the proposed algorithm applies a step of projected-projected gradient descent with backtracking line search, as proposed by Schneider and Uschmajew (2015), to its input but also to a projection of the input onto each of the lower strata to which it is considered close, and outputs a point among those thereby produced that maximally reduces the cost function. Under our assumptions on the stratified set, we prove that this algorithm produces a sequence whose accumulation points are stationary, and therefore does not follow the so-called apocalypses described by Levin, Kileel, and Boumal (2022). We illustrate the apocalypse-free property of our method through a numerical experiment on the determinantal variety.

We study a type of Riemannian gradient descent (RGD) algorithm, designed through Riemannian preconditioning, for optimization on $\mathcal{M}_k^{m\times n}$ -- the set of $m\times n$ real matrices with a fixed rank $k$. Our analysis is based on a quotient geometric view of $\mathcal{M}_k^{m\times n}$: by identifying this set with the quotient manifold of a two-term product space $\mathbb{R}_*^{m\times k}\times \mathbb{R}_*^{n\times k}$ of matrices with full column rank via matrix factorization, we find an explicit form for the update rule of the RGD algorithm, which leads to a novel approach to analysing their convergence behavior in rank-constrained optimization. We then deduce some interesting properties that reflect how RGD distinguishes from other matrix factorization algorithms such as those based on the Euclidean geometry. In particular, we show that the RGD algorithm is not only faster than Euclidean gradient descent but also does not rely on balancing techniques to ensure its efficiency while the latter does. We further show that this RGD algorithm is guaranteed to solve matrix sensing and matrix completion problems with linear convergence rate under the restricted positive definiteness property. Numerical experiments on matrix sensing and completion are provided to demonstrate these properties.

Limited-memory versions of quasi-Newton methods are efficient approaches to solving large-scale optimization problems in a Euclidean space. In particular, a quasi-Newton symmetric rank-one update used in a trust-region setting has proven to be an effective method. In this paper, a limited-memory Riemannian symmetric rank-one trust-region method with a restart strategy is proposed by combining the intrinsic representation of tangent vectors with a recently proposed efficient solver for the Euclidean limited-memory symmetric rank-one trust-region subproblem. The global convergence is established under the assumptions that the vector transport is isometric and the cost function is Lipschitz continuously differentiable. The computational and spatial complexity are analyzed, and detailed implementations are described. The performance of the proposed method is compared to limited-memory Riemannian BFGS method and other state-of-the-art methods using problems from matrix completion, independent component analysis, phase retrieval, and blind deconvolution. The proposed method is novel for problems on Riemannian and Euclidean spaces.

We compare two first-order low-rank optimization algorithms, namely PGD (Schneider and Uschmajew, 2015), which has been proven to be apocalypse-prone (Levin et al., 2021), and its apocalypse-free version PGDR obtained by equipping PGD with a suitable rank reduction mechanism (Olikier et al., 2022). Here an apocalypse refers to the situation where the stationarity measure goes to zero along a convergent sequence whereas it is nonzero at the limit. The comparison is conducted on two simple examples of apocalypses, the original one (Levin et al., 2021) and a new one. We also present a potential side effect of the rank reduction mechanism of PGDR and discuss the choice of the rank reduction parameter.

This paper considers the optimization problem in the form of $\min_{X \in \mathcal{F}_v} f(x) + \lambda \|X\|_1,$ where $f$ is smooth, $\mathcal{F}_v = \{X \in \mathbb{R}^{n \times q} : X^T X = I_q, v \in \mathrm{span}(X)\}$, and $v$ is a given positive vector. The clustering models including but not limited to the models used by $k$-means, community detection, and normalized cut can be reformulated as such optimization problems. It is proven that the domain $\mathcal{F}_v$ forms a compact embedded submanifold of $\mathbb{R}^{n \times q}$ and optimization-related tools including a family of computationally efficient retractions and an orthonormal basis of any normal space of $\mathcal{F}_v$ are derived. An inexact accelerated Riemannian proximal gradient method that allows adaptive step size is proposed and its global convergence is established. Numerical experiments on community detection in networks and normalized cut for image segmentation are used to demonstrate the performance of the proposed method.

In this paper, we develop a drift homotopy implicit particle filter method. The methodology of our approach is to adopt the concept of drift homotopy in the resampling procedure of the particle filter method for solving the nonlinear filtering problem, and we introduce an implicit particle filter method to improve the efficiency of the drift homotopy resampling procedure. Numerical experiments are carried out to demonstrate the effectiveness and efficiency of our drift homotopy implicit particle filter.

We consider the problem of minimizing a differentiable function with locally Lipschitz continuous gradient on the real determinantal variety, and present a first-order algorithm designed to find stationary points of that problem. This algorithm applies steps of steepest descent with backtracking line search on the variety, as proposed by Schneider and Uschmajew (2015), but by taking the numerical rank into account to perform suitable rank reductions. We prove that this algorithm produces sequences of iterates the accumulation points of which are stationary, and therefore does not follow the so-called apocalypses described by Levin, Kileel, and Boumal (2021).

Phylogenomics is a central area of research in biology that involves complicated, multi-step analyses. There are numerous tools needed to execute these analyses and many require a level of technical skills that impedes both the scope and depth of the questions asked. The combination of disjointed steps and command-line interfaces, unnecessary by modern standards, can increase the barrier to entry and labor needed to perform a full study, as well as introduce opportunities for human error. CloudForest (treescaper.github.io) is a software suite encompassing a streamlined set of phylogenomic tools, with its interface provided by the Galaxy framework. At the center of CloudForest is a custom application, TreeScaper, which performs most of the provided analyses or offloads the work to remote HPC resources. CloudForest features include phylogenomic inference via remote submission to CIPRES (phylo.org), calculation of distances between trees in a set, detection of communities of trees, and tree space visualization using non-linear dimensionality reduction. By incorporating Galaxy, CloudForest provides users with exportable workflows and histories, and access to the work of many other developers via the Galaxy Tool Shed. Visualization is included via the Galaxy visualization plugin feature. CloudForest is provided as a Docker image that may be used locally, or deployed on an organization's infrastructure with HPC backing. Support for containerized use, combined with the workflow-standardization of Galaxy, allows for improved reproducibility, speed, and accuracy, in turn encouraging collaboration, review, and the exploration of novel questions.

In this paper we analyze an indirect approach, called the Neighborhood Pattern Similarity approach, to solve the so-called role extraction problem of a large-scale graph. The method is based on the preliminary construction of a node similarity matrix which allows in a second stage to group together, with an appropriate clustering technique, the nodes that are assigned to have the same role. The analysis builds on the notion of ideal graphs where all nodes with the same role, are also structurally equivalent.

Low-rank matrix completion is the problem of recovering the missing entries of a data matrix by using the assumption that the true matrix admits a good low-rank approximation. Much attention has been given recently to exploiting correlations between the column/row entities to improve the matrix completion quality. In this paper, we propose preconditioned gradient descent algorithms for solving the low-rank matrix completion problem with graph Laplacian-based regularizers. Experiments on synthetic data show that our approach achieves significant speedup compared to an existing method based on alternating minimization. Experimental results on real world data also show that our methods provide low-rank solutions of similar quality in comparable or less time than the state-of-the-art method.

Community detection plays an important role in understanding and exploiting the structure of complex systems. Many algorithms have been developed for community detection using modularity maximization or other techniques. In this paper, we formulate the community detection problem as a constrained nonsmooth optimization problem on the compact Stiefel manifold. A Riemannian projected proximal gradient method is proposed and used to solve the problem. To the best of our knowledge, this is the first attempt to use Riemannian optimization for community detection problem. Numerical experimental results on synthetic benchmarks and real-world networks show that our algorithm is effective and outperforms several state-of-art algorithms.

A central focus of evolutionary biology is inferring the historical relationships among species and using this context to learn about how evolution has shaped diverse organisms. These historical relationships are represented by phylogenetic trees, and the methods used to infer these trees have been an active area of research for several decades. Despite this attention, phylogenetic workflows have changed little, even though extraordinary advances have occurred in the scale and pace at which genomic data have been collected in the past 20 years. Modern phylogenomic datasets have also raised fascinating new questions. Why do different parts of a genome often support different relationships among species? How are these different signals distributed across chromosomes? We developed a new computational framework, CloudForest, to tackle such questions. CloudForest is flexible, efficient, and tightly integrates a diverse set of tools. Here, we briefly describe the architecture of CloudForest, including the advantages it provides, and use it to investigate the distribution of phylogenetic signal along the entire X chromosome of 24 cat (Felidae) species.

Variation across inferred gene trees is arguably the most consistent and striking observation from empirical phylogenomic studies, yet many unanswered questions remain about the causes of this variation. One important reason these questions persist is because the field lacks robust, efficient, and reproducible workflows for investigating this variation. To meet this need, we are developing a portable cyberinfrastructure framework called CloudForest. CloudForest will provide researchers with a set of streamlined and integrated tools to explore the structure of large phylogenetic tree sets, including those generated by phylogenomic studies. CloudForest will meet many of the outstanding needs of empirical phylogenomic studies, such as (1) visualizing variation across gene trees, (2) revealing structure in sets of trees (forests), (3) conducting hypothesis tests regarding the causes of gene-tree variation, and (4) detecting genes that may have outlying (and potentially aberrant) histories. By addressing these challenges in a consistent way across computing platforms that range from a desktop computer to university high-performance computing systems to commercial cloud computing, CloudForest will allow biologists to make efficient use of any computational resources at their disposal with workflows appropriate for addressing a variety of important, unresolved questions in both evolutionary biology and other applied fields.