Direct Search Methods on Reductive Homogeneous Spaces

TL;DR

This paper extends direct search optimization methods to reductive homogeneous spaces, enabling efficient optimization over complex manifolds like Lie groups and Grassmannians, with practical algorithms demonstrated.

Contribution

It introduces new direct search algorithms tailored for reductive homogeneous spaces, including a tangent space-based method and a tangent-free extension.

Findings

Successful implementation of probabilistic descent on reductive homogeneous spaces

Algorithms applicable to manifolds like Grassmannians and Stiefel manifolds

Potential applications in image processing and data mining

Abstract

Direct search methods are mainly designed for use in problems with no equality constraints. However, there are many instances where the feasible set is of measure zero in the ambient space and no mesh point lies within it. There are methods for working with feasible sets that are (Riemannian) manifolds, but not all manifolds are created equal. In particular, reductive homogeneous spaces seem to be the most general space that can be conveniently optimized over. The reason is that a `law of motion' over the feasible region is also given. Examples include \Rn{n} and it's linear subspaces, Lie groups and coset manifolds such as Grassmannians and Stiefel manifolds. These are important arenas for optimization, for example, in the areas of image processing and data mining. We demonstrate optimization procedures over general reductive homogeneous spaces utilizing maps from the tangent space to the manifold. A concrete implementation of the probabilistic descent direct search method is shown. This is then extended to a procedure that works solely with the manifold elements, eliminating the need for the use of the tangent space.

Tables

Table 1. Table 1: The tangent spaces, metrics and exponential maps for some real Lie groups and matrix manifolds.

Name	$ℳ$	${𝒯}_pℳ$	$h({\omega }_1,{\omega }_2)$	${\mathrm{𝐄𝐱𝐩}}_p\left(\omega \right)$	Notes
Grassmannian	$𝙶(n,k)$	$p^T\omega =0$	$\text{Tr}\left({\omega }_1^T{\omega }_2\right)$	$\left[pV\cos \left(\Theta \right)+U\sin \left(\Theta \right)\right]V^T$	$\omega =U\Theta V^T$, art:edelman-1999
Stiefel manifold	$𝚅(n,k)$	$p^T\omega =-{\omega }^{T}p$	$\text{Tr}\left({\omega }_1^T\left[I-\frac{1}{2}p{p}^T\right]{\omega }_2\right)$		$QR=p_{⟂}{p}_{⟂}^T\omega$, art:edelman-1999
General Linear Group	${\mathrm{𝙶𝙻}}^{+}\left(n\right)$	$\mathrm{𝙶𝙻}\left(n\right)$	$\text{Tr}\left({\omega }_1^T{\omega }_2\right)$	$p\exp \left(\omega \right)$	\pbox20cm$p\in {\mathrm{𝙶𝙻}}^{+}\left(n\right)$ means
$det\left(p\right)>0$
Special Linear Group	$\mathrm{𝚂𝙻}\left(n\right)$	$Tr\left(\omega \right)=0$	$\text{Tr}\left({\omega }_1^T{\omega }_2\right)$	$p\exp \left(\omega \right)$	\pbox20cm${\mathrm{𝐄𝐱𝐩}}_p\left(\omega \right)$ is not
surjective
Special Orthogonal Group	$\mathrm{𝚂𝙾}\left(n\right)$	$\omega =-{\omega }^T$	$\text{Tr}\left({\omega }_1^T{\omega }_2\right)$	$p\exp \left(\omega \right)$	See art:gallier2000 for $\exp \left(\omega \right)$
Special Euclidean Group	$\mathrm{𝚂𝙴}\left(n\right)$		$\text{Tr}\left({\omega }_1^T{\omega }_2\right)$	$p\exp \left(\omega \right)$	\pbox20cm$A=-A^T,v\in {ℝ}^n$
See art:gallier2000 for $\exp \left(\omega \right)$
\pbox20cmSymmetric Positive Definite
Group	$\mathrm{𝚂𝙿𝙳}\left(n\right)$	$\omega ={\omega }^T$	$\text{Tr}\left(p^{-1/2}{\omega }_1{p}^{-1}{\omega }_2{p}^{-1/2}\right)$	$p^{1/2}\exp \left(p^{-1/2}\omega p^{-1/2}\right)p^{1/2}$	art:pennec-2006
Unipotent Group	$\mathrm{𝚄𝙿}\left(n\right)$	$\omega =U,U_{ii}=0$	$\text{Tr}\left({\omega }_1^T{\omega }_2\right)$	$\exp \left(\omega \right)={\sum }_{k=0}^n{\omega }^k/k!$	$U$ is upper triangular

Equations

$$M=\left[\begin{array}[]{cc}O&\mathbf{v}\\ \mathbf{0}^{T}&1\end{array}\right]\mbox{,}$$

$$M_{1}=M_{2}\left[\begin{array}[]{cc}O&\mathbf{0}\\ \mathbf{0}^{T}&1\end{array}\right]\mbox{,}$$

$$\left[\mathfrak{so}(3),\mathfrak{so}(3)\right]\subseteq \mathfrak{so}(3)\text{mbox},\left[\mathfrak{so}(3),\mathfrak{t}(3)\right]\subseteq \mathfrak{t}(3)\text{mbox}and\left[\mathfrak{t}(3),\mathfrak{t}(3)\right]\subseteq \mathfrak{so}(3)\text{mbox},$$

$$\tau=\left[\begin{array}[]{cc}0_{3}&\mathbf{v}\\ \mathbf{0}^{T}&0\end{array}\right]\mbox{,}$$

$$\omega$$

$$\left[\mathfrak{h},\mathfrak{h}\right]\subseteq \mathfrak{h}\text{mbox},\left[\mathfrak{h},\mathfrak{m}\right]\subseteq \mathfrak{m}\text{mbox}and\left[\mathfrak{m},\mathfrak{m}\right]\subseteq \mathfrak{h}\text{mbox},$$

$$\omega$$

$$\left[\mathfrak{h},\mathfrak{h}\right]\subseteq \mathfrak{h}\text{mbox}and\left[\mathfrak{h},\mathfrak{m}\right]\subseteq \mathfrak{m}\text{mbox}but\left[\mathfrak{m},\mathfrak{m}\right]⊈\mathfrak{h}\text{mbox}.$$

$$M=\left[\begin{array}[]{cc}O&Q\mathbf{v}\\ \mathbf{0}^{T}&1\end{array}\right]\mbox{ and }O=U\left[\begin{array}[]{cc}W&0\\ 0&V\end{array}\right]U^{T}$$

$$N=\left[\begin{array}[]{cc}O&\mathbf{0}\\ \mathbf{0}^{T}&1\end{array}\right]\mbox{.}$$

$$q=\left[\begin{array}[]{cc}0_{n+m}&Q\mathbf{v}\\ \mathbf{0}^{T}&0\end{array}\right]\mbox{,}$$

$$\underset{p\in (G,H)}{\min }f(p)\text{mbox}subjectto\mathbf{g}(p)\le 0$$

$$\underset{\omega \in \mathfrak{m}}{\min }f\circ \exp (\omega )\text{mbox}.$$

$$\tau=\left[\begin{array}[]{cc}0_{n}&\mathbf{v}\\ \mathbf{0}^{T}&0\end{array}\right]\mbox{ and }\exp[\tau]=\left[\begin{array}[]{cc}I_{n}&\mathbf{v}\\ \mathbf{0}^{T}&1\end{array}\right]\mbox{.}$$

$$\tau_{1}+\tau_{2}=\left[\begin{array}[]{cc}0_{n}&\mathbf{v}_{1}+\mathbf{v}_{2}\\ \mathbf{0}^{T}&0\end{array}\right]\mbox{ and }\exp[\tau_{1}+\tau_{2}]=\left[\begin{array}[]{cc}I_{n}&\mathbf{v}_{1}+\mathbf{v}_{2}\\ \mathbf{0}^{T}&1\end{array}\right]\mbox{.}$$

$$\Vert \exp [\omega ]-e\Vert \le \Vert \omega \Vert e^{\Vert \omega \Vert }$$

$$S(W)=\underset{i<j}{\max }\arccos ({\mathbf{w}}_i\cdot {\mathbf{w}}_j)$$

$$\underset{W,0\le H}{\min }\Vert X-WH{\Vert }_F^2\text{mbox}subjecttoS(W)\le ϵ\text{mbox}and\Vert {\mathbf{w}}_i\Vert =1\text{mbox}.$$

$$H_0=\arg \underset{0\le H}{\min }\Vert X-W_{0}H{\Vert }_F^2\text{mbox}where{\epsilon }_0=\Vert X-W_0{H}_0{\Vert }_F^2\text{mbox}.$$

$${\mathbf{w}}_{i,j}^{\pm }=\cos ({\varphi }_{i,j})\left[\cos ({\hat{\theta }}_{i,j})\bar{\mathbf{x}}+\sin ({\hat{\theta }}_{i,j}){\mathbf{v}}_{i,j}\right]\pm \sin ({\varphi }_{i,j}){\mathbf{u}}_{i,j}$$

Full text

††journal: JOTA

1 Introduction

Direct search methods are typically designed to work in a full-dimensional subset of $\mathbb{R}^{n}$. So, if we only have (at most) inequality constraints present, these methods are potentially useful procedures to employ, e.g., art:audet-2004 ; art:audet-2006 ; art:coope-2000 ; art:dennis-2004 ; tech:finkel-2004 ; art:price-2003 . This is especially true for black-box optimization or cases where derivative information is not available. However the presence of even one nonlinear equality constraint defining a feasible set of measure zero creates difficulties. In this case the probability of a mesh point being feasible is zero. It is shown how this difficulty can be avoided in some particularly convenient settings.

The mesh adaptive direct search (MADS) algorithms art:abramson-2009 ; art:audet-2006 provide typical illustrations. The probability that a point on the mesh used in MADS also lies in the feasible set is zero. So every point that MADS considers will likely be infeasible. Filter art:audet-2004 ; art:fletcher-2002 and progressive barrier art:audet-2009 techniques may be able to alleviate this situation if one is willing to violate the equality constraints somewhat. Also, the augmented Lagrangian pattern search algorithm in art:lewis-2002 extends the method in art:conn-1991 to non-smooth problems with equality constraints. We have in mind something quite different from these methods.

A recent example of the type of problem we will be addressing is art:audet-2015 . The method presented there handles linear equalities in the MADS algorithm art:audet-2006 . The method presented here will deal with subsets of $\mathbb{R}^{n}$ that can have a much more complicated structure.

Optimization is naturally carried out in $\mathbb{R}^{n}$, over Lie groups book:Stillwell-2008 and within other geometric objects like Grassmannians art:edelman-1999 . A common theme running through these is that of reductive homogeneous (often symmetric) spaces book:Sharpe-1997 . Globally these have a simple structure and convenient ways to move over them. Exploring how to optimize within them has been a worthwhile pursuit art:edelman-1999 ; art:Smith-1994 .

Following Sharpe book:Sharpe-1997 , Euclidean geometry was generalized through two routes: Riemannian geometry and Klein geometries. Euclidean geometry has: 1) a space $\mathbb{R}^{n}$, 2) with a symmetry Lie group, the special Euclidean group, that allows us to reach any point in $\mathbb{R}^{n}$ from any other point, and 3) a metric. Riemannian geometry generalizes Euclidean geometry by dropping the symmetry Lie group condition allowing arbitrary curvature of our underlying space. Locally the manifold looks like $\mathbb{R}^{n}$. However the symmetry group provides a convenient way to navigate within a manifold. A Klein geometry generalizes the symmetry groups that act on the underlying spaces: it specifies a space and a way to move in that space.

Generalizing optimization over $\mathbb{R}^{n}$ has tended to focus on manifolds with a Riemannian structure book:absil-2008 . Importantly, Sard’s theorem tells us that equality constraints can typically be taken to define a Riemannian manifold. The problem here is that the global structure of this type of manifold can be complicated in which case just moving over them without wandering off of their surface may be difficult.

We will extend direct search methods to constrained problems when the constraints are specified by a Klein geometry. The Klein geometry will implicitly define a subset of the ambient space $\mathbb{R}^{n}$ of our optimization variables. This subset is the feasible set, possibly further constrained with inequality constraints. The advantage of using a Klein geometry is that it additionally gives us information on how to move over the feasible set.

The current paper is laid out as follows. Section 2 provides a brief overview of the necessary differential geometry. Lie groups and Lie algebras are introduced in Section 3. Section 4 illustrates our methodology through a series of examples. We return for a look at Lie groups and algebras in Section 5 before presenting the optimization procedure in Section 6. In Section 7 we investigate an optimization routine where one can remain strictly in the Lie group and avoid any use of the Lie algebras. This is perhaps the highlight of the current paper and represents an alternate way of doing optimization. Some data mining applications of the new procedures are looked at in Section 8. A discussion follows in Section 9. Not all possibilities are dealt with: the focus is on what would be considered the easiest spaces that would most likely be used in practice.

2 Manifolds

For us a manifold $\mathcal{M}$ can be thought of as a subset of $\mathbb{R}^{n}$ that looks locally like $\mathbb{R}^{m}$, $n>m$. To formalize this, let $\mathcal{I}$ be an index set and for each $i\in\mathcal{I}$ let $\mathcal{V}_{i}\subset\mathcal{M}$ be an open subset where $\bigcup_{i\in\mathcal{I}}\mathcal{V}_{i}=\mathcal{M}$. Additionally, let there be homeomorphisms $\phi_{i}:\mathcal{V}_{i}\rightarrow\mathbb{R}^{m}$. Each pair $(\mathcal{V}_{i},\phi_{i})$ is called a chart and the set $\{(\mathcal{V}_{i},\phi_{i})\}_{i\in\mathcal{I}}$ is called an atlas. We now have local coordinates etched onto $\mathcal{M}$.

When two open subsets $\mathcal{V}_{i}$ and $\mathcal{V}_{j}$ of $\mathcal{M}$ overlap there will be a homeomorphic transition function $\tau_{ij}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ on the overlap $\mathcal{V}_{i}\cap\mathcal{V}_{j}$ defined by $\tau_{ij}=\phi_{j}\circ\phi_{i}^{-1}$. The transition functions allow us to stitch together the charts into a coherent whole. This then defines a topological manifold.

Extra structure can be placed on the manifold by placing additional requirements on $\mathcal{M}$. A typical requirement is that the transition functions are differentiable. In this case we have a differentiable manifold. If the $\tau_{ij}$ are all $\mathcal{C}^{k}$, $k\geq 1$, then $\mathcal{M}$ is said to be a $\mathcal{C}^{k}$-manifold. If $\mathcal{M}$ is a $\mathcal{C}^{\infty}$-manifold then it is called smooth. Every $\mathcal{C}^{k}$-manifold can be made smooth book:hirsch-1997 .

At every point $p\in\mathcal{M}\subset\mathbb{R}^{n}$ a tangent space $\mathcal{T}_{p}\mathcal{M}$ can be attached. If an inner product $h:\mathcal{T}_{p}\mathcal{M}\times\mathcal{T}_{p}\mathcal{M}\rightarrow\mathbb{R}$ is defined for every $p\in\mathcal{M}$ and varies smoothly over $\mathcal{M}$ then we have a Riemannian manifold. The function $h(\cdot,\cdot)$ is called a metric and with it notions of lengths and angles in $\mathcal{T}_{p}\mathcal{M}$ are available over $\mathcal{M}$. We will need surprisingly little direct use of any Riemannian structure.

Manifolds often enter an optimization problem through equality constraints.

Definition 1 (Regular level sets)

Let $\mathbf{g}:\mathbb{R}^{n+m}\rightarrow\mathbb{R}^{m}$. For the set $\mathcal{M}=\{\mathbf{x}:\mathbf{g}(\mathbf{x})=\mathbf{c}\}$ assume $\nabla\mathbf{g}(\mathbf{x})$ is full rank. Then $\mathcal{M}$ is a regular level set of $\mathbf{g}(\mathbf{x})$. Additionally, the null space of $\nabla\mathbf{g}(\mathbf{x})$ coincides with the tangent space to $\mathcal{M}$.

If the function defining the equality constraints is sufficiently differentiable then by Sard’s theorem almost every level set is a (Riemannian) manifold.

Theorem 2.1 (Regular level set manifolds and Sard’s theorem book:helgason )

Let $\mathbf{g}:\mathbb{R}^{n+m}\rightarrow\mathbb{R}^{m}$ be a $\mathcal{C}^{2}$ function. Then every regular level set of $\mathbf{g}(\mathbf{x})$ is an $n$-dimensional manifold. Further, if $\mathbf{g}(\mathbf{x})$ is a $\mathcal{C}^{n+1}$ function then almost every level set is a regular level set. A regular level set acquires a metric from its embedding in $\mathbb{R}^{n+m}$ by restricting the standard Euclidean metric to the tangent space $\mathcal{T}_{p}\mathcal{M}$ given by $\mathrm{null}(\nabla\mathbf{g}(\mathbf{p}))$ for all $p\in\mathcal{M}$.

On a Riemannian manifold $\mathcal{M}$ a geodesic is the shortest path between two (nearby) points on $\mathcal{M}$ and generalizes the notion of a straight line. The exact way this is calculated for a Riemannian manifold is not important for us. What is important though is when any arbitrary point $a\in\mathcal{M}$ can be reached from every other point $b\in\mathcal{M}$ by a geodesic. In such a case the manifold $\mathcal{M}$ is said to be geodesically complete. This has many nice properties which are covered in more detail in misc:dreisigmeyer-2006 . Not all of the manifolds considered herein will be geodesically complete. When they are not considerably more effort may be required to implement a practical optimization routine.

3 Lie Groups and Lie Algebras

A Lie group $G$ is a smooth manifold with the additional structure of being a group. For $g,h\in G$ the function $\phi(g,h)=g^{-1}h$ must additionally be smooth. A matrix Lie group is defined as any subgroup of the group of general linear transformations $\mathtt{GL}(n)$ of $\mathbb{R}^{n}$, the Lie group of invertible $n$-by-$n$ matrices. We will always work with matrix Lie groups and simply say Lie group instead of matrix Lie group.

A major advantage of having the group structure is the ability to easily move over the manifold. Every element $g\in G$ defines a diffeomorphism of $G$ by the operation of left translation defined by $L_{g}(h)=g\odot h$ for every $h\in G$. Left translation has moved the element $h$ to the new point $g\odot h$ on the manifold. Often $g\odot h$ is simply matrix multiplication but this is not always the case. For instance symmetric positive definite matrices would have $g\odot h=g^{1/2}hg^{1/2}$. We will drop the $\odot$ with the understanding that the correct multiplication is used when it differs from the usual matrix multiplication.

A path in $G$ is given by $L_{h_{1}}\circ L_{h_{2}}\circ\cdots\circ L_{h_{n}}(e)=h_{1}\cdots h_{n}=L_{h_{1}\cdots h_{n}}(e)$. The path-connected identity component of a Lie group $G$ is the subgroup of those elements of the full group that can be connected to the identity element by a path. We always assume that this is the component that is being used unless stated otherwise.

Since a Lie group is a manifold it will have a tangent space $\mathcal{T}_{e}\mathcal{M}$ at the identity. $\mathcal{T}_{e}\mathcal{M}$ is called the Lie algebra $\mathfrak{g}$ of the Lie group $G$. A Lie algebra is a vector space and is a much easier object to work with in a typical optimization procedure. For example, we can lay out a mesh for MADS in a Lie algebra. For a matrix group the elements of the Lie algebra are also given by matrices. The Lie bracket on the Lie algebra is defined as $[\omega_{1},\omega_{2}]=\omega_{1}\omega_{2}-\omega_{2}\omega_{1}$ for $\omega_{1},\omega_{2}\in\mathfrak{g}$.

There is a mapping $\exp:\mathfrak{g}\rightarrow G$ called the exponential mapping. For matrix Lie groups this corresponds to the matrix exponential. The following is a key fact for Lie groups:

Theorem 3.1 (Path from the identity book:Stillwell-2008 )

Every element of a path-connected matrix Lie group is given by a finite product of exponents of elements in the corresponding Lie algebra. That is, if $g\in G$ then $g=\exp[\omega_{1}]\cdots\exp[\omega_{n}]$ for some $\omega_{i}\in\mathfrak{g}$.

The exponential map $\exp[t\omega]$, $t\in\mathbb{R}$ and $\omega\in\mathfrak{g}$, is a geodesic that sometimes, but not always, corresponds to the Riemannian geodesic. When it doesn’t correspond we will use the Lie group exponential. An important special case is when the Lie group $G$ is compact (as a manifold). Then the exponential map is also a Riemannian geodesic and the manifold is geodesically complete.

Given a Lie group $G$ and a subgroup $H\subset G$ we can form the coset $G/H$ where every element in $G/H$ is an equivalency class $gH=\{gh:h\in H\}$ for every $g\in G$. We identify $g_{1},g_{2}\in G$ as the same element if $g_{1}H=g_{2}H$ as sets. The space $G/H$ is also a smooth manifold, called the coset manifold.

4 Examples

It is important at this point to lay out what is needed to effectively perform an optimization procedure. What is it about $\mathbb{R}^{n}$ that needs to be retained and what can be thrown away? There are three vital things that need to be done:

move in a specified direction; 2. 2.

a specified distance, and; 3. 3.

evaluate a function at every point.

If we can do these we can optimize. If we can do these efficiently we can optimize efficiently.

For Riemannian geometries all of the three requirements are satisfied. We can move a specified length along a specified direction on the manifold. Having functions defined on a Riemannian manifold also causes no difficulties. The difference from $\mathbb{R}^{n}$ is that the path we travel along now, our ’direction’, is not necessarily a straight line. Instead it’s a geodesic, the straightest possible line we could hope for and still remain on the manifold. But this is the difficulty: numerically traveling along the geodesic without straying is hard.

There are some cases of a Riemannian manifold where it is easy to stay on the geodesic. $\mathbb{R}^{n}$ is the canonical example where the geodesics are simply straight lines and the Euclidean metric gives us our sense of distance. However, it is possible for other Riemannian manifolds to have closed-form solutions for their geodesics and this makes optimizing over them a relatively straightforward procedure.

If the surface we are optimizing over has a simple global structure then we can hope that performing the optimization will be easy. A particular example of this is when we know the Lie group of symmetries of the manifold. In this case the group elements specify how to move over the manifold. Below we examine this through a few examples of progressive complexity. There will be results and definitions that, while expressed in a specific example, are generic for the spaces and Lie groups we’ll consider.

4.1 The Special Euclidean Group $\mathtt{SE}(3)$

Our first case will be $\mathbb{R}^{3}$. Here the orientation preserving symmetry group is given by the special Euclidean group $\mathtt{SE}(3)$ that includes all of the rotations and translations of $\mathbb{R}^{3}$. That is, the space looks the same under any transformation in $\mathtt{SE}(3)$. The typical way we would represent these transformations is by augmenting a vector $\mathbf{x}\in\mathbb{R}^{3}$ with a final $1$ to have $\hat{\mathbf{x}}^{T}=[\ \mathbf{x}^{T}\ 1\ ]$. A matrix representation $M$ of a group element $g\in\mathtt{SE}(3)$ is then given by

[TABLE]

where $O\in\mathtt{SO}(3)$ is a $3$-by-$3$ orthogonal matrix with $\det(O)=1$ and $\mathbf{0},\mathbf{v}\in\mathbb{R}^{3}$. Then an explicit representation of the action of $g\in\mathtt{SE}(3)$ on $\mathbf{x}\in\mathbb{R}^{3}$ is given by $M\hat{\mathbf{x}}$.

With $\mathtt{SE}(3)$ we can reach any point $\mathbf{y}\in\mathbb{R}^{3}$ from any other point $\mathbf{x}\in\mathbb{R}^{3}$. However this is not unique since we can always first perform a rotation around $\mathbf{x}$ and then translate from $\mathbf{x}$ to $\mathbf{y}$. By ‘rotation around $\mathbf{x}$’ we mean translating $\mathbf{x}$ to the origin, performing a rotation and then translating back to $\mathbf{x}$. There’s no reason to consider these (infinite number of) different ways of going from $\mathbf{x}$ to $\mathbf{y}$ as being different. Instead we form equivalency classes where $M_{1}\sim M_{2}$ when

[TABLE]

for $O\in\mathtt{SO}(3)$. Notice that the matrix that multiplies $M_{2}$ on the right in (2) is also a representation of $O\in\mathtt{SO}(3)$ since the ‘extra’ $1$ is irrelevant. It is also a member of $\mathtt{SE}(3)$ where there is no translation. We see that $\mathtt{SO}(3)\subset\mathtt{SE}(3)$ and additionally that $\mathtt{SE}(3)=\mathtt{SO}(3)\ltimes\mathtt{T}(3)$, the semi-direct product of rotations and translations in $\mathbb{R}^{3}$. From this we can form the coset manifold $\mathtt{SE}(3)/\mathtt{SO}(3)=\mathtt{T}(3)\simeq\mathbb{R}^{3}$, getting back a copy of what we started with. The optimization proceeds in the coset manifold.

A Klein geometry $(G,H)$ consists of a Lie group $G$, a subgroup $H\subset G$ and a space (manifold) $\mathcal{M}=G/H$ over which $G$ acts transitively: for every $p,q\in\mathcal{M}$ there is a $g\in G$ such that $gp=q$. Alternatively, a Klein geometry can be thought of as a space and a law for moving over this space. For the Klein geometry $(\mathtt{SE}(3),\mathtt{SO}(3))$ the following conditions for the Lie algebras $\mathfrak{se}(3)=\mathfrak{so}(3)\oplus\mathfrak{t}(3)$ hold:

[TABLE]

where $[\cdot,\cdot]$ is the Lie bracket. The first condition is that $\mathfrak{so}(3)$ is itself a Lie algebra. The second condition tells us that the $\mathfrak{t}(3)$ is invariant under the actions of $\mathfrak{so}(3)$. The final condition states that the Lie bracket of two infinitesimal translations is a rotation. As these conditions are sequentially imposed we move from a homogeneous space to a reductive homogeneous space to a symmetric space.

We are still missing a sense of moving along a straight line in $\mathbb{R}^{3}$ by a specified amount. One thing that makes reductive homogeneous spaces convenient to work with is that if the Klein geometry is given by $(G,H)$, and $\mathfrak{g}=\mathfrak{h}\oplus\mathfrak{m}$ with $[\mathfrak{h},\mathfrak{h}]\subseteq\mathfrak{h}$ and $[\mathfrak{h},\mathfrak{m}]\subseteq\mathfrak{m}$, then we can ignore the action of $\mathfrak{h}$ and use only $\mathfrak{m}$ in our optimization. What this means is we can use the same exponential map defined on $\mathfrak{g}$ and simply restrict it to $\mathfrak{m}$ to find the exponential map to $G/H$. It is the reductive condition that allows this.

In the case of $\mathfrak{t}(3)$ the matrices in the algebra are of the form

[TABLE]

and $\exp[\tau]=I_{4}+\tau$. This explicitly shows that $\mathtt{T}(3)\simeq\mathbb{R}^{3}$ through the identification $\tau\leftrightarrow\mathbf{v}$. We can travel along a straight line from from the point $T_{0}\in\mathtt{T}(3)$ using $T_{0}\exp[t\tau]$, $t\in\mathbb{R}$. As $t$ increases we are traveling further along this geodesic.

A symmetric space is always geodesically complete, which does not necessarily hold for the other homogeneous spaces book:Sternberg-2012 . Optimizing over a geodesically complete space is easier since this allows one to use a global pullback procedure from the manifold to the tangent space as covered in misc:dreisigmeyer-2006 . This means we can pullback the objective function $f:\mathcal{M}\rightarrow\mathbb{R}$, and any inequality constraints, from the manifold $\mathcal{M}$ to the Lie algebra $\mathfrak{m}$ using the function composition $f\circ\exp:\mathfrak{m}\rightarrow\mathbb{R}$. Since the Lie algebra is a normed vector space and the exponential map is smooth, any convergence results about a direct search method in $\mathbb{R}^{n}$ automatically hold without modification in the Lie algebra. In this case optimizing over the underlying space becomes a standard optimization in Euclidean space apart from a more complicated exponential function.

4.2 The Grassmannian $\mathtt{G}(4,2)$

Our second example will still be a symmetric space but more abstract. Here we will look at the Grassmannian $\mathtt{G}(4,2)$ of all two-dimensional planes in $\mathbb{R}^{4}$. We start with the Lie group $\mathtt{O}(4)$, which is the symmetry group of the $3$-sphere $\mathcal{S}^{3}\subset\mathbb{R}^{4}$. A fixed plane in $\mathbb{R}^{4}$ can be specified by the subspace normal to it. The subtlety for $\mathtt{G}(4,2)$ is that we don’t want to differentiate between any orthonormal bases that span the plane or the subspace normal to it. Then $\mathtt{G}(4,2)\simeq\mathtt{O}(4)/(\mathtt{O}(2)\times\mathtt{O}(2))$. The first $\mathtt{O}(2)$ ignores any changes in the orthonormal basis for the subspace normal to the plane; the second ignores any changes in the orthonormal basis for the plane.

For $\mathtt{G}(4,2)$, $\mathfrak{o}(4)=\mathfrak{h}\oplus\mathfrak{m}$ where

[TABLE]

with $A$ and $B$ skew-symmetric, $\eta\in\mathfrak{h}$ and $\mu\in\mathfrak{m}$. Since

[TABLE]

$\mathtt{G}(4,2)$ is a symmetric space and hence geodesically complete. Then the map $\exp[\mu]$, $\mu\in\mathfrak{m}$, allows us to reach any point in $\mathtt{G}(4,2)$ from any other point. As for $\mathbb{R}^{3}$, we can travel along a geodesic from a point $G_{0}\in\mathtt{G}(4,2)$ using $G_{0}\exp[t\mu]$, $t\in\mathbb{R}$.

Here the exponential mapping from the Lie algebra to the manifold is significantly more complex than that for the $\mathbb{R}^{3}$ case art:edelman-1999 . Table 1 gives the solution for $\mathtt{G}(n,k)$ along with some other Lie groups and reductive homogeneous spaces. Finding a way to efficiently calculate the exponential map for a reductive homogeneous space is the key to being able to efficiently optimize on it. For a Klein geometry $(G,H)$ it is enough to know the exponential for the Lie group $G$ if the space $G/H$ is a reductive homogeneous space. For us this exponential map is always the matrix exponential.

4.3 The Stiefel Manifold $\mathtt{V}(4,2)$

Our next example is the Stiefel manifold $\mathtt{V}(4,2)$ of two-dimensional planes in $\mathbb{R}^{4}$ where the orthonormal bases used to span a two-dimensional plane are distinguished. The orthonormal bases that span the normal space to a plane are not important however. So $\mathtt{V}(4,2)\simeq\mathtt{O}(4)/\mathtt{O}(2)$. Now $\mathfrak{o}(4)=\mathfrak{h}\oplus\mathfrak{m}$ is given by

[TABLE]

where $A$ and $B$ are skew-symmetric, $\eta\in\mathfrak{h}$ and $\mu\in\mathfrak{m}$.

This is our first example of a reductive homogeneous space that is not a symmetric space:

[TABLE]

So we can’t say that $\mathtt{V}(4,2)$ is geodesically complete by appealing to it being a symmetric space. But it is compact and can be made Riemannian so the Hopf-Rinow theorem tells us it is geodesically complete book:helgason . The same as for $\mathbb{R}^{3}$ and $\mathtt{G}(4,2)$, a geodesic from a point $V_{0}\in\mathtt{V}(4,2)$ is given by $V_{0}\exp[t\mu]$, $t\in\mathbb{R}$, which is both a Riemannian geodesic and a matrix exponential map. Table 1 has the closed-form solution for the geodesics.

4.4 A Fixed Subspace

This is a simple example of how equality constraints can be restated as a Klein geometry. (See art:audet-2015 also for an implementation of this sort of procedure.) Suppose we have a linear constraint $A^{T}\mathbf{x}=\mathbf{0}$ where $A^{T}:\mathbb{R}^{n+m}\rightarrow\mathbb{R}^{m}$. The subspace that $\mathbf{x}\in\mathbb{R}^{n+m}$ must lie in is given by $\mathrm{null}(A^{T})$ and $A$ gives the normal directions to the manifold. Let $Q\in\mathbb{R}^{(n+m)\times n}$ be a full rank matrix such that $A^{T}Q=0$ so that Q spans the null space. The linear subspace represented by $\mathrm{span}(Q)$ defines the manifold $\mathcal{Q}$.

Assume now that we additionally have $A^{T}A=I_{m}$ and $Q^{T}Q=I_{n}$, which can always be done. The symmetry group associated with $\mathcal{Q}$ is the subgroup of $\mathtt{SE}(n+m)$ given by matrices of the form

[TABLE]

where $U=[A\ Q]\in\mathtt{SO}(n+m)$, $W\in\mathtt{SO}(m)$, $V\in\mathtt{SO}(n)$ and $\mathbf{v}\in\mathbb{R}^{n}$. It’s easy to check that $M:\mathcal{Q}\rightarrow\mathcal{Q}$ and that all such $M$ form a Lie group $G$.

Now we have a space (manifold) given by $\mathcal{Q}$ and a law for moving over $\mathcal{Q}$ given by $G$. This defines a Klein geometry. The last piece we would like is the subgroup $H\subset G$ that will stabilize $\mathbf{0}\in\mathcal{Q}$: for all $h\in H$ we must have $h\mathbf{0}=\mathbf{0}$. Then our Klein geometry can be represented as $(G,H)$. We form equivalence classes the same as in (2) with the $O$ as in (17). Let $H\subset G$ be the subgroup of all matrices of the form

[TABLE]

Then $H$ stabilizes $\mathbf{0}\in\mathcal{Q}$ and we have our Klein geometry as $(G,H)$. Our manifold is now given by $\mathcal{Q}=G/H$. The tangent vectors $q\in\mathcal{T}_{\mathbf{0}}\mathcal{Q}$ have the form

[TABLE]

and $\exp[q]=I+q$. Any point $q\in\mathcal{Q}$ can be reached from $\mathbf{0}$ by a path of the form $\exp[q_{1}]\cdots\exp[q_{n}]\mathbf{0}$. This is how a path in the subspace would be constructed during an optimization procedure.

The example is illustrative but not practical. No one in their right mind would do this particular optimization problem as above. But the methodology presented is generic for Klein geometries associated with a reductive homogeneous space. As such it does become a practical method as the Klein geometry becomes more complicated.

5 A Second Look at Lie Groups and Lie Algebras

The Klein geometries considered above can all be represented as matrices. Since they’re reductive homogeneous spaces we can remain within the space using the matrix exponential of $\exp:\mathfrak{m}\rightarrow G/H$ by restricting the exponential mapping of $\mathfrak{g}$ to the subspace $\mathfrak{m}$. For a fixed $\omega\in\mathfrak{m}$ moving in a specific direction is given by $\exp[t\omega]$, where $t\in\mathbb{R}$. It’s easier to work with the Lie algebra since it’s a vector space. For example, one can place a mesh in $\mathfrak{g}$ and perform MADS there.

Theorem 3.1 does not tell us that the exponential map is surjective but we do have that:

Theorem 5.1 (Exponential map injectivity radius book:Stillwell-2008 )

The exponential map $\exp:\mathfrak{g}\rightarrow G$ is a diffeomorphism in a neighborhood of $0\in\mathfrak{g}$ given by $\|\omega\|<\rho$ for some $\rho>0$ where $\omega\in\mathfrak{g}$. The $\rho$ is called the injectivity radius of the exponential map. Similarly, the matrix logarithm map $\log:G\rightarrow\mathfrak{g}$ is a diffeomorphism in a neighborhood of $e\in G$ given by $\|g-e\|<\epsilon$ for some $\epsilon>0$ where $g\in G$.

The $\omega_{i}$ in Theorem 3.1 would all lie within the injectivity radius of the exponential map. Determining the injectivity radius of the exponential map is important since it determines an upper bound on allowable step sizes in any optimization routine. This is not always an easy task. One may be able to appeal to the group being compact or the space being symmetric to show geodesic completeness in which case Theorem 3.1 holds with a single $\exp[\omega]$. Alternatively, one can develop ‘tricks’ on a case by case basis.

To illustrate the latter consider the special linear group $\mathtt{SL}(n)$ where $\exp:\mathfrak{sl}(n)\rightarrow\mathtt{SL}(n)$ is not surjective. A polar decomposition of a matrix $A\in\mathtt{SL}(n)$ shows that $A=QS$ where $Q\in\mathtt{SO}(n)$, $S\in\mathtt{SPD}(n)$ with $\det(S)=1$, and both $Q$ and $S$ are unique. From Table 1 we see that every $\omega\in\mathfrak{sl}(n)$ can be represented as $\alpha+\sigma$ where $\alpha\in\mathfrak{so}(n)$, $\sigma\in\mathfrak{spd}(n)$ with $\mbox{Tr}(\sigma)=0$, and $\alpha$ and $\sigma$ are also in $\mathfrak{sl}(n)$.

Let $\mathfrak{m}$ be the traceless members of $\mathfrak{spd}(n)$ and $\mathfrak{h}\equiv\mathfrak{so}(n)$. Then $\mathfrak{sl}(n)=\mathfrak{h}\oplus\mathfrak{m}$ and $\mathtt{SL}(n)/\mathtt{SO}(n)$ is a symmetric space. To pick a general member of $\mathtt{SL}(n)$ first select a $\mu\in\mathfrak{m}$ and map this to the base space by $S=\exp[\mu]$. This is a surjective map (since $\mathcal{M}$ is a symmetric space and hence geodesically complete), so we can reach any $S\in\mathtt{SPD}(n)$, $\det(S)=1$, from the identity element. At $S$ an element $\eta\in\mathfrak{h}$ is mapped into $\mathtt{SO}(n)$ by $Q=\exp[\eta]$. This is also a surjective map. By forming the product $QS$ we see that every member of $\mathtt{SL}(n)$ can be reached by this procedure.

Unfortunately this method can not be employed in general. The reason it works so well here is that when we decompose a matrix $A\in\mathtt{SL}(n)$ into the polar form $A=QS$ both $Q$ and $S$ are also in $\mathtt{SL}(n)$. So we are always working in $\mathtt{SL}(n)$ and do not have to worry about whether the product $QS$ lies in $\mathtt{SL}(n)$ or not.

6 Algorithm and Convergence

The most general optimization problem considered here is stated as

[TABLE]

with $(G,H)$ being a Klein geometry and $\mathfrak{g}=\mathfrak{h}\oplus\mathfrak{m}$. Any inequality constraints are handled purely within the DSM and will not concern us here. They are only allowed if the DSM can handle inequality constraints. In this section we will initially present a procedure that works with a generic DSM. This is then specialized to the probabilistic descent method art:gratton-2015 in Section 7. Section 7 also develops a new method of working directly on manifolds. This new methodology relies on certain properties of the probabilistic descent method.

In the simplest case of a surjective exponential map the procedure is done by laying out the mesh in the Lie subalgebra $\mathfrak{m}\subset\mathfrak{g}$. Since the Lie algebra is a vector space this is a trivial mapping. The exponential map takes the chosen mesh points and maps them into the reductive homogeneous space $\mathcal{M}\subset G$. These points are then fed back into the DSM. (See Figure 1.) Convergence results are the same as whatever are available with the chosen DSM since we have done a global pullback of the objective function (and any inequality constraints) from $\mathcal{M}$ to $\mathfrak{m}$ misc:dreisigmeyer-2006 . In this case, ignoring the inequality constraints, (20) becomes

[TABLE]

With the surjective exponential map the optimization occurs in the fixed tangent (vector) space $\mathfrak{m}$. This is why convergence results from a DSM carry over without modification. When the exponential map is not surjective more care must be taken. What we desire is to work in a fixed tangent space after a finite number of steps in the optimization. If this occurs then convergence results are again the same as the DSM employed. This situation is shown in Figure 2. The only difference is that we may move over $\mathcal{M}$, $g_{k-1}$ and $g_{k}$ may differ, after picking out new elements from $\mathfrak{m}$ and mapping them to $\mathcal{M}$. If $g_{k-1}\neq g_{k}$ then the tangent spaces differ and we essentially start a new optimization with a new mesh.

Recall that $\rho$ is the injectivity radius of the exponential map. For non-surjective exponential maps a maximum step size $0<s_{max}\leq\rho$ is used. This guarantees that every point on $\mathcal{M}$ in a neighborhood of $g_{k}$ can be reached from $g_{k}$. In order to have convergence results we also use a second radius $0<s_{fix}<s_{max}$ such that we only move to a new point $\exp[\omega]=g_{k}\in\mathcal{M}$ when the point $\omega\in\mathfrak{m}$ lies in the annulus defined by $s_{fix}$ and $s_{max}$. (See Figure 3.) Otherwise we set $g_{k}=g_{k-1}$ and work in the same tangent space.

Suppose $g_{k}$ becomes fixed after a finite number of steps in the optimization. Since we then work in a single vector space where the exponential map is diffeomorphic all the convergence results carry over from the DSM. We are now doing a pullback procedure around $g_{k}$ but only locally on $\mathcal{M}$. If $g_{k}$ never becomes fixed the optimization does not converge. This is reasonable since for every $g_{k}$ after a finite number of steps we will move to a new point on $\mathcal{M}$ that is outside a local neighborhood of $g_{k}$ determined by $s_{fix}$.

A brief note on how we find the tangent space at each point on $\mathcal{M}$. Every tangent space looks the same as $\mathfrak{m}$ attached at $e\in\mathcal{M}$. If we are at a given point $g\in\mathcal{M}$ we can work with $\mathcal{T}_{g}\mathcal{M}$ by performing the automorphism $g:\mathcal{M}\rightarrow\mathcal{M}$. That is, we first pick a point $\omega\in\mathfrak{m}$, map this to a point $\exp[\omega]\in\mathcal{M}$ and then move that point by $L_{g}(\exp[\omega]):\mathfrak{m}\rightarrow\mathcal{M}$. This is equivalent to $\exp:\mathcal{T}_{g}\mathcal{M}\rightarrow\mathcal{M}$ and is computationally convenient.

The final optimization routine is given by Procedure 1. The hardest part of the procedure can be determining what the injectivity radius is if the exponential map is not surjective. See the probabilistic descent DSM art:gratton-2015 in Procedure 2 for a concrete example.

7 Working only with the Group

The Lie algebra occupies a central position in the optimization procedure presented in Procedure 1. In this way there is no major difference between what was presented here and the method developed in misc:dreisigmeyer-2006 . What is different is the ease that the manifold can be moved over because the symmetry group of the underlying space is known. As we stated previously, the reductive condition on a Klein geometry’s space $G/H$ is what allows this.

But if the group of symmetries is known and the elements provide a way of moving over the space we may ask: Why do we need the Lie algebra? The main reason for the use of a Lie algebra is that as a vector space it is easy to lay out a mesh or pick out a random element. If this convenience can be pushed down onto the group itself then the algebra can be dispensed with. That creates an alternate way of optimizing in reductive homogeneous spaces.

Let’s examine what this means in $\mathbb{R}^{n}$. Euclidean space is different since it doesn’t matter if one looks at the tangent space or the manifold since they both are globally identical. In subsection 4.1 we saw that a point in the Lie algebra and a point on the Lie group can be represented by, respectively,

[TABLE]

For moving in the group, if we work with the Lie algebra addition is used (followed by exponentiation) while the Lie group uses multiplication. This gives

[TABLE]

(Of course normally we dispense with all of this and just use $\mathbf{v}_{1}+\mathbf{v}_{2}$ not caring if it is being done in the Lie algebra or the Lie group.) The reason the tangent space and the manifold are interchangeable here is that $\exp[\tau_{1}+\tau_{2}]=\exp[\tau_{1}]\exp[\tau_{2}]$ which is due to the fact that $[\tau_{i},\tau_{j}]=0$ for any $\tau_{i}$ and $\tau_{j}$. A Lie algebra where the Lie bracket is trivial is called abelian and has an associated commutative (abelian) connected Lie group. What we wish to do here is to work strictly with the $g=\exp[\tau]$ elements in the Lie group and dispense with the elements $\tau$ in the Lie algebra. For abelian Lie groups (algebras) such as $\mathbb{R}^{n}$ and the $n$-dimensional tori this is trivial. In fact, $\mathbb{R}$ and the $1$-dimensional torus (the circle) effectively exhaust the abelian Lie groups since any such group is the direct product of a finite number of copies of $\mathbb{R}$ and the $1$-dimensional torus.

For abelian Lie groups it’s the case that $\exp[\omega_{1}+\omega_{2}]=\exp[\omega_{1}]\exp[\omega_{2}]$ and so summation in the Lie algebra directly translates to multiplication in the Lie group. It follows that picking some point on a mesh in the Lie algebra, which is a sum, easily translates to the Lie group. This can be extended in a straightforward way whenever sums in the Lie algebra can be easily translated to products in the Lie group regardless if the group is abelian or not. The difficulties arise when this isn’t possible.

A probabilistic descent DSM art:gratton-2015 seems to be particularly well suited to the methodology we have in mind. We present an implementation for a reductive homogeneous space in Procedure 2. This relies on the Lie algebra and the exponential map.

In Procedure 2 it’s important to notice that the work is being carried out around the $g_{k}^{\prime}\in\mathcal{M}$ at each step $k$ even if the current tangent space is being fixed at $g_{k}\in\mathcal{M}$. We fix the working tangent space in order to trivially carry over convergence results. It is only moved when there’s a danger of leaving the area where the exponential map is a diffeomorphism. However, if the procedure is carried out directly on the reductive homogeneous space without reference to $\mathfrak{m}$ we can freely move over $\mathcal{M}$ knowing that at some point the tangent space could be fixed at some point if there is convergence to a solution. The only restriction is that our move is never larger than the injectivity radius of the logarithmic map $\log:\mathcal{M}\rightarrow\mathfrak{m}$. Then we have the same convergence results as if we worked with $\mathfrak{m}$.

Assume that we have a black-box matrix generator $\mathbf{B}_{\mathcal{M}}^{r}$ that can return a random $g\in\mathcal{M}$ satisfying $\|g-e\|<r$, $0<r<R$, $R$ a constant. Further, $\mathbf{B}_{\mathcal{M}}^{r}$ must be able to generate all such $g$ that satisfy the norm condition. (Note that $R$ can be the injectivity radius of the logarithmic map but this is not a strict requirement.) The probabilistic descent algorithm that works solely on $\mathcal{M}$ is presented in Procedure 3.

We still need to specify the black boxes $\mathbf{B}_{\mathcal{M}}^{r}$ in Procedure 3. The obvious thing to use is the exponential map $\exp:\mathfrak{m}\rightarrow\mathcal{M}$ and pick a random $\omega\in\mathfrak{m}$. For that we can use the inequality

[TABLE]

where $\omega\in\mathfrak{m}$ satisfies $\|\omega\|<\rho$. However, the exponential map is being used here only as a convenience and it is not crucial that this is the $\mathbf{B}_{\mathcal{M}}^{r}$ in Procedure 3. So the algorithm is still independent of $\mathfrak{m}$. (As a note, we can use the condition $\|\omega\|<\rho$ since $\exp$ and $\log$ are inverses of each other. So working only in the region of $\mathfrak{m}$ where $\exp$ is a diffeomorphism is the same as working in the region on $\mathcal{M}$ where $\log$ is a diffeomorphism.) Any other generator is suitable though. For example, the unitary, orthogonal and symplectic groups are examined in art:diaconis1987 ; art:mezzadri-2007 ; art:stewart-1980 . If we can work with orthogonal groups then the Euclidean group can also be worked with. Generating random symmetric positive definite, unipotent, and general and special linear group members are more straightforward prospects. The closed form solutions for the exponential maps of the Grassmannian, Stiefel manifold and unipotent group in Table 1 makes the exponential mapping route more computationally tractable for these manifolds. Of course the norm condition still needs to be taken into account.

One thing we need to consider is the distribution of the random matrices on the manifold $\mathcal{M}$. Looking at Figure 4 the pullback procedure is to take the function $f:\mathcal{M}\rightarrow\mathbb{R}$ and define a new function on a vector space $\mathbb{V}$ by defining $f\circ\phi:\mathbb{V}\rightarrow\mathbb{R}$ where $\phi:\mathbb{V}\rightarrow\mathcal{M}$. As long as $\phi$ is sufficiently smooth this creates no difficulties and the optimization can proceed on $\mathbb{V}$ with all of the convergence results of the employed DSM. For a DSM ‘sufficiently smooth’ will typically mean a $\mathcal{C}^{2}$, $\mathcal{C}^{1}$ or Lipschitz function that is at least locally surjective.

For the probabilistic descent method there are conditions on the distribution of the random vectors in the vector space $\mathbb{V}$. These are easily enforced, but while the optimization is carried out in $\mathbb{V}$ the final solution actually lies on $\mathcal{M}$. Given the freedom in choosing $\phi$ it’s not obvious how placing restrictions on the vectors chosen in $\mathbb{V}$ carry over to restrictions on the elements of $\mathcal{M}$. About all that can be stated is that the probability density function of choosing an element must not vanish in any direction from the identity element $e\in\mathcal{M}$.

8 Applications

The procedures developed herein have been applied to two data mining algorithms with positive results. The first data mining application is concerned with dimensionality reduction of manifold-valued data and involves optimizing over Grassmannians misc:dreisigmeyer-2017a . Also, the objective function is not differentiable and a derivative-free optimization method is the only appropriate procedure.

The second application is an extension of the semi-nonnegative matrix factorization and optimizes over hyper-spheres misc:dreisigmeyer-2017b . This additionally involves the use of optional search steps. Equality constraints on vector lengths make the restriction to a manifold implicit. These are naturally enforced by employing the exponential and logarithmic maps from the tangent space to the manifold and vice versa. That is, at the expense of a more complicated objective function the difficulties of working with equality constraints are removed. There is also an inequality constraint on the hyper-area defined by the columns of a matrix and it is convenient to treat this as black-box.

Both of the data mining applications have optimization problems that are naturally stated on manifolds. The exponential and logarithmic maps were used rather than finding points on the manifolds directly. Modifying the solution method for the modified semi-nonnegative matrix factorization in misc:dreisigmeyer-2017b to work solely on the hypersphere $\mathcal{S}^{k-1}$, $k\in\mathbb{N}$, is relativity straightforward. We show how this can be done presently.

Let $X\in\mathbb{R}^{n\times m}$ be the data matrix where there are $m$ data points in $\mathbb{R}^{n}$, and choose a $k\in\mathbb{N}$. Typically $k\ll n$. We’re looking for a decomposition $X\approx WH$ where $W\in\mathbb{R}^{n\times k}$ will have normalized columns, $\|\mathbf{w}_{i}\|=1$ for $i=1,\ldots,k$, and $H\in\mathbb{R}^{k\times m}$ will have nonnegative entries, denoted by $0\leq H$. The measure of model tightness is the maximum geodesic distance

[TABLE]

over $\mathcal{S}^{n-1}$ between any two columns of $W$. The optimization problem is

[TABLE]

Take $\widehat{X}$ to be the data set $X$ with columns normalized to unit length and $\bar{\mathbf{x}}$ be the Karcher mean of $\widehat{X}$ proc:krakowski-2007 . Find an initial $W_{0}=[\mathbf{w}_{0,1}|\cdots|\mathbf{w}_{0,k}]$ where $S(W_{0})\leq\epsilon$. $H_{0}$ is found by solving the nonnegative least squares (NNLS) problem

[TABLE]

Set

• •

$i_{max}$ as the number of iterations to run the algorithm and initialize $i=0$,

• •

$\alpha_{max}=1$ as the maximum step size,

• •

$\alpha_{0}=\alpha_{max}$ as the initial step size,

• •

$\theta=1/2$ as the step size decrease,

• •

$\gamma=2$ as the step size increase, and

• •

$\rho(\alpha)=10^{-3}\alpha^{2}$ as the forcing function.

At each step of the algorithm the vectors $\bar{\mathbf{x}}$ and $\mathbf{w}_{i,j}$, $j=1,\ldots,k$, determine a two-dimensional subspace. There is a normalized vector $\mathbf{v}_{i,j}$ in this subspace that satisfies $\bar{\mathbf{x}}\cdot\mathbf{v}_{i,j}=0$, and an angle $0<\theta_{i,j}<\pi/2$ such that $\mathbf{w}_{i,j}=\cos(\theta_{i,j})\bar{\mathbf{x}}+\sin(\theta_{i,j})\mathbf{v}_{i,j}$. A contraction of $\mathbf{w}_{i,j}$ is a reduction in the value of $\theta_{i,j}$ by some multiplicative factor $0<\delta_{i,j}<1$. [Below we’ll use $\delta=\delta_{i,j}\equiv 0.99$.] A dilation of $\mathbf{w}_{i,j}$ is an increase in the value of $\theta_{i,j}$ by some multiplicative factor $1<\gamma_{i,j}$. Consider an additional normalized random vector $\mathbf{u}_{i,j}$ where $\bar{\mathbf{x}}\cdot\mathbf{u}_{i,j}=\mathbf{u}_{i,j}\cdot\mathbf{v}_{i,j}=0$, and an angle $0<\phi_{i,j}<\pi/2$. A perturbation of $\mathbf{w}_{i,j}$ is given by

[TABLE]

where $\hat{\theta}_{i,j}$ will be in some neighborhood of $\theta_{i,j}$.

The optimization routine is as follows:

Step 1

Set $i\leftarrow i+1$.

Step 2 (Optional Search Steps)

Try the following two searches.

Step 2a (optimal solution with contraction)

Set $W^{\prime}=XH^{\dagger}_{i-1}$, $H^{\dagger}_{i-1}$ the pseudo-inverse of $H_{i-1}$, and normalize the columns to unit length. Until $S(W^{\prime})\leq\epsilon$ and $0<\theta_{j}^{\prime}<\pi/2$, iteratively contract the $\mathbf{w}^{\prime}_{j}$ by a factor of $\delta=0.99$. Then find $H^{\prime}$ by solving the NNLS problem. If $\varepsilon_{i-1}-\varepsilon^{\prime}>\rho(\alpha_{i-1})$ set $\varepsilon_{i}\leftarrow\varepsilon^{\prime}$, $W_{i}\leftarrow W^{\prime}$, $H_{i}\leftarrow H^{\prime}$, $\alpha_{i}=\min(\alpha_{max},\gamma\alpha_{i-1})$ and Goto Step 5.

Step 2b (dilation)

Form $\mathbf{w}^{\prime}_{j}$ by dilating $\mathbf{w}_{i-1,j}$ by a random factor $1<\gamma^{\prime}_{j}<1+\alpha_{i-1}$ under the restriction $0<\theta_{j}^{\prime}<\pi/2$. If $S(W^{\prime})\leq\epsilon$ and $\varepsilon_{i-1}-\varepsilon^{\prime}>\rho(\alpha_{i-1})$ set $\varepsilon_{i}\leftarrow\varepsilon^{\prime}$, $W_{i}\leftarrow W^{\prime}$, $H_{i}\leftarrow H^{\prime}$, $\alpha_{i}=\min(\alpha_{max},\gamma\alpha_{i-1})$ and Goto Step 5.

Step 3 (Poll Steps)

Pick random angles $0<\theta_{j}^{\prime}<\pi/2$ where $|\theta_{j}^{\prime}-\theta_{i-1,j}|<\alpha_{i-1}$, and $0<\phi_{j}^{\prime}<\alpha_{i-1}$.

Step 3a

Form $\mathbf{w}^{\prime+}_{j}$ by perturbing $\mathbf{w}_{i-1,j}$. If $S(W^{\prime})\leq\epsilon$ and $0<\arccos(\mathbf{w}^{\prime+}_{j}\cdot\bar{\mathbf{x}})<\pi/2$ then find $H^{\prime}$ and $\varepsilon^{\prime}$ by solving the NNLS problem. If $\varepsilon_{i-1}-\varepsilon^{\prime}>\rho(\alpha_{i-1})$ set $\varepsilon_{i}\leftarrow\varepsilon^{\prime}$, $W_{i}\leftarrow W^{\prime}$, $H_{i}\leftarrow H^{\prime}$, $\alpha_{i}=\min(\alpha_{max},\gamma\alpha_{i-1})$ and Goto Step 5.

Step 3b

Form $\mathbf{w}^{\prime-}_{j}$ by perturbing $\mathbf{w}_{i-1,j}$. If $S(W^{\prime})\leq\epsilon$ and $0<\arccos(\mathbf{w}^{\prime-}_{j}\cdot\bar{\mathbf{x}})<\pi/2$ then find $H^{\prime}$ and $\varepsilon^{\prime}$ by solving the NNLS problem. If $\varepsilon_{i-1}-\varepsilon^{\prime}>\rho(\alpha_{i-1})$ set $\varepsilon_{i}\leftarrow\varepsilon^{\prime}$, $W_{i}\leftarrow W^{\prime}$, $H_{i}\leftarrow H^{\prime}$, $\alpha_{i}=\min(\alpha_{max},\gamma\alpha_{i-1})$ and Goto Step 5.

Step 4

Set $\epsilon_{i}\leftarrow\epsilon_{i-1}$, $W_{i}\leftarrow W_{i-1}$, $H_{i}\leftarrow H_{i-1}$ and $\alpha_{i}=\theta\alpha_{i-1}$.

Step 5

If $i=i_{max}$ return $W_{i}$ and $H_{i}$, otherwise Goto Step 1.

In the above method the tangent space to $\mathcal{S}^{k-1}$ is never used. Rather, we move directly over $\mathcal{S}^{k-1}$ or pick new points on $\mathcal{S}^{k-1}$. This avoids the need for an exponential (logarithmic) map from (to) the tangent space to (from) the manifold. The search steps also demonstrate the importance of using the probabilistic descent method in art:gratton-2015 . When we let $W=XH^{\dagger}$ and normalize the columns of $W$, the results of mapping these points on $\mathcal{S}^{k-1}$ to the tangent space will most likely not lie on any a priori existing mesh.

9 Conclusions

If an optimization problem has a feasible set with enough symmetry then it may be possible to find a ‘law of motion’ over the set. This will have a group structure as well as a notion of distance. These two requirements, along with an initial feasible point, are all that are needed for performing an optimization routine. Using the group structure the initial feasible point can be moved over the manifold a specified distance. This movement allows us to (locally) reconstruct the entire feasible set. However a difficulty with direct search methods like MADS is the need for a mesh. The mesh implicitly restricts the group elements and the freedom to reconstruct the feasible set.

When working with nonlinear manifolds, even those with a high degree of symmetry, placing an appropriate mesh on the feasible set can be a daunting challenge. One way to avoid the problem is to have the mesh be in a tangent space to the manifold. This does work but the underlying symmetries of the feasible set are lost. The probabilistic descent methodart:gratton-2015 removes the need for a mesh. Because of this one is now free to move over the manifold using the ‘law of motion’ with only a restriction on the distance traveled. Removing the reliance on a mesh allows the symmetries of the feasible set to be fully utilized.

Disclaimer Any opinions and conclusions expressed herein are those of the author and do not necessarily represent the views of the U.S. Census Bureau. The research in this paper does not use any confidential Census Bureau information. This was authored by an employee of the US national government. As such, the Government retains a nonexclusive, royalty-free right to publish or reproduce this article, or to allow others to do so, for Government purposes only.

Acknowledgements The author would like to thank the editor and reviewers for their helpful comments and suggestions.