Constructing Packings in Grassmannian Manifolds
via Alternating Projection
This paper describes a numerical method for finding good packings in Grassmannian manifolds equipped with various metrics. This investigation also encompasses packing in projective spaces. In each case, producing a good packing is equivalent to constructing a matrix that has certain structural and spectral properties. By alternately enforcing the structural condition and then the spectral condition, it is often possible to reach a matrix that satisfies both. One may then extract a packing from this matrix.
This approach is both powerful and versatile. In cases where experiments have been performed, the alternating projection method yields packings that compete with the best packings recorded. It also extends to problems that have not been studied numerically. For example, it can be used to produce packings of subspaces in real and complex Grassmannian spaces equipped with the Fubini–Study distance; these packings are valuable in wireless communications. One can prove that some of the novel configurations constructed by the algorithm have packing diameters that are nearly optimal.
Key words and phrases:Combinatorial optimization, packing, projective spaces, Grassmannian spaces, Tammes’ Problem
2000 Mathematics Subject Classification:Primary: 51N15, 52C17
Let us begin with the standard facetious example. Imagine that several mutually inimical nations build their capital cities on the surface of a featureless globe. Being concerned about missile strikes, they wish to locate the closest pair of cities as far apart as possible. In other words, what is the best way to pack points on the surface of a two-dimensional sphere?
This question, first discussed by the Dutch biologist Tammes [Tam30], is the prototypical example of packing in a compact metric space. It has been studied in detail for the last 75 years. More recently, researchers have started to ask about packings in other compact spaces. In particular, several communities have investigated how to arrange subspaces in a Euclidean space so that they are as distinct as possible. An equivalent formulation is to find the best packings of points in a Grassmannian manifold. This problem has applications in quantum computing and wireless communications. There has been theoretical interest in subspace packing since the 1960s [Tót65], but the first detailed numerical study appears in a 1996 paper of Conway, Hardin, and Sloane [CHS96].
The aim of this paper is to describe a flexible numerical method that can be used to construct packings in Grassmannian manifolds equipped with several different metrics. The rest of this introduction provides a formal statement of abstract packing problems, and it offers an overview of our approach to solving them.
1.1. Abstract Packing Problems
Although we will be working with Grassmannian manifolds, it is more instructive to introduce packing problems in an abstract setting. Let be a compact metric space endowed with the distance function . The packing diameter of a finite subset is the minimum distance between some pair of distinct points drawn from . That is,
In other words, the packing diameter of a set is the diameter of the largest open ball that can be centered at each point of the set without encompassing any other point. (It is also common to study the packing radius, which is half the diameter of this ball.) An optimal packing of points is an ensemble that solves the mathematical program
where returns the cardinality of a finite set. The optimal packing problem is guaranteed to have a solution because the metric space is compact and the objective is a continuous function of the ensemble .
This article focuses on a feasibility problem closely connected with optimal packing. Given a number , the goal is to produce a set of points for which
This problem is notoriously difficult to solve because it is highly nonconvex, and it is even more difficult to determine the maximum value of for which the feasibility problem is soluble. This maximum value of corresponds with the diameter of an optimal packing.
1.2. Alternating Projection
We will attempt to solve the feasibility problem (1.1) in Grassmannian manifolds equipped with a number of different metrics, but the same basic algorithm applies in each case. Here is a high-level description of our approach.
First, we show that each configuration of subspaces is associated with a block Gram matrix whose blocks control the distances between pairs of subspaces. Then we prove that a configuration solves the feasibility problem (1.1) if and only if its Gram matrix possesses both a structural property and a spectral property. The overall algorithm consists of the following steps.
Choose an initial configuration and construct its matrix.
Alternately enforce the structural condition and the spectral condition in hope of reaching a matrix that satisfies both.
Extract a configuration of subspaces from the output matrix.
In our work, we choose the initial configuration randomly and then remove similar subspaces from it with a simple algorithm. One can imagine more sophisticated approaches to constructing the initial configuration.
Flexibility and ease of implementation are the major advantages of alternating projection. This article demonstrates that appropriate modifications of this basic technique allow us to construct solutions to the feasibility problem in Grassmannian manifolds equipped with various metrics. Some of these problems have never been studied numerically, and the experiments point toward intriguing phenomena that deserve theoretical attention. Moreover, we believe that the possibilities of this method have not been exhausted and that it will see other applications in the future.
Alternating projection does have several drawbacks. It may converge very slowly, and it does not always yield a high level of numerical precision. In addition, it may not deliver good packings when the ambient dimension or the number of subspaces in the configuration is large.
1.3. Motivation and Related Work
This work was motivated by applications in electrical engineering. In particular, subspace packings solve certain extremal problems that arise in multiple-antenna communication systems [ZT02, HMR00, LJSH04]. This application requires complex Grassmannian packings that consist of a small number of subspaces in an ambient space of low dimension. Our algorithm is quite effective in this parameter regime. The resulting packings fill a significant gap in the literature, since existing tables consider only the real case [Slo04a]. See Section 6.1 for additional discussion of the wireless application.
The approach to packing via alternating projection was discussed in a previous publication [TDJS05], but the experiments were limited to a single case. We are aware of several other numerical methods that can be used to construct packings in Grassmannian manifolds [CHS96, Tro01, ARU01]. These techniques rely on ideas from nonlinear programming.
1.4. Historical Interlude
The problem of constructing optimal packings in various metric spaces has a long and lovely history. The most famous example may be Kepler’s Conjecture that an optimal packing of spheres in three-dimensional Euclidean space111The infinite extent of a Euclidean space necessitates a more subtle definition of an optimal packing. locates them at the points of a face-centered cubic lattice. For millennia, greengrocers have applied this theorem when stacking oranges, but it has only been established rigorously within the last few years [Hal04]. Packing problems play a major role in modern communications because error-correcting codes may be interpreted as packings in the Hamming space of binary strings [CT91]. The standard reference on packing is the magnum opus of Conway and Sloane [CS98]. Classical monographs on the subject were written by L. Fejes Tóth [Tót64] and C. A. Rogers [Rog64].
The idea of applying alternating projection to feasibility problems first appeared in the work of von Neumann [vN50]. He proved that an alternating projection between two closed subspaces of a Hilbert space converges to the orthogonal projection of the initial iterate onto the intersection of the two subspaces. Cheney and Goldstein subsequently showed that an alternating projection between two closed, convex subsets of a Hilbert space always converges to a point in their intersection (provided that the intersection is nonempty) [CG59]. This result does not apply in our setting because one of the constraint sets we define is not convex.
1.5. Outline of Article
Here is a brief overview of this article. In Section 2, we develop a basic description of Grassmannian manifolds and present some natural metrics. Section 3 explains why alternating projection is a natural algorithm for producing Grassmannian packings, and it outlines how to apply this algorithm for one specific metric. Section 4 gives some theoretical upper bounds on the optimal diameter of packings in Grassmannian manifolds. Section 5 describes the outcomes of an extensive set of numerical experiments and explains how to apply the algorithm to other metrics. Section 6 offers some discussion and conclusions. Appendix A explores how our methodology applies to Tammes’ Problem of packing on the surface of a sphere. Finally, Appendix B contains tables and figures that detail the experimental results.
2. Packing in Grassmannian Manifolds
This section introduces our notation and a simple description of the Grassmannian manifold. It presents several natural metrics on the manifold, and it shows how to represent a configuration of subspaces in matrix form.
We work in the vector space . The symbol denotes the complex-conjugate transpose of a vector (or matrix). We equip the vector space with its usual inner product . This inner product generates the norm via the formula .
The -dimensional identity matrix is ; we sometimes omit the subscript if it is unnecessary. A square matrix is positive semidefinite when its eigenvalues are all nonnegative. We write to indicate that is positive semidefinite.
A square, complex matrix is unitary if it satisfies . If in addition the entries of are real, the matrix is orthogonal. The unitary group can be presented as the collection of all unitary matrices with ordinary matrix multiplication. The real orthogonal group can be presented as the collection of all real orthogonal matrices with the usual matrix multiplication.
Suppose that is a general matrix. The Frobenius norm is calculated as , where the trace operator sums the diagonal entries of the matrix. The spectral norm is denoted by ; it returns the largest singular value of . Both these norms are unitarily invariant, which means that whenever and are unitary.
2.2. Grassmannian Manifolds
The (complex) Grassmannian manifold is the collection of all -dimensional subspaces of . This space is isomorphic to a quotient of unitary groups:
To understand the equivalence, note that each orthonormal basis from can be split into vectors, which span a -dimensional subspace, and vectors, which span the orthogonal complement of that subspace. To obtain a unique representation for the subspace, it is necessary to divide by isometries that fix the subspace and by isometries that fix its complement. It is evident that is always isomorphic to .
Similarly, the real Grassmannian manifold is the collection of all -dimensional subspaces of . This space is isomorphic to a quotient of orthogonal groups:
If we need to refer to the real and complex Grassmannians simultaneously, we write .
In the theoretical development, we concentrate on complex Grassmannians since the development for the real case is identical, except that all the matrices are real-valued instead of complex-valued. A second reason for focusing on the complex case is that complex packings arise naturally in wireless communications [LJS03].
When each subspace has dimension , the Grassmannian manifold reduces to a simpler object called a projective space. The elements of a projective space can be viewed as lines through the origin of a Euclidean space. The standard notation is . We will spend a significant amount of attention on packings of this manifold.
2.3. Principal Angles
Suppose that and are two subspaces in . These subspaces are inclined against each other by different principal angles. The smallest principal angle is the minimum angle formed by a pair of unit vectors drawn from . That is,
The second principal angle is defined as the smallest angle attained by a pair of unit vectors that is orthogonal to the first pair, i.e.,
The remaining principal angles are defined analogously. The sequence of principal angles is nondecreasing, and it is contained in the range . We only consider metrics that are functions of the principal angles between two subspaces.
Let us present a more computational definition of the principal angles [BG73]. Suppose that the columns of and form orthonormal bases for the subspaces and . More rigorously, is a matrix that satisfies and . The matrix has an analogous definition. Next we compute a singular value decomposition of the product :
where and are unitary matrices and is a nonnegative, diagonal matrix with nonincreasing entries. The matrix of singular values is uniquely determined, and its entries are the cosines of the principal angles between and :
This definition of the principal angles is most convenient numerically because singular value decompositions can be computed efficiently with standard software. We also note that this definition of the principal angles does not depend on the choice of matrices and that represent the two subspaces.
2.4. Metrics on Grassmannian Manifolds
Grassmannian manifolds admit many interesting metrics, which lead to different packing problems. This section describes some of these metrics.
The chordal distance between two -dimensional subspaces and is given by
The values of this metric range between zero and . The chordal distance is the easiest to work with, and it also yields the most symmetric packings [CHS96].
The spectral distance is
The values of this metric range between zero and one. As we will see, this metric promotes a special type of packing called an equi-isoclinic configuration of subspaces.
The geodesic distance is
This metric takes values between zero and . From the point of view of differential geometry, the geodesic distance is very natural, but it does not seem to lead to very interesting packings [CHS96], so we will not discuss it any further.
Grassmannian manifolds support several other interesting metrics, some of which are listed in [BN02]. In case we are working in a projective space, i.e., , all of these metrics reduce to the acute angle between two lines or the sine thereof. Therefore, the metrics are equivalent up to a monotonically increasing transformation, and they promote identical packings.
2.5. Representing Configurations of Subspaces
Suppose that is a collection of subspaces in . Let us develop a method for representing this configuration numerically. To each subspace , we associate a (nonunique) matrix whose columns form an orthonormal basis for that subspace, i.e., and . Now collate these matrices into a configuration matrix
In the sequel, we do not distinguish between the configuration and the matrix .
The Gram matrix of is defined as the matrix . By construction, the Gram matrix is positive semidefinite, and its rank does not exceed . It is best to regard the Gram matrix as an block matrix comprised of blocks, and we index it as such. Observe that each block satisfies
In particular, each diagonal block is an identity matrix. Meanwhile, the singular values of the off-diagonal block equal the cosines of the principal angles between the two subspaces and .
Conversely, let be an block matrix with each block of size . Suppose that the matrix is positive semidefinite, that its rank does not exceed , and that its diagonal blocks are identity matrices. Then we can factor where is a configuration matrix. That is, the columns of form orthogonal bases for different -dimensional subspaces of .
As we will see, each metric on the Grassmannian manifold leads to a measure of “magnitude” for the off-diagonal blocks on the Gram matrix . A configuration solves the feasibility problem (1.1) if and only if each off-diagonal block of its Gram matrix has sufficiently small magnitude. So solving the feasibility problem is equivalent to producing a Gram matrix with appropriate properties.
3. Alternating Projection for Chordal Distance
In this section, we elaborate on the idea that solving the feasibility problem is equivalent with constructing a Gram matrix that meets certain conditions. These conditions fall into two different categories: structural properties and spectral properties. This observation leads naturally to an alternating projection algorithm for solving the feasibility problem. The algorithm alternately enforces the structural properties and then the spectral properties in hope of producing a Gram matrix that satisfies them all. This section illustrates how this approach unfolds when distances are measured with respect to the chordal metric. In Section 5, we describe adaptations for other metrics.
3.1. Packings with Chordal Distance
Suppose that we seek a packing of subspaces in equipped with the chordal distance. If is a configuration of subspaces, its packing diameter is
Given a parameter , the feasibility problem elicits a configuration that satisfies
We may rearrange this inequality to obtain a simpler condition:
In fact, we may formulate the feasibility problem purely in terms of the Gram matrix. Suppose that the configuration satisfies (3.1) with parameter . Then its Gram matrix must have the following six properties:
Each diagonal block of is an identity matrix.
for each .
is positive semidefinite.
has rank or less.
has trace .
Some of these properties are redundant, but we have listed them separately for reasons soon to become apparent. Conversely, suppose that a matrix satisfies Properties 1–6. Then it is always possible to factor it to extract a configuration of subspaces that solves (3.1). The factorization of can be obtained most easily from an eigenvalue decomposition of .
3.2. The Algorithm
Observe that Properties 1–3 are structural properties. By this, we mean that they constrain the entries of the Gram matrix directly. Properties 4–6, on the other hand, are spectral properties. That is, they control the eigenvalues of the matrix. It is not easy to enforce structural and spectral properties simultaneously, so we must resort to half measures. Starting from an initial matrix, our algorithm will alternately enforce Properties 1–3 and then Properties 4–6 in hope of reaching a matrix that satisfies all six properties at once.
To be more rigorous, let us define the structural constraint set
Although the structural constraint set evidently depends on the parameter , we will usually eliminate from the notation for simplicity. We also define the spectral constraint set
Both constraint sets are closed and bounded, hence compact. The structural constraint set is convex, but the spectral constraint set is not.
To solve the feasibility problem (3.1), we must find a matrix that lies in the intersection of and . This section states the algorithm, and the succeeding two sections provide some implementation details.
Algorithm 1 (Alternating Projection).
A Hermitian matrix
The maximum number of iterations
A matrix that belongs to and whose diagonal blocks are identity matrices
Determine a matrix that solves
Determine a matrix that solves
If , return to Step 2.
Define the block-diagonal matrix .
Return the matrix
The iterates generated by this algorithm are not guaranteed to converge in norm. Therefore, we have chosen to halt the algorithm after a fixed number of steps instead of checking the behavior of the sequence of iterates. We discuss the convergence properties of the algorithm in the sequel.
The scaling in the last step normalizes the diagonal blocks of the matrix but preserves its inertia (i.e., numbers of negative, zero, and positive eigenvalues). Since is a positive-semidefinite matrix with rank or less, the output matrix shares these traits. It follows that the output matrix always admits a factorization where is a configuration matrix. Property 3 is the only one of the six properties that may be violated.
3.3. The Matrix Nearness Problems
To implement Algorithm 1, we must solve the matrix nearness problems in Steps 2 and 3. The first one is straightforward.
Let be an Hermitian matrix. With respect to the Frobenius norm, the unique matrix in nearest to has diagonal blocks equal to the identity and off-diagonal blocks that satisfy
It is rather more difficult to find a nearest matrix in the spectral constraint set. To state the result, we define the plus operator by the rule .
Let be an Hermitian matrix whose eigenvalue decomposition is with the eigenvalues arranged in nonincreasing order: . With respect to the Frobenius norm, a matrix in closest to is given by
where the scalar is chosen so that
This best approximation is unique provided that .
The nearest matrix described by this theorem can be computed efficiently from an eigenvalue decomposition of . (See [GVL96] for computational details.) The value of is uniquely determined, but one must solve a small rootfinding problem to solve it. The bisection method is an appropriate technique since the plus operator is nondifferentiable. We omit the details, which are routine.
Given an Hermitian matrix , denote by the vector of eigenvalues arranged in nonincreasing order. Then we may decompose for some unitary matrix .
Finding the matrix in closest to is equivalent to solving the optimization problem
First, we fix the eigenvalues of and minimize with respect to the unitary part of its eigenvalue decomposition. In consequence of the Hoffman–Wielandt Theorem [HJ85], the objective function is bounded below:
Equality holds if and only if and are simultaneously diagonalizable by a unitary matrix. Therefore, if we decompose , the objective function attains its minimal value whenever . Note that the matrix may not be uniquely determined.
We find the optimal vector of eigenvalues for the matrix by solving the (strictly) convex program
This minimization is accomplished by an application of Karush–Kuhn–Tucker theory [Roc70]. In short, the top eigenvalues of are translated an equal amount, and those that become negative are set to zero. The size of the translation is chosen to fulfill the third condition (which controls the trace of ). The entries of the optimal are nonincreasing on account of the ordering of .
Finally, the uniqueness claim follows from the fact that the eigenspace associated with the top eigenvectors of is uniquely determined if and only if . ∎
3.4. Choosing an Initial Configuration
The success of the algorithm depends on adequate selection of the input matrix . We have found that the following strategy is reasonably effective. It chooses random subspaces and adds them to the initial configuration only if they are sufficiently distant from the subspaces that have already been chosen.
Algorithm 4 (Initial Configuration).
The ambient dimension , the subspace dimension , and the number of subspaces
An upper bound on the similarity between subspaces
The maximum number of random selections
A matrix from whose off-diagonal blocks also satisfy
Initialize and .
Increment . If , print a failure notice and stop.
Pick a matrix whose range is a uniformly random subspace in .
If for each , then increment .
If , return to Step 2.
Form the matrix .
Return the Gram matrix .
To implement Step 3, we use the method developed in [Ste80]. Draw a matrix whose entries are iid complex, standard normal random variables, and perform a decomposition. The first columns of the unitary part of the decomposition form an orthonormal basis for a random -dimensional subspace.
The purpose of the parameter is to prevent the starting configuration from containing blocks that are nearly identical. The extreme case places no restriction on the similarity between blocks. If is chosen too small (or if we are unlucky in our random choices), then this selection procedure may fail. For this reason, we add an iteration counter to prevent the algorithm from entering an infinite loop. We typically choose values of very close to the maximum value.
3.5. Theoretical Behavior of Algorithm
It is important to be aware that packing problems are typically difficult to solve. Therefore, we cannot expect that our algorithm will necessarily produce a point in the intersection of the constraint sets. One may ask whether we can make any guarantees about the behavior of Algorithm 1. This turns out to be difficult. Indeed, there is potential that an alternating projection algorithm will fail to generate a convergent sequence of iterates [Mey76]. Nevertheless, it can be shown that the sequence of iterates has accumulation points and that these accumulation points satisfy a weak structural property.
In practice, the alternating projection algorithm seems to converge, but a theoretical justification for this observation is lacking. A more serious problem is that the algorithm frequently requires as many as 5000 iterations before the iterates settle down. This is one of the major weaknesses of our approach.
For reference, we offer the best theoretical convergence result that we know. The distance between a matrix and a compact collection of matrices is defined as
It can be shown that the distance function is Lipschitz, hence continuous.
Theorem 5 (Global Convergence).
Suppose that Algorithm 1 generates an infinite sequence of iterates . This sequence has at least one accumulation point.
Every accumulation point lies in .
Every accumulation point satisfies
Every accumulation point satisfies
The existence of an accumulation point follows from the compactness of the constraint sets. The algorithm does not increase the distance between successive iterates, which is bounded below by zero. Therefore, this distance must converge. The distance functions are continuous, so we can take limits to obtain the remaining assertions. ∎
A more detailed treatment requires the machinery of point-to-set maps, and it would not enhance our main discussion. Please see the appendices of [TDJS05] for additional information.
4. Bounds on the Packing diameter
To assay the quality of the packings that we produce, it helps to have some upper bounds on the packing diameter. If a configuration of subspaces has a packing diameter close to the upper bound, that configuration must be a nearly optimal packing. This approach allows us to establish that many of the packings we construct numerically have packing diameters that are essentially optimal.
Theorem 6 (Conway–Hardin–Sloane [Chs96]).
The packing diameter of subspaces in the Grassmannian manifold equipped with chordal distance is bounded above as
If the bound is met, all pairs of subspaces are equidistant. When , the bound is attainable only if . When , the bound is attainable only if .
The complex case is not stated in [CHS96], but it follows from an identical argument. We refer to (4.1) as the Rankin bound for subspace packings with respect to the chordal distance. The reason for the nomenclature is that the result is established by embedding the chordal Grassmannian manifold into a Euclidean sphere and applying the classical Rankin bound for sphere packing [Ran47].
It is also possible to draw a corollary on packing with respect to the spectral distance; this result is novel. A subspace packing is said to be equi-isoclinic if all the principal angles between all pairs of subspaces are identical [LS73].
We have the following bound on the packing diameter of subspaces in the Grassmannian manifold equipped with the spectral distance.
If the bound is met, the packing is equi-isoclinic.
We refer to (4.2) as the Rankin bound for subspace packings with respect to spectral distance.
The power mean inequality (equivalently, Hölder’s inequality) yields
For angles between zero and , equality holds if and only if . It follows that
If the second inequality is met, then all pairs of subspaces are equidistant with respect to the chordal metric. Moreover, if the first inequality is met, then the principal angles between each pair of subspaces are constant. Together, these two conditions imply that the packing is equi-isoclinic. ∎
An upper bound on the maximum number of equi-isoclinic subspaces is available. Its authors do not believe that it is sharp.
Theorem 8 (Lemmens–Seidel [Ls73]).
The maximum number of equi-isoclinic -dimensional subspaces of is no greater than
Similarly, the maximum number of equi-isoclinic -dimensional subspaces of does not exceed
Our approach to packing is experimental rather than theoretical, so the real question is how Algorithm 1 performs in practice. In principle, this question is difficult to resolve because the optimal packing diameter is unknown for almost all combinations of and . Whenever possible, we compared our results with the Rankin bound and with the “world record” packings tabulated by N. J. A. Sloane and his colleagues [Slo04a]. In many cases, the algorithm was able to identify a nearly optimal packing. Moreover, it yields interesting results for packing problems that have not received numerical attention.
In the next subsection, we describe detailed experiments on packing in real and complex projective spaces. Then, we move on to packings of subspaces with respect to the chordal distance. Afterward, we study the spectral distance and the Fubini–Study distance.
5.1. Projective Packings
Line packings are the simplest type of Grassmannian packing, so they offer a natural starting point. Our goal is to produce the best packing of lines in . In the real case, Sloane’s tables allow us to determine how much our packings fall short of the world record. In the complex setting, there is no comparable resource, so we must rely on the Rankin bound to gauge how well the algorithm performs.
Let us begin with packing in real projective spaces. We attempted to construct configurations of real lines whose maximum absolute inner product fell within of the best value tabulated in [Slo04a]. For pairs with and , we computed the putatively optimal value of the feasibility parameter from Sloane’s data and equation (3.2). In each of 10 trials, we constructed a starting matrix using Algorithm 4 with parameters and . (Recall that the value of determines the maximum number of random subspaces that are drawn when trying to construct the initial configuration.) We applied alternating projection, Algorithm 1, with the computed value of and the maximum number of iterations . (Our numerical experience indicates that increasing the maximum number of iterations beyond 5000 does not confer a significant benefit.) We halted the iteration in Step 4 if the iterate exhibited no off-diagonal entry with absolute value greater than . After 10 trials, we recorded the largest packing diameter attained, as well as the average value of the packing diameter. We also recorded the average number of iterations the alternating projection required per trial.
Table B delivers the results of this experiment. Following Sloane, we have reported the degrees of arc subtended by the closest pair of lines. We believe that it is easiest to interpret the results geometrically when they are stated in this fashion. All the tables and figures related to packing are collated at the back of this paper for easy comparison.
According to the table, the best configurations produced by alternating projection consistently attain packing diameters tenths or hundredths of a degree away from the best configurations known. The average configurations returned by alternating projection are slightly worse, but they usually fall within a degree of the putative optimal. Moreover, the algorithm finds certain configurations with ease. For the pair , fewer than 1000 iterations are required on average to achieve a packing within 0.001 degrees of optimal.
A second observation is that the alternating projection algorithm typically performs better when the number of points is small. The largest errors are all clustered at larger values of . A corollary observation is that the average number of iterations per trial tends to increase with the number of points.
There are several anomalies that we would like to point out. The most interesting pathology occurs at the pair . The best packing diameter calculated by alternating projection is about worse than the optimal configuration, and it is also worse than the best packing diameter computed for the pair . From Sloane’s tables, we can see that the (putative) optimal packing of 19 lines in is actually a subset of the best packing of 20 lines. Perhaps the fact that this packing is degenerate makes it difficult to construct. A similar event occurs (less dramatically) at the pair . The table also shows that the algorithm performs less effectively when the number of lines exceeds 20.
In complex projective spaces, this methodology does not apply because there are no tables available. In fact, we only know of one paper that contains numerical work on packing in complex projective spaces [ARU01], but it gives very few examples of good packings. The only method we know for gauging the quality of a complex line packing is to compare it against an upper bound. The Rankin bound for projective packings, which is derived in Section 4, states that every configuration of lines in either or satisfies the inequality
This bound is attainable only for rare combinations of and . In particular, the bound can be met in only if . In the space , attainment requires that . Any arrangement of lines that meets the Rankin bound must be equiangular. These optimal configurations are called equiangular tight frames. See [SJ03, HP04, TDJS05, STDJ07] for more details.
We performed some ad hoc experiments to produce configurations of complex lines with large packing diameters. For each pair , we used the Rankin bound to determine a lower limit on the feasibility parameter . Starting matrices were constructed with Algorithm 4 using values of ranging between 0.9 and 1.0. (Algorithm 4 typically fails for smaller values of .) For values of the feasibility parameter between the minimal value and twice the minimal value, we performed 5000 iterations of Algorithm 1, and we recorded the largest packing diameter attained during these trials.
Table B compares our results against the Rankin bound. We see that many of the complex line configurations have packing diameters much smaller than the Rankin bound, which is not surprising because the bound is usually not attainable. Some of our configurations fall within a thousandth of a degree of the bound, which is essentially optimal.
Table B contains a few oddities. In , the best packing diameter computed for is worse than the packing diameter for . This configuration of 25 lines is an equiangular tight frame, which means that it is an optimal packing [TDJS05, Table 1]. It seems likely that the optimal configurations for the preceding values of are just subsets of the optimal arrangement of 25 lines. As before, it may be difficult to calculate this type of degenerate packing. A similar event occurs less dramatically at the pair and at the pairs and .
Figure 1 compares the quality of the best real projective packings from [Slo04a] with the best complex projective packings that we obtained. It is natural that the complex packings are better than the real packings because the real projective space can be embedded isometrically into the complex projective space. But it is remarkable how badly the real packings compare with the complex packings. The only cases where the real and complex ensembles have the same packing diameter occur when the real configuration meets the Rankin bound.
5.2. The Chordal Distance
Emboldened by this success with projective packings, we move on to packings of subspaces with respect to the chordal distance. Once again, we are able to use Sloane’s tables for guidance in the real case. In the complex case, we fall back on the Rankin bound.
For each triple , we determined a value for the feasibility parameter from the best packing diameter Sloane recorded for subspaces in , along with equation (3.2). We constructed starting points using the modified version of Algorithm 4 with , which represents no constraint. (We found that the alternating projection performed no better with initial configurations generated from smaller values of .) Then we executed Algorithm 1 with the calculated value of for 5000 iterations.
Table B demonstrates how the best packings we obtained compare with Sloane’s best packings. Many of our real configurations attained a squared packing diameter within of the best value Sloane recorded. Our algorithm was especially successful for smaller numbers of subspaces, but its performance began to flag as the number of subspaces approached 20.
Table B contains several anomalies. For example, our configurations of subspaces in yield worse packing diameters than the configuration of 17 subspaces. It turns out that this configuration of 17 subspaces is optimal, and Sloane’s data show that the (putative) optimal arrangements of 11 to 16 subspaces are all subsets of this configuration. This is the same problem that occurred in some of our earlier experiments, and it suggests again that our algorithm has difficulty locating these degenerate configurations precisely.
The literature contains very few experimental results on packing in complex Grassmannian manifolds equipped with chordal distance. To our knowledge, the only numerical work appears in two short tables from [ARU01]. Therefore, we found it valuable to compare our results against the Rankin bound for subspace packings, which is derived in Section 4. For reference, this bound requires that every configuration of subspaces in satisfy the inequality
This bound cannot always be met. In particular, the bound is attainable in the complex setting only if . In the real setting, the bound requires that . When the bound is attained, each pair of subspaces in is equidistant.
We performed some ad hoc experiments to construct a table of packings in equipped with the chordal distance. For each triple , we constructed random starting points using Algorithm 4 with (which represents no constraint). Then we used the Rankin bound to calculate a lower limit on the feasibility parameter . For this value of , we executed the alternating projection, Algorithm 1, for 5000 iterations.
The best packing diameters we obtained are listed in Table B. We see that there is a remarkable correspondence between the squared packing diameters of our configurations and the Rankin bound. Indeed, many of our packings are within of the bound, which means that these configurations are essentially optimal. The algorithm was less successful as approached , which is an upper bound on the number of subspaces for which the Rankin bound is attainable.
Figure 2 compares the packing diameters of the best configurations in real and complex Grassmannian spaces equipped with chordal distance. It is remarkable that both real and complex packings almost meet the Rankin bound for all where it is attainable. Notice how the real packing diameters fall off as soon as exceeds . In theory, a complex configuration should always attain a better packing diameter than the corresponding real configuration because the real Grassmannian space can be embedded isometrically into the complex Grassmannian space. The figure shows that our best arrangements of 17 and 18 subspaces in are actually slightly worse than the real arrangements calculated by Sloane. This indicates a failure of the alternating projection algorithm.
5.3. The Spectral Distance
Next, we consider how to compute Grassmannian packings with respect to the spectral distance. This investigation requires some small modifications to the algorithm, which are described in the next subsection. Afterward, we provide the results of some numerical experiments.
5.3.1. Modifications to Algorithm
To construct packings with respect to the spectral distance, we tread a familiar path. Suppose that we wish to produce a configuration of subspaces in with a packing diameter . The feasibility problem requires that
where . This leads to the convex structural constraint set
The spectral constraint set is the same as before. The next proposition shows how to find the matrix in closest to an initial matrix. In preparation, define the truncation operator for numbers, and extend it to matrices by applying it to each component.
Let be an Hermitian matrix. With respect to the Frobenius norm, the unique matrix in nearest to has a block identity diagonal. If the off-diagonal block has a singular value decomposition , then
To determine the off-diagonal block of the solution matrix , we must solve the optimization problem
The Frobenius norm is strictly convex and the spectral norm is convex, so this problem has a unique solution.
Let return the vector of decreasingly ordered singular values of a matrix. Suppose that has the singular value decomposition . The constraint in the optimization problem depends only on the singular values of , and so the Hoffman–Wielandt Theorem for singular values [HJ85] allows us to check that the solution has the form .
To determine the singular values of the solution, we must solve the (strictly) convex program
An easy application of Karush–Kuhn–Tucker theory [Roc70] proves that the solution is obtained by truncating the singular values of that exceed . ∎
5.3.2. Numerical Results
To our knowledge, there are no numerical studies of packing in Grassmannian spaces equipped with spectral distance. To gauge the quality of our results, we compare them against the upper bound of Corollary 7. In the real or complex setting, a configuration of subspaces in with respect to the spectral distance must satisfy the bound
In the real case, the bound is attainable only if , while attainment in the complex case requires that [LS73]. When a configuration meets the bound, the subspaces are not only equidistant but also equi-isoclinic. That is, all principal angles between all pairs of subspaces are identical.
We performed some limited ad hoc experiments in an effort to produce good configurations of subspaces with respect to the spectral distance. We constructed random starting points using the modified version of Algorithm 4 with , which represents no constraint. (Again, we did not find that smaller values of improved the performance of the alternating projection.) From the Rankin bound, we calculated the smallest possible value of the feasibility parameter . For values of ranging from the minimal value to twice the minimal value, we ran the alternating projection, Algorithm 1, for 5000 iterations, and we recorded the best packing diameters that we obtained.
Table B displays the results of our calculations. We see that some of our configurations essentially meet the Rankin Bound, which means that they are equi-isoclinic. It is clear that alternating projection also succeeds reasonably well for this packing problem.
The most notable pathology in the table occurs for configurations of 8 and 9 subspaces in . In these cases, the algorithm always yielded arrangements of subspaces with a zero packing diameter, which implies that two of the subspaces intersect nontrivially. Nevertheless, we were able to construct random starting points with a nonzero packing diameter, which means that the algorithm is making the initial configuration worse. We do not understand the reason for this failure.
Figure 3 makes a graphical comparison between the real and complex subspace packings. On the whole, the complex packings are much better than the real packings. For example, every configuration of subspaces in nearly meets the Rankin bound, while just two of the real configurations achieve the same distinction. In comparison, it is curious how few arrangements in come anywhere near the Rankin bound.
5.4. The Fubini–Study Distance
When we approach the problem of packing in Grassmannian manifolds equipped with the Fubini–Study distance, we are truly out in the wilderness. To our knowledge, the literature contains neither experimental nor theoretical treatments of this question. Moreover, we are not presently aware of general upper bounds on the Fubini–Study packing diameter that we might use to assay the quality of a configuration of subspaces. Nevertheless, we attempted a few basic experiments. The investigation entails some more modifications to the algorithm, which are described below. Afterward, we go over our experimental results. We view this work as very preliminary.
5.4.1. Modifications to Algorithm
Suppose that we wish to construct a configuration of subspaces whose Fubini–Study packing diameter exceeds . The feasibility condition is
where . This leads to the structural constraint set
Unhappily, this set is no longer convex. To produce a nearest matrix in , we must solve a nonlinear programming problem. The following proposition describes a numerically favorable formulation.
Let be an Hermitian matrix. Suppose that the off-diagonal block has singular value decomposition . Let , and find a (real) vector that solves the optimization problem
In Frobenius norm, a matrix from that is closest to has a block-identity diagonal and off-diagonal blocks
We use to denote the componentwise exponential of a vector. One may establish that the optimization problem is not convex by calculating the Hessian of the objective function.
To determine the off-diagonal block of the solution matrix , we must solve the optimization problem
We may reformulate this problem as
A familiar argument proves that the solution matrix has the same left and right singular vectors as . To obtain the singular values of the solution, we consider the mathematical program
Change variables to complete the proof. ∎
5.4.2. Numerical Experiments
We implemented the modified version of Algorithm 1 in Matlab, using the built-in nonlinear programming software to solve the optimization problem required by the proposition. For a few triples , we ran 100 to 500 iterations of the algorithm for various values of the feasibility parameter . (Given the exploratory nature of these experiments, we found that the implementation was too slow to increase the number of iterations.)
The results appear in Table B. For small values of , we find that the packings exhibit the maximum possible packing diameter , which shows that the algorithm is succeeding in these cases. For larger values of , we are unable to judge how close the packings might decline from optimal.
Figure 4 compares the quality of our real packings against our complex packings. In each case, the complex packing is at least as good as the real packing, as we would expect. The smooth decline in the quality of the complex packings suggests that there is some underlying order to the packing diameters, but it remains to be discovered.
To perform large-scale experiments, it will probably be necessary to tailor an algorithm that can solve the nonlinear programming problems more quickly. It may also be essential to implement the alternating projection in a programming environment more efficient than Matlab. Therefore, a detailed study of packing with respect to the Fubini–Study distance must remain a topic for future research.
6.1. Subspace Packing in Wireless Communications
Configurations of subspaces arise in several aspects of wireless communication, especially in systems with multiple transmit and receive antennas. The intuition behind this connection is that the transmitted and received signals in a multiple antenna system are connected by a matrix transformation, or matrix channel.
Subspace packings occur in two wireless applications: noncoherent communication and in subspace quantization. The noncoherent application is primarily of theoretical interest, while subspace quantization has a strong impact on practical wireless systems. Grassmannian packings appear in these situations due to an assumption that the matrix channel should be modeled as a complex Gaussian random matrix.
In the noncoherent communication problem, it has been shown that, from an information-theoretic perspective, under certain assumptions about the channel matrix, the optimum transmit signal corresponds to a packing in where corresponds to the minimum of the number of transmit and receive antennas and corresponds to the number of consecutive samples over which the channel is constant [ZT02, HM00]. In other words, the number of subspaces is determined by the system configuration, while is determined by the carrier frequency and the degree of mobility in the propagation channel.
On account of this application, several papers have investigated the problem of finding packings in Grassmannian manifolds. One approach for the case of is presented in [HM00]. This paper proposes a numerical algorithm for finding line packings, but it does not discuss its properties or connect it with the general subspace packing problem. Another approach, based on discrete Fourier transform matrices, appears in [HMR00]. This construction is both structured and flexible, but it does not lead to optimal packings. The paper [ARU01] studies Grassmannian packings in detail, and it contains an algorithm for finding packings in the complex Grassmannian manifold equipped with chordal distance. This algorithm is quite complex: it uses surrogate functionals to solve a sequence of relaxed nonlinear programs. The authors tabulate several excellent chordal packings, but it is not clear whether their method generalizes to other metrics.
The subspace quantization problem also leads to Grassmannian packings. In multiple-antenna wireless systems, one must quantize the dominant subspace in the matrix communication channel. Optimal quantizers can be viewed as packings in , where is the dimension of the subspace and is the number of transmit antennas. The chordal distance, the spectral distance, and the Fubini–Study distance are all useful in this connection [LHJ05, LJ05]. This literature does not describe any new algorithms for constructing packings; it leverages results from the noncoherent communication literature. Communication strategies based on quantization via subspace packings have been incorporated into at least one recent standard [Wir05].
We have shown that the alternating projection algorithm can be used to solve many different packing problems. The method is easy to understand and to implement, even while it is versatile and powerful. In cases where experiments have been performed, we have often been able to match the best packings known. Moreover, we have extended the method to solve problems that have not been studied numerically. Using the Rankin bounds, we have been able to show that many of our packings are essentially optimal. It seems clear that alternating projection is an effective numerical algorithm for packing.
Appendix A Tammes’ Problem
The alternating projection method can also be used to study Tammes’ Problem of packing points on a sphere [Tam30]. This question has received an enormous amount of attention over the last 75 years, and extensive tables of putatively optimal packings are available [Slo04b]. This appendix offers a brief treatment of our work on this problem.
a.1. Modifications to Algorithm
Suppose that we wish to produce a configuration of points on the unit sphere with a packing diameter . The feasibility problem requires that
where . This leads to the convex structural constraint set
The spectral constraint set is the same as before. The associated matrix nearness problem is trivial to solve.
Let be a real, symmetric matrix. With respect to Frobenius norm, the unique matrix in closest to has a unit diagonal and off-diagonal entries that satisfy
a.2. Numerical Results
Tammes’ Problem has been studied for 75 years, and many putatively optimal configurations are available. Therefore, we attempted to produce packings whose maximum inner product fell within of the best value tabulated by N. J. A. Sloane and his colleagues [Slo04b]. This resource draws from all the experimental and theoretical work on Tammes’ Problem, and it should be considered the gold standard.
Our experimental setup echoes the setup for real projective packings. We implemented the algorithms in Matlab, and we performed the following experiment for pairs with and . First, we computed the putatively optimal maximum inner product using the data from [Slo04b]. In each of 10 trials, we constructed a starting matrix using Algorithm 4 with parameters and . Then, we executed the alternating projection, Algorithm 1, with the calculated value of and the maximum number of iterations set to . We stopped the alternating projection in Step 4 if the iterate contained no off-diagonal entry greater than and proceeded with Step 6. After 10 trials, we recorded the largest packing diameter attained, as well as the average value of the packing diameter. We also recorded the average number of iterations the alternating projection required during each trial.
Table B provides the results of this experiment. The most striking feature of Table B is that the best configurations returned by alternating projection consistently attain packing diameters that fall hundredths or thousandths of a degree away from the best packing diameters recorded by Sloane. If we examine the maximum inner product in the configuration instead, the difference is usually on the order of or , which we expect based on our stopping criterion. The average-case results are somewhat worse. Nevertheless, the average configuration returned by alternating projection typically attains a packing diameter only several tenths of a degree away from optimal.
A second observation is that the alternating projection algorithm typically performs better when the number of points is small. The largest errors are all clustered at larger values of . A corollary observation is that the average number of iterations per trial tends to increase with the number of points. We believe that the explanation for these phenomena is that Tammes’ Problem has a combinatorial regime, where solutions have a lot of symmetry and structure, and a random regime, where the solutions have very little order. The algorithm typically seems to perform better in the combinatorial regime, although it fails for certain unusually structured ensembles.
This claim is supported somewhat by theoretical results for . Optimal configurations have only been established for and . Of these, the cases are trivial. The cases fall from the vertices of various well-known polyhedra. The cases are degenerate, obtained by leaving a point out of the solutions for . The remaining cases involve complicated constructions based on graphs [EZ01]. The algorithm was able to calculate the known optimal configurations to a high order of accuracy, but it generally performed slightly better for the nondegenerate cases.
On the other hand, there is at least one case where the algorithm failed to match the optimal packing diameter, even though the optimal configuration is highly symmetric. The best arrangement of 24 points on locates them at vertices of a very special polytope called the 24-cell [Slo04b]. The best configuration produced by the algorithm has a packing diameter worse. It seems that this optimal configuration is very difficult for the algorithm to find. Less dramatic failures occurred at pairs , , , , and . But in each of these cases, our best packing declined more than a tenth of a degree from the best recorded.
Appendix B Tables and Figures
Our experiments resulted in tables of packing diameters. We did not store the configurations produced by the algorithm. The Matlab code that produced these data is available on request from firstname.lastname@example.org.
These tables and figures are intended only to describe the results of our experiments; it is likely that many of the packing diameters could be improved with additional effort. In all cases, we present the results of calculations for the stated problem, even if we obtained a better packing by solving a different problem. For example, a complex packing should always improve on the corresponding real packing. If the numbers indicate otherwise, it just means that the complex experiment yielded an inferior result. As a second example, the optimal packing diameter must not decrease as the number of points increases. When the numbers indicate otherwise, it means that running the algorithm with more points yielded a better result than running it with fewer. These failures may reflect the difficulty of various packing problems.
List of Tables
- 1 Packing in real projective spaces
- 2 Packing in complex projective spaces
- 3 Packing in real Grassmannians with chordal distance
- 4 Packing in complex Grassmannians with chordal distance
- 5 Packing in Grassmannians with spectral distance
- 6 Packing in Grassmannians with Fubini–Study distance
- 7 Packing on spheres
List of Figures
|Packing diameters (Degrees)||Iterations|
|NJAS||Best of 10||Error||Avg. of 10||Error||Avg. of 10|
|Packing diameters (Degrees)|
|Squared Packing diameters|