I want a smooth orthogonalization process

Question

The following question is related to research I am doing on reinforcement learning on manifolds.

I have a set of basis vectors $\boldsymbol{B} = \{\boldsymbol{b}_1,\dots,\boldsymbol{b}_k\}$ that span a $k$-dimensional linear subspace of $\mathbb{R}^n$ ($\boldsymbol{b}_i \in \mathbb{R}^n$ $\forall \ i \in \{1,\dots,k\}$ and $k<n$). The basis vectors can change. I want a map that converts the basis vectors to an orthonormal basis that spans the same subspace as the original basis vectors with the requirement that the map is smooth.

I have considered the Gram–Schmidt process, but it seems to have problems. For example, consider two basis vectors in $\mathbb{R}^2$. Imagine the first vector points in the positive x-axis direction, and the second is a rotation of the first by a small angle of + or - $\epsilon$. Depending on the sign of the rotation, I will get a second orthonormal vector pointing in the direction of the either positive or negative y-axis. Therefore, slightly changing $\boldsymbol{B}$ can lead to very different orthonormal vectors. Are there other methods that don't have problems like this?

EDIT: I am looking for a practical solution to an applied problem. It is ok to assume that the basis vectors are always linearly independent and have nonzero magnitude.

Answer 1

This is not possible as you formulate it. Consider the following two bases in $\mathbb R^3$:

\begin{gather*} (e_1,\ e_1+\epsilon\, e_2) \\ (e_1,\ e_1+\epsilon\, e_3). \end{gather*}

These are close together by your definition. However, there is no orthonormal basis of the first subspace that is close to any orthonormal basis of the second subspace (in the topology of the Stiefel manifold, or in any reasonable sense of closeness).

This is because, regardless of $\epsilon$, the first basis spans $(e_1,e_2)$ and the second basis spans $(e_1, e_3)$. These are different subspaces, corresponding to different points on the Grassmanian, and so the fibers over them in the Stiefel manifold are separated by a nonzero distance.

Answer 2

You want a function $f$ that assigns to each linearly independent $k$-tuple of vectors $B$ an orthonormal $k$-tuple of vectors $f(B)$ such that these span the same subspace. As others have pointed out, you cannot expect the function $f$ to be continuous at points where it is not defined. You cannot even expect it to be uniformly continuous. Think of the case when $k=1$.

On a different note, you might find that you prefer something else to Gram-Schmidt for a different reason. Write $B$ as a $k$ by $n$ matrix. The problem is, for a rank $k$ matrix $B$, to find an invertible $k$ by $k$ matrix $M$ such that the matrix $MB$ satisfies $(MB)(MB)^t=I$. Gram-Schmidt says take $M$ to be the unique lower-triangular matrix with positive diagonal entries that satisfies this. Another way is to consider the positive-definite symmetric matrix $BB^t$ and let $M$ be the inverse of the unique positive-definite symmetric matrix $A$ such that $A^2=BB^t$. I could imagine that for your application you might prefer this to G-M, either because it tends to distort the basis less in some sense or because it works for unordered bases.

Answer 3

The Gram–Schmidt process is locally smooth, i.e. there is a small neighborhood around every basis where you don't encounter the problem you mention and where the orthonromal basis from produced depends smoothly on the input data. The problem that you describe appears only when you consider it globally, as was noted in other answers and comments.

There are other algorithms for QR decomposition apart from GS based on Householder or Givens rotations and they have some advantages (numerical stability, paralellizability). The appropriate search term would be "thin QR decomposition" or "reduced QR decomposition".

Please note that QR decompositions actually produce orthogonal basis with a rather strong property: first $l$-tuples of vectors $\{q_1, \dotsc, q_l\}$ generate the same subspace as $\{b_1, \dotsc, b_l\}$ for all $l\leq k$.

The method proposed by Tom Goodwillie is sometimes called "symmetric orthogonalization" or "SVD based orthogonalization". If you consider SVD of $B = U\Sigma V^t$, then the orthogonal basis is given by rows of $UV^t$. This solution has a remarkable property that it is the nearest "orthogonal" matrix to $B$. You can find more under the name polar decomposition.

Kind of midway between these two approaches is QR decomposition with pivoting or rank revealing QR decomposition which tackles the problem of $BB^t$ having small eigenvalues (which corresponds to the basis being close to a degenerate one).

Answer 4

I can not answer the question but maybe the following perspectives helps clarify the obstruction. Let $k \leq n$ and consider $$ \operatorname{GL}(n) \overset{\operatorname{span}}{\longrightarrow} Gr(k, n) \overset{\approx}{\longleftrightarrow} O(n) / (O(k) \times O(n - k)) \overset{\pi_{n,k}}{\longleftarrow} O(n). $$ $\operatorname{GL}(n)$ is the group of invertible matrices on $\mathbb{R}^n$, i.e. this represents your iniital bases. $\operatorname{Gr}(k, n)$ is the Grassmannian, the differentiable manifold of $k$-dimensional subspaces of $\mathbb{R}^n$. One way to construct the Grassmannian is to start with the orthogonal group $O(n)$, i.e. the orthogonal matrices on $\mathbb{R}^n$, and take the quotient by the embedding $(O(k) \times O(n - k)) \longrightarrow O(n)$ given by the matrix which applies the $O(k)$ part to the first $k$ dimensions and the $O(n - k)$ part to the remaining $n - k$ dimensions. Basically, an element of $O(n) / (O(k) \times O(n - k))$ is an equivalence classes of orthogonal matrices on $\mathbb{R}^n$ which leave invariant some $k$-dimensional subspace. What would be natural to do is to invert the projection $\pi_{n,k}$, so that $\pi_{n,k}^{-1} \circ {\approx} \circ \operatorname{span}$ gives you your desired map. Of course $\pi_{n,k}$ is not injective, so you can not invert it. But maybe one can study this map and find an atlas of natural subsets and a natural choice of $\pi_{n,k}^{-1}$ on each, somehow.

Here I don't understand well what happens so I'm guessing, but I think $O(n) / (O(k) \times O(n - k)) \overset{\pi_{n,k}}{\longleftarrow} O(n)$ is indeed a covering and therefore the universal covering space of the Grassmannian $\operatorname{Gr}(k, n)$ should give you information about the obstruction to inverting $\pi_{n,k}$. I read that the universal covering space is the "oriented Grassmannian" and it is a double cover. This would suggest that perhaps by splitting $O(n)$ into two parts (orientation). Then for each half you should be able to choose some inverse of $\pi_{n,k}$ and obtain two orthogonalization procedures, one which produces positive orientation and one which produces negative orientation bases.