Barycentric Subspace Analysis (BSA)

• Barycentric Subspace Analysis (BSA) is a dimensionality reduction method on metric and Riemannian spaces that generates subspaces using weighted reference points.
• It generalizes traditional affine spans and principal geodesic analysis by optimizing accumulated unexplained variance across hierarchically nested subspaces.
• Recent extensions adapt BSA to network-valued data and barycentric interpolation on manifolds, offering interpretable and efficient representations for complex datasets.

Searching arXiv for the cited BSA papers to ground the article in current preprints and confirm metadata. arXiv search query: "(Pennec, 2016) Barycentric Subspace Analysis on Manifolds" Barycentric Subspace Analysis (BSA) is a dimensionality reduction method on metric and Riemannian spaces in which the approximating “subspace” is generated by points rather than by vectors. For a fixed set of anchors, the associated barycentric subspace is the locus of weighted Fréchet or Karcher means of those anchors under appropriate weight constraints. In Euclidean spaces, this construction recovers affine spans or convex hulls; on manifolds, it generalizes geodesic subspaces, supports hierarchically nested flags, and yields a global optimization criterion based on accumulated unexplained variance (AUV). Later work has adapted the framework to unlabeled network-valued data through spectral graph spaces, while related barycentric constructions on manifolds of linear subspaces have been used for parametric reduced-order modeling (Pennec, 2016, Maignant et al., 31 Jul 2025, Oulghelou et al., 2020).

1. Foundational definition and barycentric coordinates

The basic object in BSA is a subspace generated by reference points. For $k+1$ reference points, barycentric weights are taken in the projective set

$$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$

with normalized weights

$$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$

The framework allows both convex weights and general signed affine weights. This point-based formulation is the defining distinction between BSA and PCA-type methods generated by vectors (Pennec, 2016).

In Euclidean space, barycentric subspaces coincide with ordinary affine spans. For affinely independent points $\{p_0,\dots,p_k\}\subset\mathbb{R}^n$,

$$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$

equivalently characterized by the barycentric balance

$\sum_{i=0}^k\lambda_i(p_i-x)=0.$

This Euclidean identity is the prototype for the manifold construction.

On a geodesically complete Riemannian manifold $(\mathcal{M},g)$, the formulation is transferred from linear differences to logarithmic displacements. For reference points $\{x_0,\dots,x_k\}$, one removes their cut loci and works on

$$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$

where $\log_x(x_i)$ is well-defined and smooth. The first weighted moment is then

$$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$0

A central consequence is that the construction depends on distances and points rather than on a fixed tangent space. The original formulation explicitly notes that this makes it extensible to geodesic spaces which are not Riemannian and, in particular, to stratified spaces where previous generalizations of PCA cannot span several strata (Pennec, 2016).

2. Exponential, Fréchet, and Karcher barycentric subspaces

The most direct manifold analogue of an affine span is the Exponential Barycentric Subspace (EBS). Given references $$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$1,

$$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$2

If

$$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$3

then the EBS condition is $$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$4 for some $$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$5. The affine span in a manifold is defined as the metric completion of the EBS,

$$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$6

Fréchet and Karcher barycentric subspaces refine this critical-point definition by imposing minimization of a weighted variance. With

$$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$7

the Fréchet barycentric subspace is the set of global minimizers over $$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$8 and $$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$9, while the Karcher barycentric subspace is defined analogously with local minima. On $$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$0, critical points satisfy

$$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$1

Accordingly, the inclusion

$$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$2

holds on $$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$3 (Pennec, 2016).

The local structure of EBS is controlled by the rank of $$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$4 and the Hessian

$$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$5

The nondegenerate EBS generated by $$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$6 affinely independent points is a stratified space of dimension $$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$7 on $$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$8; on an $$\underline{\lambda}_i=\frac{\lambda_i}{\sum_{j=0}^k\lambda_j},\qquad \sum_{i=0}^k\underline{\lambda}_i=1.$$9-dimensional stratum, $\{p_0,\dots,p_k\}\subset\mathbb{R}^n$0 has exactly $\{p_0,\dots,p_k\}\subset\mathbb{R}^n$1 vanishing singular values. Positive points, where $\{p_0,\dots,p_k\}\subset\mathbb{R}^n$2, form the nondegenerate Karcher barycentric subspace. Equivalent characterizations use

$\{p_0,\dots,p_k\}\subset\mathbb{R}^n$3

so that

$\{p_0,\dots,p_k\}\subset\mathbb{R}^n$4

equivalently the smallest eigenvalue of $\{p_0,\dots,p_k\}\subset\mathbb{R}^n$5 vanishes, or the $\{p_0,\dots,p_k\}\subset\mathbb{R}^n$6th eigenvalue of $\{p_0,\dots,p_k\}\subset\mathbb{R}^n$7 vanishes.

Two geometric limits are especially important. First, when normalized weights are restricted to be nonnegative, the resulting barycentric simplex contains the references, their geodesic edges, and the Fréchet mean. Second, if the references coalesce toward a basepoint along tangent directions, barycentric subspaces converge to restricted geodesic submanifolds. This establishes the formal statement that EBS generalizes geodesic subspaces rather than merely coexisting with them. In constant-curvature spaces the geometry becomes explicit: on spheres, $\{p_0,\dots,p_k\}\subset\mathbb{R}^n$8 is the great $\{p_0,\dots,p_k\}\subset\mathbb{R}^n$9-subsphere containing the references, while on hyperboloids it is the $$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$0-dimensional hyperboloid induced by the ambient Minkowski hyperplane containing the references (Pennec, 2016).

3. Flags, accumulated unexplained variance, and the PCA generalization

BSA is not only a family of point-generated subspaces; it is also a flag optimization framework. In Euclidean space, a flag is a strictly increasing sequence of affine subspaces

$$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$1

For data $$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$2, the unexplained variance of a subspace $$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$3 is

$$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$4

where $$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$5, and the accumulated unexplained variance is

$$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$6

The original theory shows that Euclidean PCA minimizes this AUV over flags of nested affine subspaces. Under simple multiplicity of the first $$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$7 eigenvalues of the covariance matrix, the flag generated by the mean and the first $$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$8 principal directions is the AUV minimizer (Pennec, 2016).

The manifold generalization replaces linear subspaces by flags of affine spans generated by ordered reference points: $$\mathrm{Aff}(p_0,\dots,p_k)=\left\{x=\sum_{i=0}^k\lambda_i p_i\ :\ \sum_{i=0}^k\lambda_i=1\right\},$$9 For a submanifold $\sum_{i=0}^k\lambda_i(p_i-x)=0.$0,

$\sum_{i=0}^k\lambda_i(p_i-x)=0.$1

and the BSA problem becomes

$\sum_{i=0}^k\lambda_i(p_i-x)=0.$2

This is the point-based analogue of PCA on manifolds. Its principal output is a hierarchy of properly embedded barycentric subspaces rather than a single tangent-space linearization.

A common misconception is that BSA is only a reparameterization of tangent PCA or principal geodesic analysis. The original formulation distinguishes them sharply. Tangent PCA unfolds data into one tangent space and performs Euclidean PCA there. Principal Geodesic Analysis minimizes squared distances to subspaces spanned by geodesics through a point with tangents in a linear subspace. BSA instead uses multiple reference points, need not pass through the mean, admits forward and backward nested sequences, and optimizes a global AUV criterion across dimensions. The limit theorem shows that PGA-type restricted geodesic submanifolds arise when the reference points coalesce, so the relationship is one of generalization, not equivalence. The same source also contrasts BSA with Geodesic PCA and Principal Nested Spheres, emphasizing that BSA naturally accommodates stratified geodesic spaces and can span several strata (Pennec, 2016).

4. Algorithms, projection, and numerical issues

Computationally, BSA can be formulated either as a continuous optimization problem over reference points or as a sample-limited problem in which the references are chosen from the dataset. The sample-limited form is algorithmically simple but combinatorial: exhaustive enumeration scales as $\sum_{i=0}^k\lambda_i(p_i-x)=0.$3, which is practical only for small $\sum_{i=0}^k\lambda_i(p_i-x)=0.$4. Continuous optimization over $\sum_{i=0}^k\lambda_i(p_i-x)=0.$5 instead shifts the cost to repeated evaluations of logarithm maps, projection solvers, and derivatives of $\sum_{i=0}^k\lambda_i(p_i-x)=0.$6 (Pennec, 2016).

Affine independence is a structural prerequisite. For points $\sum_{i=0}^k\lambda_i(p_i-x)=0.$7, one requires that for each $\sum_{i=0}^k\lambda_i(p_i-x)=0.$8, the $\sum_{i=0}^k\lambda_i(p_i-x)=0.$9 vectors $(\mathcal{M},g)$0 are linearly independent; on spheres and hyperboloids this is equivalent to $(\mathcal{M},g)$1. Once references are fixed, computing the EBS amounts to solving

$(\mathcal{M},g)$2

Two standard routes are described: fixing $(\mathcal{M},g)$3 and solving $(\mathcal{M},g)$4 by Newton’s method on the manifold, with Jacobian

$(\mathcal{M},g)$5

or fixing $(\mathcal{M},g)$6 and extracting $(\mathcal{M},g)$7 from the right singular vectors associated to vanishing singular values.

Projection onto a barycentric subspace is defined by

$(\mathcal{M},g)$8

In general manifolds, this is handled by iterative geodesic descent. In constant-curvature spaces, explicit ambient formulas are available. On the sphere,

$(\mathcal{M},g)$9

with the sign chosen to minimize geodesic distance. On the hyperboloid, the ambient Minkowski projection is normalized back to the hyperboloid. These explicit formulas are one reason spheres and hyperbolic spaces serve as the canonical examples in the original development.

Optimization over flags can proceed forward, by adding one point at a time, or backward, by starting from a high-dimensional affine span and removing points. The theory also notes several recurrent numerical issues: non-uniqueness, since multiple tuples can parametrize the same subspace; cut loci, where $\{x_0,\dots,x_k\}$0 becomes multivalued; sensitivity to initialization in the nonconvex AUV objective; and ill-conditioning when curvature is small or references coalesce. For nonnegative weights inside a strongly convex geodesic ball with radius

$\{x_0,\dots,x_k\}$1

the weighted Fréchet mean is unique, and the barycentric simplex is the graph of a $\{x_0,\dots,x_k\}$2-dimensional differentiable function of the normalized weights. This provides the main local existence and uniqueness guarantee used by convex variants of the construction (Pennec, 2016).

5. Network-valued data and spectral graph spaces

A recent extension applies BSA to unlabeled network-valued data by changing the ambient space rather than the core barycentric principle. The setting is the spectral graph space

$\{x_0,\dots,x_k\}$3

where symmetric adjacency matrices are quotiented by orthogonal conjugation. This replaces equivalence classes of isomorphic networks by equivalence classes of cospectral networks. For a graph $\{x_0,\dots,x_k\}$4 with adjacency matrix $\{x_0,\dots,x_k\}$5, the relevant coordinates are the sorted eigenvalues

$\{x_0,\dots,x_k\}$6

The quotient distance becomes

$\{x_0,\dots,x_k\}$7

so distances, geodesics, and the log map reduce to Euclidean operations on sorted spectra (Maignant et al., 31 Jul 2025).

Within $\{x_0,\dots,x_k\}$8, barycentric subspaces become convex polytopes in the sorted cone. For anchors $\{x_0,\dots,x_k\}$9, the barycentric subspace is the locus of barycenters over affine weights satisfying sorted-eigenvalue feasibility. Equivalently,

$$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$0

If the anchor spectra are affinely independent, the dimension equals $$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$1. In the convex variant, one also imposes $$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$2, yielding the convex hull of the anchor spectra intersected with the sorted cone.

Projection of a network with spectrum $$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$3 onto such a subspace is a quadratic program in the barycentric weights: $$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$4 plus the linear inequalities enforcing sortedness, and $$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$5 in convex BSA. This yields an explicitly interpretable representation when anchors are selected among the observed networks. The paper terms this sample-limited BSA and emphasizes that projections are then barycentric combinations of actual networks rather than linear directions that may correspond to non-existent topologies.

The empirical comparison with tangent PCA is framed primarily around interpretability, but it also reports reconstruction errors. On a two-parameter dataset with $$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$6 and $$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$7, the squared projection error is $$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$8 for tangent PCA, $$\mathcal{M}^*=\mathcal{M}\setminus\bigcup_{i=0}^k C(x_i),$$9 for sample-limited BSA, and $\log_x(x_i)$0 for convex BSA. On a clustered dataset with $\log_x(x_i)$1 and $\log_x(x_i)$2, backward convex BSA selected one anchor per cluster and showed a clear elbow at a 2D subspace. On a real airlines dataset with $\log_x(x_i)$3 and $\log_x(x_i)$4 macro-regions, backward convex BSA suggested a 1D subspace, with anchors corresponding to easyJet and Swiss and projections interpreted as a ranking by centralization level. These experiments support the claim that BSA can preserve existing topologies and produce interpretable feature subspaces when structural variability is discrete and nonlinear (Maignant et al., 31 Jul 2025).

The same framework has clear limitations. Cospectral equivalence is coarser than isomorphism, so non-isomorphic graphs can collide. Anchor selection is sensitive: poor anchors reduce both interpretability and accuracy. Reconstruction in the adjacency domain may introduce small self-loops. The paper presents convexity constraints, diversity heuristics, backward selection, and possible future use of Laplacian or normalized Laplacian spectra as mitigations rather than complete resolutions (Maignant et al., 31 Jul 2025).

6. Relation to barycentric interpolation on manifolds of subspaces

BSA is conceptually linked to later work on barycentric interpolation of linear subspaces, although the objectives are different. Both rely on barycenters, or Fréchet and Karcher means, defined by a Riemannian metric on a manifold of subspaces; in both, principal angles and invariance under right-orthogonal transformations are foundational, and barycentric combinations seek subspaces minimizing weighted sums of squared Riemannian distances. The reduced-order modeling application in fluid dynamics makes this explicit on

$\log_x(x_i)$5

which is smoothly equivalent to the Grassmann manifold $\log_x(x_i)$6 (Oulghelou et al., 2020).

In that setting, the chosen quotient geometry is not the canonical Grassmann geodesic metric. For $\log_x(x_i)$7 with $\log_x(x_i)$8 nonsingular, if

$\log_x(x_i)$9

then

$$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$00

The Karcher condition for a weighted barycenter yields the fixed-point equation

$$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$01

with iteration

$$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$02

This construction is used to interpolate POD subspaces across Reynolds numbers and assemble parametric reduced-order Navier–Stokes models.

The conceptual similarity to BSA is direct, but the scope is different. BSA builds nested families of barycentric subspaces for exploratory data analysis and manifold dimension reduction. The reduced-order modeling paper uses a single Karcher barycenter to interpolate across a parameter domain, producing a parametric subspace for projection-based ROM; it is not a hierarchical data-analytic method. Numerically, the barycentric PROM matches the accuracy of ITSGM on the reported cylinder and lid-driven cavity tests while reducing online update time by about three orders of magnitude: for the cylinder, operator assembly is approximately $$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$03 for the barycentric PROM versus approximately $$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$04 for ITSGM, and for the lid-driven cavity it is approximately $$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$05 versus approximately $$\mathcal{P}_k=\left\{\lambda=(\lambda_0:\cdots:\lambda_k)\in \mathbb{P}^k\ \text{s.t.}\ \sum_{i=0}^k\lambda_i\neq 0\right\},$$06 (Oulghelou et al., 2020).

This suggests that barycentric subspace constructions are not confined to statistical dimension reduction. The common mathematical ground is the Karcher mean and the distance induced by a chosen Riemannian structure, while the divergence lies in objective: fast parametric ROM construction in one case, and exploratory manifold data analysis with nested barycentric subspaces in the other (Oulghelou et al., 2020).