The Procrustes Problem for Orthogonal Stiefel Matrices.
In this paper we consider the Procrustes problem on the manifold of orthogonal Stiefel matrices. Given matrices ${\cal A}\in {\Bbb R}^{m\times k},$ ${\cal B}\in {\Bbb R}^{m\times p},$ $m\ge p \ge k,$ we seek the minimum of $\|{\cal A}-{\cal B}Q\|^2$ for all matrices $Q\in {\Bbb R}^{p\times k},$ $Q^TQ=I_{k\times k}$. We introduce a class of relaxation methods for generating sequences of approximations to a minimizer and offer a geometric interpretation of these methods. Results of numerical experiments illustrating the convergence of the methods are given.
Introduction
We begin by defining the set OSt(p; k) of orthogonal Stiefel matrices:
which is a compact submanifold of dimension
of the manifold
O(p) of all p \Theta p orthogonal matrices which has dimension 1
(trace A T
2 denote the standard Frobenius norm in R m\Thetak .
The Procrustes problem for orthogonal Stiefel matrices is to minimize
for all Q 2 OSt(p; k).
Problem (1.2) can be simplified by performing the singular value decomposition
of the matrix B 2 R m\Thetap . Let
where ~
Due to the fact that the last rows of
are zeros we will simplify (1.3) by introducing new notations. We denote
assuming from now on that all oe
A 2 R p\Thetak to be a matrix composed of the first p rows of ~
A. Consequently the
Procrustes minimization on the set of orthogonal Stiefel matrices is :
For given A 2 R p\Thetak and diagonal \Sigma 2 R p\Thetap minimize
for all Q 2 OSt(p; k).
A. W. BOJANCZYK AND A. LUTOBORSKI
The original formulations of the Procrustes problem can be found in [1], [2]. We
may write (1.4) explicitly as
The Procrustes problem has been solved analytically in the orthogonal case when
see [8]. In this case Q 2 O(p) and we have
Provided that the singular value decomposition of \SigmaA is the minimizer
in (1.6) is then
The functional P[A; \Sigma] in (1.5) is a sum of two functionals in Q : the bilinear
functional trace(Q T \Sigma 2 Q) and the linear functional \Gamma2trace(Q T \SigmaA). It is well
known how to minimize each of the functionals separately.
The minimum value of the bilinear functional is equal to the sum of squares of
the k smallest diagonal entries of \Sigma. This result is due to Ky Fan, [9]. The linear
functional is minimized when trace(Q T \SigmaA) is maximized. The maximum of this
trace is given by the sum of the singular values of the matrix \Sigma T A. This upper
bound on the trace functional has been established by J.Von Neumann in [13], see
also [8].
Separate minimization of the quadratic and the linear part are well understood
since we know both the analytical solutions and robust numerical methods. The
analytical solution of the orthogonal Procrustes problem for Stiefel matrices is not
known to the best of our knowledge and constitutes a major challenge.
It will be useful to interpret the minimization (1.4) geometrically. To do that we
define an eccentric Stiefel manifold OSt[\Sigma](p;
The eccentric Stiefel manifold OSt[\Sigma](p; k) is an image of the orthogonal Stiefel
manifold OSt(p; under the linear mapping Q \Gamma! \SigmaQ. The image of a sphere
kg of radius
k in R p\Thetak under this mapping is an ellipsoid in
R p\Thetak of which OSt[\Sigma](p; k) is a subset. The eccentric Stiefel manifold is a compact
set contained in a larger ball in R p\Thetak centered at 0 and of radius
We note that OSt[\Sigma](p; 1) is a standard ellipsoid in R
and OSt[I](p;
Therefore if
min
then a point \SigmaQ   is the projection of A onto the eccentric Stiefel manifold OSt[\Sigma](p; k).
Due to the compactness of the manifold a projection \SigmaQ   exists. The big difficulty
which we face in the task of computing the minimizer Q   is the fact that the
manifold OSt[\Sigma](p; k) is not a convex set.
PROCRUSTES PROBLEM FOR STIEFEL MATRICES 3
2. Notations
Elementary plane rotation by an angle OE is represented by
sin OE cos OE
Elementary plane reflection about the line with slope tan OE is
' cos OE sin OE
sin
For we introduce the following submatrices of
consists of the m-th row of Q
consists of the m and n-th rows of Q
consists of the rows complementary to Q [m;n]
consists of the entries on the intersections of the
m-th and n-th rows and columns of Q:
A plane rotation by an angle OE in the (k; l)-plane in R p is represented by a matrix
G k;l (OE) such that
G k;l (OE) [k;l]
A plane reflection R k;l (OE) in the (k; l)-plane is defined similarily by means of R(OE).
J k;l (p) is the set of all plane rotations and reflections in the (k; l)-plane. J (p) is
the set of plane rotations and reflections in all planes. Clearly
J k;l (p) ae J (p) ae O(p)
3. Relaxation methods for the Procrustes problem
The Stiefel manifold OSt(p; k) is the admissible set for the minimizer of the
functional P in (1.4). This manifold however is not a vector space which poses
severe restrictions on how the succesive approximations can be obtained from the
previous ones. Additive corrections are not admissible, but the Stiefel manifold
is closed with respect to left multiplication by an orthogonal matrix R 2 O(p).
Thus RQ, where Q 2 OSt(p; k), is an admissible approximation. Consequently, we
restrict our considerations to a class of minimization methods which construct the
approximations -
Q to the minimizer Q   by the rule
where Q and -
respectively the current and the next approximations to the
minimizer.
In what follows we will consider only relaxation minimization methods which
seek for the minimizer of the functional P, according to (3.1) with
and M is the dimension of the manifold OSt(p; k). Each R i 2 O(p),
depends on a single parameter whose value results from a scalar
minimization problem. We will refer to the left multiplication by R in (3.1) as to
a sweep. Our relaxation method consists of repeated applications of sweeps which
produce a minimizing sequence for the problem (1.4).
4 A. W. BOJANCZYK AND A. LUTOBORSKI
We will choose matrices R i to be orthogonally similar to a plane rotation or
reflection. Different choices of similarities will lead to different relaxation methods.
We set
and define
may depend on the current approximation to the
. It is the choice of P i that fully determines the relaxation method
(3.1),(3.3),(3.4). The selection of the parameter ff in (3.4) will result from the scalar
minimization
ff
The matrix R i can be viewed as a plane rotation or reflection in a plane spanned
by a pair of columns of the matrix P i . The indices (r; s) of this pair of columns
are selected according to an ordering N of a set of pairs D. The ordering
is an bijection, where D oe
This inclusion guarantees that D contains at least M distinct pairs necessary to
construct an arbitrary Q 2 OSt(p; k) as a product of matrices R i .
It is clear that relaxation methods satisfying (3.5) will always produce a nonincreasing
sequence of the values P[A; \Sigma](Q i ).
If I p\Thetap in (3.4), then R and the sweep (3.1) has the following
particularly simple form
The relaxation method defined by (3.6) will be refered to as a left-sided relaxation
method or LSRM.
If
i is the orthogonal complement of Q i , then
and hence
I k\Thetak
Thus by induction the sweep (3.1) has the form
I k\Thetak
The relaxation method defined by (3.7) will be refered to as a right-sided relaxation
method or RSRM.
A specific type of a right-sided relaxation method was investigated by H. Park
in [11]. The method in [11] is based on the concepts discussed earlier by Ten Berge
and Knol [2] where the Procrustes problem for orthogonal Stiefel matrices, called
there unbalanced problem, is solved by means of iteratively solving a sequence of orthogonal
Procrustes problems, called balanced problems. The relaxation approach
new left-sided relaxation method which is the topic of our paper.
PROCRUSTES PROBLEM FOR STIEFEL MATRICES 5
Another interesting aspect of this approach is a clear geometric interpretation of
the relaxation step. For the study of other minimization methods on submanifolds
of spaces of matrices see [10] and [12] .
4. Planar Procrustes problem
We will now present the left-sided relaxation method. Without loss of generality
let us assume that the planes (r; s) in which transformations operate are chosen in
the row cyclic order, in the way analogous to that used in the cyclic Jacobi method
for the SVD computation, see [5].
In this case N : D
is given by N
Q in (3.6) has the following form
r
Y
where J r;s 2 J r;s (p). Let Q (r;s) be the current approximation in the sweep.
The next approximation to the minimizer is Q . The selection of the
parameter ff results from the scalar minimization
Our main goal now is to show how to find ff in (4.2).
Consider the functional J \Gamma! (4.2), where for simplicity of
notation we omitted all indices. Without loss of generality we assume that N
diag
where G(ff) is a plane rotation (the case of reflection is similar and can be treated
in a completely analogous way). The minimization in (4.2) is precisely the minimization
of
' a 11 a 12 \Delta \Delta \Delta a 1k
a 21 a 22 \Delta \Delta \Delta a 2k
G(ff)
Let UQ \GammaV T
Q be the SVD decomposition of Q [1;2] such that
. Note that the last columns of the matrix B in (4.4)
are always approximated by zero columns, and thus the minimization of f(ff) is
equivalent to minimization restricted to the first two columns, that is
'' c \Gammas
6 A. W. BOJANCZYK AND A. LUTOBORSKI
is the angle of the plane rotation ( or
reflection ) UQ .
We may write (4.5) explicitly as
We now denote
By completing the squares we may represent
F (OE) in the following form:
sin OE
for . Thus the minimization of the functional (4.9) is equivalent
to the following problem, For given C 2 R 2\Theta1 and diagonal Z 2 R 2\Theta2 minimize
for all q 2 OSt(2; 1),
The minimization problem of the type (4.10) will be called a planar Procrustes
problem. Such problem has to be solved on each step of our relaxation method and
is geometrically equivalent to projecting C onto an ellipse. In the next section we
consider two different iterative methods for finding the projection.
5. Projection on an ellipse
The geometrical formulation of the planar Procrustes problem (4.10) is very
simple. Given a point C and an ellipse
ae
oe
in R 2 we want to find a point which is a
projection of C onto E .
This can be achieved in a variety of ways. We describe the classical projection of
a point onto an ellipse due to Appolonius leading to a scalar fourth order algebraic
equation and a method of iterated reflections based on the reflection property of
the ellipse.
PROCRUSTES PROBLEM FOR STIEFEL MATRICES 7
5.1. The hyperbola of Appolonius. Recall the construction of a normal to an
ellipse from a point, see [14], due to Appolonius. With the given ellipse E in (5.1)
and with the point C we associate the hyperbola H in the following way. The
equation of the normal to the ellipse at
z
z
Clearly CS is the normal to E iff the point
Equivalently CS is the normal to E iff the point S satisfies
z
for . (5.3) is a quadratic, in x 1 and x 2 , equation of the
hyperbola H which can be written as
is the center of H with coordinates
F

Figure

1. Hyperbola of Appolonius
The hyperbola H has asymptotes parallel to the axes of the ellipse and passes
through the origin and the point C. H degenerates to a line when the point C is
on one of the axes of the ellipse. This nongeneric, trivial case will not be analyzed
here. H can also be characterized as a locus of the centers of the conics in the
pencil generated by the ellipse and an arbitrary circle centered at C .
To find the coordinates of the projection point S we have to intersect the hyperbola
of Appolonius with the ellipse that is to solve a system of two quadratic
equations (5.1) and (5.4). Using (5.4) to eliminate x 2 from (5.1) we obtain a fourth
8 A. W. BOJANCZYK AND A. LUTOBORSKI
order polynomial equation in x 1
For any specific numerical values of the coefficients this equation can be easily
solved symbolically. A simpler, purely numerical alternative is to solve the system
using Newton's method.
Another alternative is to reduce the system to a scalar equation. Assume that
is in the first quadrant and that
S be the projection of
setting sin OE) in (5.3) and next substituting
leads to the equation t, where
It is easy to see that
and since then the function g(t) is convex and
has one positive root. It can also be seen that for
we have g(t method starting from the initial approximation
will generate a decreasing, convergent sequence of approximations to the root of
5.2. Iterated reflections. Assume, as before, that C is in the first quadrant and
that
Fig.2. Every other case reduces to this one through reflections.
Let
PROCRUSTES PROBLEM FOR STIEFEL MATRICES 9
F
R

Figure

2. Iterated reflections
The reflection property of E says that (for C both inside and outside E) S is
characterized by :
Let , as on Figure 2, be given by
sin
where for C outside
tan OE
cos
(In the case C inside E , the coordinates of the points L and R can be computed
similarly).
Analogously as in the bisection method, we compute
(an intermediate point between L and R) from L and
R by setting OE
(where h \Delta; \Delta i denotes the inner product), then OE R - OE - OE M and we set L
so that OE hold we set L
We thus construct a sequence fLng ae E such that lim n!1
5.3. Some remarks on the general and planar Procrustes problems. The
planar Procrustes problem has several features which the general problem (1.4) of
projection onto the eccentric Stiefel manifold does not posses.
ffl E has the reflection properties and Appolonius normal.
The reflection properties of the ellipse do not extend to the eccentric Stiefel
manifold and in particular not even to ellipsoids in R p . The construction of an
Appolonius normal to the ellipse based on the orthogonality of the ellipse and an
A. W. BOJANCZYK AND A. LUTOBORSKI
associated hyperbola which results in a scalar equation (5.6) is also particular to
the planar problem. As a result in the case p ? k ? 1 our relaxation step, which
amounts to solving a planar Procrustes problem, cannot be directly generalized to
a higher dimensional problem.
ffl A point not belonging to E has either a unique projection onto E or a finite
number of projections.
Hence if
conv(E) then there exits a unique projection S of C onto E characterized
by
for all F 2 conv(E). A point C inside the ellipse E has a non-unique projection
C is on the major axis between points B and \GammaB. Various locations
of C and its projection(s) S are shown on Fig. 3. along with the upper part of
the evolute of the ellipse which is the locus of the centers of curvature of E . The
number of normals 2,3 or 4 that can be drawn from C depends on the position of
C relative to the evolute.
F

Figure

3. Points in the first quadrant and their projections on the ellipse.
The non-uniqueness of S for C on the segment [\GammaB; B] is reflected in non-
differentiability of the function shown on Fig. 4. In general
solutions to the Procrustes problem may form a submanifold of the Stiefel manifold.
Finally we observe that E cuts the plane into two components. However if
PROCRUSTES PROBLEM FOR STIEFEL MATRICES 11
then C is on a sphere of radius
2 in R 3\Theta2 containing the Stiefel manifold OSt[I](3; 2)
but the point kC for any k 2 Rcannot be connected to 0 with a segment intersecting
OSt(3; 2).

Figure

4. Graph of
Both observations concerning the analytical problem are reflected in the computations
and have computational implications.
6. Geometric Interpretation of Left and Right Relaxation Methods
Since the notion of the standard ellipsoid OSt[\Sigma](p; 1) in R p is very intuitive we
will now interpret the minimization problem (1.11) treating matrices in R p\Thetak as k-tuples
of vecors in R p . Let be a given k-tuple of vectors in R p .
be the current approximation to the minimizer. Clearly
the points \Sigmaq i all belong to the ellipsoid OSt[\Sigma](p; 1). Thus the minimization of
P[A; \Sigma](Q) can be interpreted as finding points \Sigmaq
i on the ellipsoid, where q
are
orthonormal vectors, that best match, as measured by P[A; \Sigma](Q), the given vectors
a i in R p .
The relaxation method described in Section 3 can be interpreted as follows. Pick
an orthonormal basis in R p . In the next sweep rotate the current set of vectors q i
as a frame, in planes spanned by all pairs of the vectors from the current basis.
In the left-sided relaxation method the basis is the cannonical basis and is the
same for all sweeps. All relaxation steps are exactly the same, and all amount to
solving a planar Procrustes problem.
In the right-sided relaxation method the basis consists of two subsets and changes
from sweep to sweep. The first subset of the basis consists of the columns of the
current approximation Q and the second subset consists of the columns of the
orthogonal complement Q ? of Q.
Working only with the columns of Q is equivalent to the so-called balanced
Procrustes problem studied by Park in [11] which can be solved by means of an
A. W. BOJANCZYK AND A. LUTOBORSKI
SVD computation. The relaxation step in [11] for the balanced problem consists of
computing the SVD of the 2 \Theta 2 matrix (\Sigmaq r In our relaxation
setting, the relaxation step in the right-sided relaxation method requires solving
the scalar minimization problem (3.5),
min
which leads to a linear equation in tangent of ff and is equivalent to the 2 \Theta 2 SVD
computation in [11]. Each of these 2 \Theta 2 steps is a rotation of vectors \Sigmaq r and
\Sigmaq s in the plane spanned by q r and q s so the rotated vectors on the ellipsoid best
approximate the two given vectors a r and a s .
However, as the columns of Q do not span the whole space R p , it might happen
that
k g, and hence it might not be possible
to generate a sequence of approximations that will converge to Q   . In order to
overcome this problem the matrix Q is extended by its orthogonal complement
subproblems in [11] involving vectors from both subsets are referred to
in [11] as unbalanced subproblems. These scalar minimizations have the following
min
c
s
That is, the unbalanced subproblem is to find a vector on the ellipsoid in the plane
spanned by q r and q s closest to the given vector a r . As the intersection of this plane
and the ellipsoid is an ellipse, the unbalanced subproblem can be expressed as a
planar Procrustes problem (4.10) and any of the algorithms discussed in Section 5
can be used to solve this unbalanced problem.
Other choices of bases may be possible but the choices leading to the left and
the right sided-relaxation methods seem to be the most natural.
7. Numerical experiments
In this section we present numerical experiments illustrating the behavior of
the left and the right relaxation methods discussed in Section 3. We will start by
summarizing the left and the right relaxation methods given below in pseudocodes.
Given A 2 R p\Thetak ,
construct of sequences of Stiefel matrices approximating the minimizer of (1.4).
Algorithm LSRM:
1. Initialization:
set Maxstep,
2. Iterate sweeps:
while (r threshold and n ! Maxstep
solve planar Procrustes problem
PROCRUSTES PROBLEM FOR STIEFEL MATRICES 13Algorithm RSRM:
1. Initialization:
set Maxstep,
2. Iterate sweeps:
while (r threshold and n ! Maxstep
solve min ff
solve planar Procrustes problem
measure the cost of the two methods by the number of sweeps performed by
each of the two methods.
A sweep in the LSRM method consists of p(p+1)=2 planar Procrustes problems.
Each planar Procrustes problem requires computation of the SVD of a 2 \Theta k matrix.
This can be achieved by first computing the QR decomposition followed by a 2 \Theta 2
SVD problem. After the SVD is calculated, a projection on an ellipse has to be
determined. The cost of a sweep is approximately O(kp 2 ) floating point operations.
A sweep in the RSRM method consists of k(k computations of p \Theta 2 SVD
problems. In adition, there are k(p \Gamma Procrustes problems, each requiring
computation of the SVD of a p \Theta 2 matrix followed by computation of a projection
on an ellipse. Thus the cost of a sweep is again approximately O(kp 2 ) floating point
operations.
Surely, the precise cost of a sweep will depend on the number of iterations needed
for obtaining satisfactory projections on the resulting ellipses. For each projection,
this will depend on the location of the point being projected as well as the shape
of the ellipse. Computation of the projection will be most costly when the ellipse
is flat.
As can be seen, sweeps in the two methods may have different costs. However,
the number of sweeps performend by each of the methods will give some bases for
comparing the convergence behavior of the two methods.
We begin by illustrating the behavior of the LSRM method for finding Q   in
the Procrustes problem with
The initial approximation is I 4\Theta2 . Some intermediate values of Q are listed
in

Table

1.
14 A. W. BOJANCZYK AND A. LUTOBORSKI

Table

1. Matrices Q in a minimizing sequence generated by LSRM.
We will now present comperative numerical results for the LSRM and RSMR
methods.
Recall that the functional P is a sum of a linear and a bilinear term. We will
consider classes of examples when the functional can be approximated by its linear
or the bilinear term.
In the first class of examples the linear term dominates the bilinear term or in
other words when jjAjj ?? jj\Sigmajj. We deal here with a perturbed linear functional.
The minimum of the functional P can be approximated by the sum of singular
values of \Sigma T A.
The second class of examples consists of cases when the quadratic term dominates
the linear term, that is when jjAjj !! jj\Sigmajj. We deal here with a perturbed bilinear
functional. Then the minimum value of the functional P can be approximated by
the sum of the k smallest singular values of \Sigma.
The third class of exampless consists of cases when the functional is quadratic,
that is when A - \SigmaQ for some Q 2 OSt(p; k). The minimum of the functional is
then close to zero.
PROCRUSTES PROBLEM FOR STIEFEL MATRICES 15
In each class of examples we pick two different matrices \Sigma: one corresponding
to the ellipsoid being almost a sphere, that is when \Sigma - I, the other corresponding
to the ellipsoid being very flat in one or more planes, that is when oe 1
is large.
The algorithms were written in MATLAB 4.2 and run on an HP9000 workstation
with the machine relative precision
As the initial approximation we took I p\Thetak for the LSRM, and
I p\Thetap for the RSRM. The planar Procrustes solver used was based on
the hyperbola of Appolonius (the iterated reflections solver was giving numerically
equivalent results). Some representative results are shown in Tables 2-6.
RSRM LSRM
6

Table

2.
2.
RSRM LSRM
6

Table

3.

Table

2 illustrates the behavior of the two methods when the ellipsoid is almost
a sphere and when there exists Q such that A. That is the bilinear and
the linear terms are of comparable size. The experiments suggest that the LSRM
requires less sweeps to obtain a satisfactory approximation to the minimizer.

Table

3 illustrates the behavior of the two methods when the length of the half
of the ellipsoid's axes is approximately 1.0 and the other half is approximately 0.01.
A. W. BOJANCZYK AND A. LUTOBORSKI
In addition, there exists Q such that A. In this case the convergance of the
RSRM is particularly slow. We observed that, at least initially, the RSRM fails to
locate the minimizer in OSt(4; 2) being unable to establish the proper signs of the
entries of the matrix Q   . The LSRM on the other hand approximates the minimizer
correctly.
RSRM LSRM
6

Table

4.
2.
RSRM LSRM
28 -1.43e-02 3.16e+01 0.00e+00

Table

5.
2.

Table

4 illustrates the behavior of the two methods when the ellipsoid is almost a
sphere but now A is chosen so That is the bilinear term dominates
the linear term. In this case the minimum of the functional can be estimated by
the minimum value of the bilinear term. In Table 4 esterror denotes the difference
between the minimum value of the bilinear term and the computed value of the
functional, and
Qk where Q and -
Q are the
last and penultimate approximations to the minimizer. The experiments suggest
that the LSRM requires less sweeps to obtain a satisfactory approximation to the
minimizer.
PROCRUSTES PROBLEM FOR STIEFEL MATRICES 17

Table

5 illustrates the behavior of the two methods when the ellipsoid is almost a
sphere but now A is chosen so That is the linear term dominates the
quadratic terms. In this case the minimum of the functional can be estimated by
the minimum value of the linear term. In Table 5 esterror denotes the difference
between the minimum value of the linear term and the computed value of the
functional. The experiments suggest that the RSRM requires less sweeps to obtain
a satisfactory approximation to the minimizer.
8. Remarks on Constrained Linear Least Squares Problems
The quadratically constrained linear least squares problem
min
arises in many applications [6], [7], [3]. By changing variables this problem can be
transformed into a special Procrustes problem. The Procrustes problem is
min
Procrustes problem is
equivalent to projecting the point a onto the ellipsoid OSt[\Sigma](p; 1).
1. It is clear that the vector has the direction of the normal
vector to the ellipsoid at the point \Sigmaq. Thus if \Sigmaq   is the projection of a on the
ellipsoid then the vector a \Gamma \Sigmaq   is parallel to the vector n. Thus there exists a
scalar fi so
can obtain an equation for fi
The parameter fi can be computed by solving this equation. Then the components
of q are given by
The equation (8.4) is the so-called secular equation which characterizes the critical
points of the Lagrangian
see [8]. Thus the multiplier fi in (8.3) is the Lagrange multiplier - in (8.5).



--R



Algorithms for the regularization of ill-conditioned least squares problems
On the stationary values of a second-degree polynomial on the unit sphere
The cyclic Jacobi method for computing the principal values of a complex matrix
Least Squares with a Quadratic Constraint
Quadratically constrained least squares and quadratic problems
Matrix Computations Johns Hopkins
On a theorem of Weyl concerning eigenvalues of linear transformations
On the convergence of the Euler-Jacobi method
A parallel algorithm for the unbalanced orthogonal Procrustes problem
Optimization techniques on Riemannian manifolds
Some matrix inequalities and metrization of the matrix space
Enzyklopadie der Elementar Mathematik B.
--TR
stiefel manifolds;projections on ellipsoids;relaxation methods;procrustes problem