Anchor PCA

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Anchor PCA

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Anchor PCA

arXiv:2606.06233v1 [stat.ML] 4 Jun 2026

Benedikt Seiter ETH Zürich [email protected]

Anya Fries Seminar for Statistics ETH Zürich

Julius von Kügelgen Seminar for Statistics ETH Zürich

Jonas Peters Seminar for Statistics ETH Zürich

Abstract Principal component analysis (PCA) is one of the most widely used unsupervised dimension reduction techniques. We study PCA for data from multiple related domains. Since principal components generally differ across domains, one way to obtain a shared low-rank embedding is to perform PCA on the pooled data. However, this approach can focus on spurious directions that exhibit high variation in only a few domains. To find a robust embedding that still explains most variance in unseen but similar domains, we propose instead to focus on shared directions of variation. To this end, we introduce Anchor PCA which trades off overall explained variance with agreement between the shared and domain-specific low-rank embeddings. Anchor PCA amounts to PCA on a modified target matrix and thus can be solved efficiently. Moreover, we show that Anchor PCA recovers a maximal invariant subspace and admits a minimax reconstruction interpretation under bounded domain-specific covariance inflations. On simulated and real-world gas sensor data with temporal drift, we demonstrate, respectively, that Anchor PCA recovers the maximally invariant subspace and yields embeddings that explain more variance on unseen domains than the pooling baseline and a worst-case alternative. Taken together, these findings establish Anchor PCA as a promising approach to robust unsupervised dimension reduction from multi-domain data.

Introduction

Principal component analysis (PCA) [20] is one of the most widely used methods for unsupervised dimensionality reduction. Given a dataset, PCA finds a low-dimensional subspace that captures as much variance as possible. It is used extensively in exploratory data analysis, preprocessing, and representation learning across the natural sciences, engineering, and machine learning [e.g., 24, 27, 40, 46]. Here, we focus on the following two use cases of PCA: providing interpretable directions (that may carry physical meaning) and explaining variance. We assume that we have data collected across multiple related but heterogeneous settings—different geographic locations, time points or experimental conditions. We refer to these settings as source domains. If the principal subspaces differ substantially across the source domains, we may achieve neither of the two goals mentioned above. For example, PCA applied to the average covariance (poolPCA) could be dominated by a direction that concentrates variance in only one or two domains. This direction may explain very little variance elsewhere and thus may not yield physically relevant directions; similarly, it may not explain a lot of variance in a target domain. In this paper, we propose to focus not only on explained variance but also on finding invariant directions. This builds on work in multi-domain dimension reduction and invariant learning. In dimension reduction, recent approaches optimize worst-case performance across domains [15, 45] or enforce invariance of distributions [26]. In supervised learning, invariance of predictive relations or representations across domains is widely studied [3, 44]. In contrast, we introduce a geometric notion of invariance based on agreement of domain-specific principal subspaces, and explicitly trade this against explained variance. Similar in spirit is anchor regression [33], which introduces a soft invariPreprint.

Table 1: Eigendecomposition and reconstruction errors for the motivating example. If k = 3, the invariant directions are a = c1 and b = c2 ; u, v, and w are domain-specific directions in the c3 –c4 plane, and u⊥ , v ⊥ , w⊥ are their orthogonal complements in that plane. The domain-specific columns list eigenvectors and eigenvalues, equivalently the explained variance of each direction. The invariant directions do not need fixed ranks: a is top-ranked in domain 2, while b outranks a in domain 3. The perturbed column evaluates the same reconstruction-error curves shown later in Fig. 2 at the dotted vertical line. Rank

Domain 1

Domain 2

Domain 3

Method

direction value direction value direction value 1 2 3 4

u a b u⊥

220 140 90 25

a v b v⊥

120 90 70 10

w b a w⊥

320 120 80 10

Avg. reconstr. error (↓) Original Perturbed covariances covariances

poolPCA AnchorPCAλ=25 AnchorPCA∞

93.3 98.5 113.8

243.3 163.4 172.5

ance penalty in the supervised setting. Concretely, we propose AnchorPCAλ (with λ ∈ R>0 ∪ {∞}), a method that finds a k-dimensional subspace where the parameter λ trades off between explaining variance and finding invariant directions (e.g., for λ = 0, AnchorPCAλ reduces to poolPCA). We prove that AnchorPCA satisfies two attractive theoretical properties. First, if there exists a nontrivial m-dimensional invariant subspace S⋆ , that is, a subspace that lies in the top-k principal subspace of every domain, then AnchorPCA∞ is guaranteed to contain it, and the leading directions of AnchorPCAλ converge to S⋆ as λ → ∞ (Thm. 3.1). Second, AnchorPCAλ admits a minimax reconstruction interpretation: it is the rank-k projector that minimizes worst-case average reconstruction error over test domains that are subject to domain-specific covariance inflations (Thm. 4.1). We now introduce the setting and build intuition through the motivating example in § 1.2. § 2 introduces Anchor PCA, § 3 and § 4 contain the main theoretical results, § 5 shows learning methods for finite data, and § 6 contains simulated and real-world experiments. § A, § B, § C, and § D contain further related work, further details on Anchor PCA, proofs, and experiment details, respectively. 1.1

Setting and notation

We consider a set of E source domains indexed by e ∈ E := {1, . . . , E}. For each domain e, let xe ∈ Rp be a zero-mean random (column) vector with covariance Σe ∈ Sp+ where Sp+ := {Σ ∈ Rp×p : Σ = Σ⊤ , Σ ⪰ 0} is the set of symmetric positive semidefinite matrices and we further assume Tr(Σe ) > 0. Recall that PCA selects the eigenvectors of a covariance matrix Σ corresponding to its k largest eigenvalues, also called a top-k eigenbasis. We denote a top-k eigenbasis of Σe by We ∈ Op×k , where Op×k := {W ∈ Rp×k : W ⊤ W = Ik } is the Stiefel manifold of p × k matrices (e) with orthonormal columns. The rank-k principal projection matrices are given by Πk := We We⊤ , ⊤ which can be viewed as a linear autoencoder with encoder We and decoder We . Given the domain-specific covariance matrices (Σe )e∈E (or i.i.d. samples from the corresponding distributions), we seek to solve the following problem, which we formalize in § 3 and § 4. Find a shared k-dimensional subspace that is robust with respect to changing domains. PE PE (e) Denote the average covariance by Σ := E1 e=1 Σe and the average projector by Π := E1 e=1 Πk . ⊤ For general W ∈ Op×k , we write ΠW := W W for the projector onto the space spanned by W . p The Frobenius and operator norms of a matrix X ∈ Rp×k are defined as ∥X∥F := Tr(X ⊤ X) and ∥X∥op := sup∥v∥2 =1 ∥Xv∥2 , respectively. For a linear subspace S ≤ Rp , let ΠS denote the orthogonal projector onto S; if U has orthonormal columns spanning S, then ΠS = U U ⊤ . 1.2

Motivating example

We consider a simple example in R4 (so p = 4) with canonical basis c1 , c2 , c3 , c4 , writing a := c1 and b := c2 , where we seek a k = 3-dimensional representation. We observe E = 3 domains. Table 1 summarizes the eigenvectors and eigenvalues of the covariance matrices (we defer the details to § D.1). 2

Domain 1 Domain 2 Domain 3

poolPCA AnchorPCA = 25 AnchorPCA

1.0

0.5

v u

0.0

0.5

0.0 0.5

1.0 1.0

0.5

0.0

0.5

1.0

0.5

0.0

0.5

1.0

1.0 1.0

0.5

0.0

0.5

1.0

0.5

0.0

0.5

1.0

Figure 1: Geometric view of the motivating example. Here, we ignore the direction a, which is chosen by all methods, yielding a reduced representation in the coordinates (c3 , c4 , b). Left: the three local top-3 eigenspaces, shown as planes through the invariant line. Right: the recovered rank-3 subspaces in the same coordinates. Unlike poolPCA, both AnchorPCAλ=25 and AnchorPCA∞ contain the true invariant direction b and differ only in the third direction (λ = 25 puts more emphasis on explained variance, so AnchorPCAλ=25 is pulled toward the higher variance direction w). Informally, we consider a direction invariant if it helps explain variance in every domain (the domain acts as an ‘anchor’ specifying what we would like to be invariant against [see also 33]). In this example, span(a, b) is an invariant space: in all domains, a and b are among the top-k eigenvectors. Fig. 1 visualizes the subspaces that are chosen by the different methods. The figure shows the qualitative behavior that will recur throughout the paper. poolPCA follows pooled variance and contains the space spanned by the nuisance directions: it keeps a but misses the invariant direction b. By contrast, both AnchorPCAλ=25 and AnchorPCA∞ recover the invariant span exactly. The difference between these two lies in the third direction: AnchorPCAλ=25 puts a stronger focus on the explained variance and is thus pulled more strongly toward w, which is the direction with largest eigenvalue (see Table 1). In this paper, we focus on the subspace and not the basis vectors; that is, if in domain 3, the top-k eigenspace contained ã, b̃, instead of a, b, such that span(a, b) = span(ã, b̃), we would still regard span(a, b) an invariant subspace. The other directions u, v, and w can be considered nuisance directions (in the light of invariance).

Maximal invariant subspace and Anchor PCA

We now formally define the notion of a maximal invariant subspace and the Anchor PCA problem. 2.1

Maximal invariant subspace

Given source domains E with covariances (Σe )e∈E , a natural choice is to ignore the domain structure and solve PCA on the pooled covariance Σ. We refer to this procedure as poolPCA. Formally, a matrix W ⋆ solves rank-k poolPCA if W ⋆ ∈ arg maxW ∈Op×k Tr(W ⊤ ΣW ). Indeed, for any Pk W = [w1 , . . . , wk ] ∈ Op×k , we have Tr(W ⊤ ΣW ) = j=1 wj⊤ Σwj = Tr(ΣΠW ), which is the total variance explained by projecting onto the candidate subspace span(W ); hence, maximizing the trace recovers the k leading eigenvectors of Σ. By definition, poolPCA focuses on the average explained variance, which may result in large differences in explained variance across source domains. In this paper, we propose to focus on directions of variation that are shared among all domains. In particular, we can consider those directions that appear within the top-k principal directions in all domains. Definition 2.1 (Maximal invariant subspace). Let k < p be fixed. The maximal invariant subspace and its dimension are defined as S⋆ :=

E \

(e) Im Πk ,

e=1

m := dim(S⋆ ).

(If there are eigenvalue ties, S⋆ and Π and the AnchorPCA solutions—defined below—depend on a fixed choice of a rank-k principal projector in each domain.) The definition of S⋆ depends on k: in general, the dimension m of S⋆ is at most k, but varies according to the problem at hand. In Example 1.2, m = 2 and S⋆ = span(a, b). Equivalently, we can define S⋆ as the eigenspace of Π asPE (e) sociated with eigenvalue 1; indeed, for every unit vector v, we have v ⊤ Πv = E1 e=1 v ⊤ Πk v ≤ 1, (e) and equality holds if and only if Πk v = v for every e, that is, if and only if v ∈ S⋆ . We will use this characterization when estimating S⋆ from data (§ 5.2). 2.2

Anchor PCA

In practice, we may not expect to find any (nonzero) direction that is exactly invariant, in which case S⋆ = {0} is trivial. Furthermore, we may be interested in finding a subspace of fixed dimension k, rather than one whose dimension can vary depending on the concrete set of covariance matrices. We therefore propose selecting a rank-k subspace that minimizes disagreement with the domain-specific top-k eigenspaces. Such a space may be represented using different bases; we propose to choose one that maximizes pooled explained variance. We refer to this approach as AnchorPCA∞ . Definition 2.2 (AnchorPCA∞ ). A matrix Wk⌣,∞ ∈ Op×k solves rank-k AnchorPCA∞ if Wk⌣,∞ ∈ arg max Tr W ∈Op×k

⊤

W ΣW

subject to:

W ∈ arg min

E X

f ∈Op×k e=1 W

(e)

ΠW f − Πk

. F

(1)

The corresponding rank-k AnchorPCAλ subspace and projection matrix are given by Sk⌣,∞ := ,∞ := Wk⌣,∞ (Wk⌣,∞ )⊤ , respectively. Im(Wk⌣,∞ ) and Π⌣ k We prove in § 4 that AnchorPCA∞ satisfies an asymptotic minimax-optimality guarantee for reconstruction error under additional variance inflations of arbitrary strengths in certain directions. If we only expect mild inflations, however, AnchorPCA∞ is overly conservative. We therefore also consider the following unconstrained version in which the hard invariance constraint is relaxed to a soft penalty. Definition 2.3 (AnchorPCAλ ). For λ > 0, a matrix Wk⌣,λ ∈ Op×k solves rank-k AnchorPCAλ if ) ( E X 2 (e) ,λ ⊤ . (2) Wk⌣ ∈ arg max Tr W ΣW − λ ΠW − Π k W ∈Op×k

e=1

The corresponding rank-k AnchorPCAλ subspace and projection matrix are given by Sk⌣,λ := ,λ Im(Wk⌣,λ ) and Π⌣ := Wk⌣,λ (Wk⌣,λ )⊤ , respectively. k While AnchorPCA∞ only considers maximally invariant solutions and explained variance serves as a secondary objective, AnchorPCAλ balances the two, with the trade-off controlled by the invariance level λ. We prove in § B.2.1 that, as λ grows, the two solutions converge to each other: under suitable assumptions, we have ,λ ,∞ Π⌣ − Π⌣ →0 as λ → ∞. k k op

Thm. C.8 (ii) studies the convergence at the level of directions rather than subspaces and projections. Both Anchor PCA methods can be extended to a soft agreement score: even when S⋆ is too small or trivial, directions with larger Π-eigenvalues are still preferred, rather than requiring exact equality of eigenvalues to one. 2.3

Projection-based formulation

The optimization problems (1) and (2) allow for an equivalent, projection-based view: expanding the PE PE (e) (e) penalty in (1) gives e=1 ∥ΠW − Πk ∥2F = e=1 2k − 2 Tr(ΠW Πk ) = 2Ek − 2E Tr(W ⊤ ΠW ). Hence, the constraint is equivalent to maximizing Tr(W ⊤ ΠW ) over Op×k . Similarly, for (2), solutions are obtained by solving a standard rank-k PCA on a modified target matrix Mλ max Tr(W ⊤ Mλ W ),

W ∈Op×k

where

Mλ := Σ + 2Eλ Π.

(3)

This view has two advantages: (i) It suggests computationally feasible ways to solve the optimization problems of AnchorPCA∞ and AnchorPCAλ (see § 5). (ii) It allows for a geometric interpretation: (e) Tr(ΠW Πk ) is the sum of squared cosines of the principal angles between span(W ) and the principal subspace of domain e. Thus, AnchorPCA∞ first maximizes average overlap with the domain-specific principal subspaces and then maximizes pooled explained variance among those solutions. 2.4

Choosing an ordered basis

Both objectives (1) and (2) depend on W only through ΠW = W W ⊤ , so any orthonormal basis of a solution space defined by W forms another solution; statements about individual columns requires choosing an ordered basis. We propose to use the following convention for choosing an ordered basis (see § B.1 for details). For AnchorPCA∞ , we first order directions by decreasing agreement, that is, by decreasing Π-eigenvalue; this creates blocks of equal eigenvalues (which we refer to in the algorithm and proof sections); we then choose a basis that diagonalizes Σ restricted to the Π-eigenspace corresponding to that block and order the vectors by non-increasing Σ-eigenvalue, equivalently, explained variance. For AnchorPCAλ , we use an ordered top-k eigenbasis of Mλ , see Eq. (3), giving the usual sequential PCA interpretation for the modified covariance Mλ . Both of these conventions are implemented in our code, see § 5. We will see in § 3 and § 4 that both versions of AnchorPCA can find invariant directions and come with improved recovery and robustness properties compared to poolPCA.

Recovery guarantees

We now study when AnchorPCA recovers invariant structure. As motivated in Example 1.2, the (e) subspace target is S⋆ = ∩E e=1 Im(Πk ), the part of the local top-k principal subspaces shared by all domains. Since AnchorPCAλ and AnchorPCA∞ return rank-k subspaces, there are two natural recovery questions: whether these subspaces contain S⋆ , and whether their leading m = dim(S⋆ ) directions span S⋆ . We first prove that the AnchorPCA∞ subspace Sk⌣,∞ contains S⋆ , and that AnchorPCAλ identifies S⋆ asymptotically as λ → ∞. Theorem 3.1 (Invariant subspace identification and asymptotic convergence). Let m = dim(S⋆ ) ≤ k. ,∞ Let Π⌣ be any rank-k AnchorPCA∞ projection matrix. Then k ,∞ ΠS⋆ = ΠS⋆ , Π⌣ k

equivalently

,∞ S⋆ ⊆ Im(Π⌣ ). k

Moreover, for λ > 0, choose the AnchorPCAλ solution Wk⌣,λ = (w1⌣,λ , . . . , wk⌣,λ ) to be an ordered ,λ orthonormal eigenbasis of Mλ associated with its k largest eigenvalues. Let Π⌣ denote the m ,λ ⌣,λ ). Then, there exists λ such that for all λ > λ , the projection matrix onto span(w1⌣ , . . . , wm 0 0 span of these first m eigenvectors is unique and ⌣,λ ) → S , span(w1⌣,λ , . . . , wm ⋆

,λ Π⌣ − ΠS⋆ op → 0 m

i.e.,

as λ → ∞.

Under the ordered basis convention in § B.1, also the first m AnchorPCA∞ directions span S⋆ (see Thm. C.8). Moreover, § B.2.2 shows that under an additional block-diagonality assumption on Σ, AnchorPCAλ recovers S⋆ already for finite values of λ (this assumption is satisfied in the example of § 1.2).

Optimal reconstruction under perturbations

We now study robustness of the solutions under structured covariance inflations. It is straightforward to see that S⋆ satisfies the following robustness guarantee under eigengap conditions on the covariances Σe (see Prop. B.5 for details): S⋆ is invariant to arbitrarily large admissible domainspecific top-k variance inflations; more precisely, it remains the maximally invariant subspace under (e) perturbations of the form Σ′e = Σe + ∆e with 0 ⪯ ∆e and Im(∆e ) ⊆ Im(Πk ) for all e ∈ E. We now prove that both versions of AnchorPCA satisfy a similar guarantee under bounded versions of the same perturbations. More precisely, for a covariance matrix Σ, the expected reconstruction 5

Avg. reconstruction error

600

120

500

110 100

400

300 200 poolPCA maxRCS

0 = 2E = 150

100 0

100

AnchorPCA = 25 AnchorPCA

maxRegret norm-maxRegret

200 300 Perturbation strength

400 (ρ)

500 (e)

Figure 2: Average reconstruction error along the perturbation path Σe = Σe + ρ Πk for the example in § 1.2. As suggested by Thm. 4.1, poolPCA is best only for small perturbations (0 ≤ ρ ≲ 6.08). Among the wcPCA [15] baselines, maxRCS is best only on the tiny interval 6.08 ≲ ρ ≲ 6.15, after which norm-maxRegret is best until ρ ≲ 9.79. AnchorPCA25 dominates on the broad middle regime 9.79 ≲ ρ ≲ 371.83, and AnchorPCA∞ performs best for strong perturbations ρ ≳ 371.83. error of a rank-k projector ΠW is given by E[∥x − ΠW x∥22 ] = Tr(Σ(Ip − ΠW )). We further define, for any λ > 0, the following set of perturbed covariances n o (e) Cλ := (Σ′1 , . . . , Σ′E ) Σ′e = Σe + ∆e ∈ Sp+ , 0 ⪯ ∆e ⪯ 2Eλ Πk for all e ∈ E and obtain the following result. Theorem 4.1 (Robust reconstruction). For W ∈ Op×k , we define the average reconstruction error PE over {Σe }e∈E as R(W ; Σ1 , . . . , ΣE ) := E1 e=1 Tr(Σe (Ip − ΠW )). Then, (i), Wk⌣,λ ∈ arg min

′

max ′

W ∈Op×k (Σ1 ,...,ΣE )∈Cλ

R(W ; Σ′1 , . . . , Σ′E ) ⇐⇒ Wk⌣,λ solves rank-k AnchorPCAλ .

and, (ii), with B := arg minW f ∈Op×k

(e) 2 f − Πk ∥F e=1 ∥ΠW

Wk⌣,∞ ∈ arg min R(W ; Σ1 , . . . , ΣE )

⇐⇒

W ∈B

Wk⌣,∞ solves rank-k AnchorPCA∞ .

Thus, if we consider target domains that are subject to covariance inflations described by the set Cλ , then AnchorPCAλ minimizes the worst-case reconstruction error over such domains. We now aim to provide intuition for part (ii). The proof of Thm. 4.1 shows that max ′ ′

(Σ1 ,...,ΣE )∈Cλ

R(W ; Σ′1 , . . . , Σ′E ) =

E E X 1 X (e) Tr Σe (Ip − ΠW ) + λ ∥ΠW − Πk ∥2F . E e=1 e=1

(4)

Equation (4) shows that the worst-case reconstruction error is affine in λ, with slope equal to the projector-disagreement penalty. AnchorPCA∞ first minimizes this slope over rank-k projectors and then, among projectors with minimal slope, minimizes the nominal reconstruction error.1 Fig. 2 visualizes parts of Thm. 4.1. It shows that poolPCA is optimal only for small perturbations, whereas AnchorPCAλ and AnchorPCA∞ guard against medium-sized and large perturbations. Further details are provided in § D.2. The experiment in § 6.2 shows that robustness guarantees like the ones shown in Thm. 4.1 may be of practical relevance. 1 A single minimax problem over the unbounded union

λ>0 Cλ is typically degenerate, as can be seen from Equation (4). So the constrained formulation is the natural non-asymptotic statement.

Algorithm 1 Block-stabilized empirical AnchorPCA∞ b Σ, b k, toln Input: Π, b ⌣,∞ , S c⋆ , m c ⌣,∞ , Π b∞ Output: W k k b = Pp ρbj vbj vb⊤ with ρb1 ≥ · · · ≥ ρbp . 1: Compute an eigendecomposition Π j j=1 2: Set a ← 1 and ℓ ← 1. 3: while a ≤ p do 4: Set b ← a. 5: while b < p and ρba − ρbb+1 ≤ toln do 6: Set b ← b + 1. 7: end while dℓ ← span{b 8: Set Iℓ ← {a, . . . , b} and Eig vj : j ∈ Iℓ }. 9: Set a ← b + 1 and ℓ ← ℓ + 1. 10: end while P P d d 11: Let r be the first block index such that ℓ≤r dim(Eigℓ ) ≥ k, and set q ← k − ℓ<r dim(Eigℓ ). b d 12: For all ℓ < r, take an ordered eigenbasis of Σ within Eigℓ . b within Eig dr . 13: Take the top-q eigenvectors of Σ ⌣,∞ c 14: Concatenate these vectors to obtain Wk . c ⌣,∞ (W c ⌣,∞ )⊤ . b ⌣,∞ ← W 15: Set Π k k k c⋆ ← Eig d1 and m d1 ). 16: Set S b ∞ ← dim(Eig

Learning methods for Anchor PCA

We now describe the finite-sample methods for both AnchorPCA Pne versions. Suppose domain e ∈ E provides observations xe,1 , . . . , xe,ne ∈ Rp . Define x̄e := n1e i=1 xe,i and the empirical covariance Pne 1 ⊤ p×k b c b e , and write the (xe,i − x̄e )(xe,i − x̄e ) . Let We ∈ R be a top-k eigenbasis of Σ Σe := ne −1

i=1

(e)

b := W ce W ce⊤ , average covariance empirical versions of the domain-specific projection matrix as Π PE b PEk b (e) b b as Σ := e=1 Σe /E, and average projection matrix as Π := e=1 Πk /E, where we omit the bar compared to their population-level counterparts Σ, Π for readability. 5.1

Finite-sample implementation

cλ := Σ b + 2Eλ Π, b see We define empirical AnchorPCAλ as follows: we compute rank-k PCA of M ,λ ⌣ c cλ . For AnchorPCA∞ , however, a Equation (3), and return an ordered top-k eigenbasis W of M k direct plug-in implementation replacing population quantities by its empirical counterparts would b would not yield a consistent method: for example, due to estimation error, all eigenvalues of Π generally be distinct, creating blocks of equal eigenvalues of size one (see § 2.4). Instead, we group b consecutive Π-eigenspaces whenever eigenvalues remain within toln of the first eigenvalue of that group, where n := mine∈E ne ; for each of these groups of eigenspaces, we compute the direct sum d1 , Eig d2 , . . .; see Algorithm 1 for details. If S⋆ ̸= {0} and toln and denote the resulting spaces by Eig separates the empirical eigenvalue-1 block from the rest of the spectrum (in our experiments we use d1 consistently estimates S⋆ (see § 5.3). Thus, the algorithm toln = 0.5n−0.4 as a default choice), Eig also returns estimates of the invariant subspace and its dimension. 5.2

FindS⋆ : Identifying the invariant subspace S⋆

In some settings, the object of interest is not a full k-dimensional subspace (as estimated by AnchorPCA), but only the maximally invariant subspace S⋆ . We now describe how to estimate an unordered basis of S⋆ directly. An eigenvalue problem. Recall that S⋆ is the eigenspace of Π associated with eigenvalue 1 (§ 2). A natural approach would therefore be to identify eigenvalues equal to 1. While this characterization is exact at the population level, it cannot be applied directly to finite samples. Although bootstrapping projection matrices is feasible [e.g., 5], the null hypothesis λi (Π) = 1 lies on the boundary of the 7

parameter space, and empirical eigenvalues will have value strictly less than 1, making the use of bootstrap-based confidence intervals difficult. Moreover, bootstrapping the null or using a parametric bootstrap is challenging, as the null is composite and depends on an unknown invariant subspace. Testing the dimension of S⋆ . We instead directly test the dimension of S⋆ . If dim(S⋆ ) = m b associated with its m largest eigenvalues. Let were known, S⋆ is estimated by the eigenspace of Π mmin := max(0, p − E(p − k)). This motivates the sequential testing procedure H0,s : dim(S⋆ ) = s

HA,s : dim(S⋆ ) < s,

(5)

for s ∈ {k, k − 1, . . . , mmin + 1}, stopping at the first s for which the null is not rejected; if all nulls are rejected, estimate dim(S⋆ ) as mmin . To construct such tests, we adapt ideas from partial common principal component analysis (PCPCA) subspaces [37]. In the PCPCA setting, it is assumed that for all e ∈ E, we observe a sample of i.i.d. Gaussian observations with population positive-definite covariance Σe . Additionally, Σe is assumed to satisfy a strict eigengap between its k th and (k + 1)st eigenvalues. Under these assumptions, Wald-type tests are available for the dimension of the union of the top-k eigenspaces of (Σe )e∈E . To translate this to our setting, we consider the orthogonal (e) complement projector matrices Ip − Πk and let r be the dimension of the union of their eigenspaces. (e) Since a vector belongs to S⋆ if and only if it lies in the null space of Ip − Πk for every e ∈ E, r = p − dim(S⋆ ). We apply the PCPCA subspace test to the bottom-(p − k) eigenspaces of the (e) covariance matrices, equivalently to Im(Ip − Πk ), and thus obtain a sequential testing procedure for (5). Combining this sequential test with the eigenspace estimator above, we define FindS⋆ as b associated with its the procedure that first estimates dim(S⋆ ) and then returns the eigenspace of Π estimated number of largest eigenvalues. The test statistic and further details are given in § B.3. 5.3

Consistency of Anchor PCA

Both AnchorPCA methods are consistent under mild conditions (the full statement and proof are given in § C.4). Assume that the minimum domain sample size n := mine∈E ne → ∞, that observations are centered with finite second moments, and that each domain covariance Σe has a strict eigengap at b → Σ and Π b → Π almost surely, and AnchorPCAλ is consistent whenever Mλ itself rank k. Then Σ has a strict eigengap at rank k. For AnchorPCA∞ , consistency requires a little more care: the block tolerance used in the algorithm must shrink to zero, but slowly enough relative to the estimation error b Any tolerance of the form toln = cn−α with 0 < α < 1/2 works under finite fourth moments. of Π. Under these conditions, the full rank-k projector is consistent when the population AnchorPCA∞ solution projector is unique; moreover, if S⋆ ̸= {0}, the first block returned by AnchorPCA∞ is a consistent estimator of S⋆ . For FindS⋆ , Schott’s asymptotic theory [37] implies that the probability of correctly estimating dim(S⋆ ) converges to 1 − α, where α is the nominal level of the tests.

Experiments

Code to reproduce all simulations and figures is publicly available at https://github.com/ benesei/anchor-pca. 6.1

Recovering the invariant subspace S⋆

In this experiment, all environments have equal sample size, which we call N (which thus also d1 , of the block-stabilized equals n = mine∈E ne ). § 5.3 shows that the first grouped eigenspace, Eig AnchorPCA∞ construction consistently estimates S⋆ when the tolerance shrinks while dominating the empirical perturbation. We test this statement in a random-subspace model where S⋆ is known, using the default automatic block tolerance (tolN = 0.5N −0.4 ). Below we choose the value m = dim(S⋆ ) c⋆ . We measure subspace error by ∥Π b − ΠS ∥op . and let m b ∞ denote the dimension of the estimate S ⋆ S⋆ The distribution draws are as follows. We draw a random m-dimensional subspace S⋆ ⊂ Rp and, for each domain, a different rank-k principal subspace containing S⋆ . The covariance in that domain is then constructed so that this subspace is exactly its top-k eigenspace. We consider two DGP regimes. In the easy setting, stable and domain-specific top eigenvalues have comparable sizes; in the difficult setting, the stable directions are weaker and the domain-specific top directions are stronger. § D.3 8

Dimension recovery

1.0 c = m) Pr(m

0.6 0.4

0.8

AnchorPCA ∞

0.6

FindS

color = DGP easy hard

0.4

0.2 0.0

texture = method

kΠ S − Π S k op

0.8

Subspace error

1.0

0.2

102

103

0.0

102

N per environment

103

N per environment

Figure 3: Recovery of S⋆ by AnchorPCA∞ and FindS⋆ (§ 5.2). Gaussian data with E = 5, p = 10, k = 5, and m = dim(S⋆ ) = 2. We show mean correct-dimension probability (left) and median operator-norm projector error (right) against the true S⋆ . As predicted by § 5.3, the first grouped eigenspace of AnchorPCA∞ (solid) recovers S⋆ with growing sample size and FindS⋆ (dashed) has dimension-recovery error close to the nominal test size of 0.05. For each sample size N , the curves show the average over 100 independent distribution draws, for each of which we draw 20 independent datasets. Shaded bands show 10–90% quantiles across distribution draws (shown only for AnchorPCA∞ ).

provides details and two even more challenging configurations: a smaller invariant subspace (m = 1) and a low-domain regime with 2 = E < k = 5, both of which have smaller agreement-separation gaps ρm (Π) − ρm+1 (Π) than the main configuration. Fig. 3 is in agreement with the consistent recovery of S⋆ . In the easy setting, the correct block dimension is already recovered in about 64% of runs at N = 50 and reaches probability 1 by N = 5000; in the hard setting, the corresponding probability is about 24% at N = 50 but also reaches 1 by N = 5000. The median projector error decreases to about 0.02 in the easy setting and 0.04 in the difficult setting at N = 5000. § 5.3 does not assume Gaussianity. Indeed, § D.3.4 shows the same qualitative behavior for Gaussian-mixture samples. The dashed curves in Fig. 3 show that FindS⋆ gives a comparable subspace diagnostic in this configuration, while its dimension-recovery error remains close to 0.05 for large sample sizes, the nominal level of the individual tests, consistent with § 5.3.

6.2

Gas sensor array drift

We evaluate AnchorPCA on gas-sensor drift data [13, 41], consisting of 13,910 recordings from 16 chemical sensors, summarized by p = 128 features and grouped into ten temporal batches. We treat batches as domains: source batches are used for fitting and later batches are held out. The data are known to exhibit sensor drift, i.e., gradual changes in sensor response caused by aging, poisoning, or environmental effects, so that the same gas can induce different response distributions over time [13, 42]. The hope is that the projector agreement of AnchorPCA across early batches downweights drift-specific high-variance directions and retains directions that remain informative under future b b W )/ Tr(Σ b b ), temporal shifts. We evaluate by normalized batch explained variance 100 · Tr(W ⊤ Σ where W is fit only on source batches. Fig. 4 shows the results for the split B1–B6 → B7–B10 for k = 20. The finite-penalty method AnchorPCAλ=1 illustrates the trade-off: its mean source explained variance of 95.7% lies between those of poolPCA (98.1%) and AnchorPCA∞ (92.8%); but its mean target explained variance (80.6%) outperforms poolPCA (75.5%). AnchorPCA∞ has the weakest source fit (mainly due to B2) but yields an even higher target explained variance (87.3%) and improves over poolPCA on every heldout target batch. In § D.4, we consider alternative splits and values of k: in only two cases (eight source domains and k = 20, 30) poolPCA seems to have sufficient information to perform well on one of the target domains; in all other settings, AnchorPCA outperforms poolPCA. 9

% explained variance

100 Method AnchorPCA ∞ AnchorPCA λ = 1 norm-maxRegret poolPCA

80 60

90.2%

Source batches

Target batches

B10

Figure 4: Gas-sensor drift, source B1–B6 and target B7–B10. All methods use k = 20 and sourceonly preprocessing/fitting; norm-maxRegret is the normalized regret wcPCA baseline [15]. The finite-penalty method AnchorPCAλ=1 exhibits the source–target compromise: it retains more source explained variance than AnchorPCA∞ but yields smaller average target gains. AnchorPCA∞ gives up more source explained variance, achieves the highest mean target EV, and improves over poolPCA on every held-out target batch; the arrow marks its largest relative target gain over poolPCA (B9, 90.2%).

Discussion and future work

We introduced Anchor PCA, a novel framework for PCA on multiple domains. We prove recovery guarantees and robustness results and show its effectiveness on simulated and real data. Our invariance notion (membership in the local top-k principal subspace of every domain) depends on the choice of k and can be sensitive to near eigenvalue ties. The minimax guarantee in Thm. 4.1 covers bounded variance inflations inside local top-k eigenspaces, but not mean shifts, variance deflations, or rotations of the principal subspaces. Finally, λ needs to be chosen in practice but we believe that λ = ∞ is a reasonable choice. While we have focused on robust linear dimension reduction based on PCA, we discuss possible nonlinear extensions of AnchorPCAλ based on autoencoders in § E.

Acknowledgments The authors used GPT-5.4 (xhigh) and Claude Opus 4.6 for assistance with manuscript drafting, proofreading of human-written text and mathematical proofs, and development of simulations, figures, and diagrams, but take full responsibility for the content. JvK is supported by The Branco Weiss Fellowship—Society in Science.

References [1] Abubakar Abid, Martin J. Zhang, Vivek K. Bagaria, and James Zou. Exploring patterns enriched in a dataset with contrastive principal component analysis. Nature Communications, 9(1):2134, 2018. doi: 10.1038/s41467-018-04608-8. 15 [2] Alekh Agarwal and Tong Zhang. Minimax regret optimization for robust machine learning under distribution shift. In Proceedings of the 35th Conference on Learning Theory, volume 178 of Proceedings of Machine Learning Research, pages 2704–2729. PMLR, 2022. 40 [3] Martin Arjovsky, Léon Bottou, Ishaan Gulrajani, and David Lopez-Paz. Invariant risk minimization. arXiv preprint arXiv:1907.02893, 2019. 1, 15 [4] Krishna B. Athreya and Soumendra N. Lahiri. Measure Theory and Probability Theory. Springer Texts in Statistics. Springer, New York, 2006. doi: 10.1007/978-0-387-35434-7. 23 [5] Rudolf Beran and Muni S. Srivastava. Bootstrap tests and confidence regions for functions of a covariance matrix. The Annals of Statistics, 13(1):95–115, 1985. doi: 10.1214/aos/1176346579. 7 [6] Rajendra Bhatia. Matrix Analysis, volume 169 of Graduate Texts in Mathematics. Springer, New York, 1997. doi: 10.1007/978-1-4612-0653-8. 16 10

[7] Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, 2004. doi: 10.1017/CBO9780511804441. 28 [8] Rune Christiansen, Niklas Pfister, Martin Emil Jakobsen, Nicola Gnecco, and Jonas Peters. A causal framework for distribution generalization. IEEE Transactions on Pattern Analysis and Machine Intelligence, 44(10):6614–6630, 2022. doi: 10.1109/TPAMI.2021.3094760. 15 [9] Petros Drineas and Ilse C. F. Ipsen. Low-rank matrix approximations do not need a singular value gap. SIAM Journal on Matrix Analysis and Applications, 40(1):299–319, 2019. doi: 10.1137/18M1163658. 22 [10] Qing Feng, Meilei Jiang, Jan Hannig, and J. S. Marron. Angle-based joint and individual variation explained. Journal of Multivariate Analysis, 166:241–265, 2018. doi: 10.1016/j.jmva. 2018.03.008. 15 [11] Bernhard N. Flury. Common principal components in k groups. Journal of the American Statistical Association, 79(388):892–898, 1984. doi: 10.1080/01621459.1984.10477108. 15 [12] Bernhard N. Flury and Walter Gautschi. An algorithm for simultaneous orthogonal transformation of several positive definite symmetric matrices to nearly diagonal form. SIAM Journal on Scientific and Statistical Computing, 7(1):169–184, 1986. doi: 10.1137/0907013. 15 [13] Jordi Fonollosa, Irene Rodríguez-Luján, and Ramón Huerta. Chemical gas sensor array dataset. Data in Brief, 3:85–89, 2015. doi: 10.1016/j.dib.2015.01.003. 9, 36 [14] Francesco Freni, Anya Fries, Linus Kühne, Markus Reichstein, and Jonas Peters. Maximum risk minimization with random forests. arXiv preprint arXiv:2512.10445, 2025. 40 [15] Anya Fries, Markus Reichstein, David Blei, and Jonas Peters. Worst-case low-rank approximations. arXiv preprint arXiv:2603.11304, 2026. 1, 6, 10, 15, 32, 37 [16] Tilmann Gneiting and Adrian E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal of the American Statistical Association, 102(477):359–378, 2007. doi: 10.1198/ 016214506000001437. 40 [17] Trygve Haavelmo. The statistical implications of a system of simultaneous equations. Econometrica, 11(1):1–12, 1943. doi: 10.2307/1905714. 15 [18] Geoffrey E. Hinton and Ruslan R. Salakhutdinov. Reducing the dimensionality of data with neural networks. Science, 313(5786):504–507, 2006. doi: 10.1126/science.1127647. 39 [19] Roger A. Horn and Charles R. Johnson. Matrix Analysis. Cambridge University Press, 2 edition, 2012. doi: 10.1017/CBO9780511810817. 19, 22, 24, 25, 26, 30 [20] Ian T. Jolliffe. Principal Component Analysis. Springer Series in Statistics. Springer, New York, 2 edition, 2002. doi: 10.1007/b98835. 1 [21] Margherita Lazzaretto, Jonas Peters, and Niklas Pfister. Invariant subspace decomposition. Journal of Machine Learning Research, 26(95):1–56, 2025. 15 [22] Eric F. Lock, Katherine A. Hoadley, J. S. Marron, and Andrew B. Nobel. Joint and individual variation explained (JIVE) for integrated analysis of multiple data types. The Annals of Applied Statistics, 7(1):523–542, 2013. doi: 10.1214/12-AOAS597. 15 [23] Chaochao Lu, Yuhuai Wu, José Miguel Hernández-Lobato, and Bernhard Schölkopf. Invariant causal representation learning for out-of-distribution generalization. In International Conference on Learning Representations, 2022. 15 [24] Mirco Migliavacca, Talie Musavi, Miguel D. Mahecha, Jacob A. Nelson, Jürgen Knauer, Dennis D. Baldocchi, Oscar Perez-Priego, Rune Christiansen, Jonas Peters, et al. The three major axes of terrestrial ecosystem function. Nature, 598(7881):468–472, 2021. doi: 10.1038/ s41586-021-03939-9. 1 11

[25] Krikamol Muandet, David Balduzzi, and Bernhard Schölkopf. Domain generalization via invariant feature representation. In Sanjoy Dasgupta and David McAllester, editors, Proceedings of the 30th International Conference on Machine Learning, volume 28 of Proceedings of Machine Learning Research, pages 10–18. PMLR, 2013. 15 [26] Yotam Norman and Ron Meir. Unsupervised representation learning – an invariant risk minimization perspective. In International Conference on Learning Representations, 2026. 1, 15 [27] John Novembre, Toby Johnson, Katarzyna Bryc, Zoltán Kutalik, Adam R. Boyko, Adam Auton, Amit Indap, Karen S. King, Sven Bergmann, Matthew R. Nelson, Matthew Stephens, and Carlos D. Bustamante. Genes mirror geography within europe. Nature, 456(7218):98–101, 2008. doi: 10.1038/nature07331. 1 [28] Sinno Jialin Pan, Ivor W. Tsang, James T. Kwok, and Qiang Yang. Domain adaptation via transfer component analysis. IEEE Transactions on Neural Networks, 22(2):199–210, 2011. doi: 10.1109/TNN.2010.2091281. 15 [29] Judea Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, 2 edition, 2009. doi: 10.1017/CBO9780511803161. 15 [30] Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen. Causal inference by using invariant prediction: Identification and confidence intervals. Journal of the Royal Statistical Society: Series B, 78(5):947–1012, 2016. doi: 10.1111/rssb.12167. 15 [31] Herbert Robbins and Sutton Monro. A stochastic approximation method. The Annals of Mathematical Statistics, 22(3):400–407, 1951. doi: 10.1214/aoms/1177729586. 39 [32] Mateo Rojas-Carulla, Bernhard Schölkopf, Richard Turner, and Jonas Peters. Invariant models for causal transfer learning. Journal of Machine Learning Research, 19(36):1–34, 2018. 15 [33] Dominik Rothenhäusler, Nicolai Meinshausen, Peter Bühlmann, and Jonas Peters. Anchor regression: Heterogeneous data meet causality. Journal of the Royal Statistical Society: Series B, 83(2):215–246, 2021. doi: 10.1111/rssb.12398. 1, 3, 15, 32 [34] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams. Learning representations by back-propagating errors. Nature, 323(6088):533–536, 1986. doi: 10.1038/323533a0. 39 [35] Samira Samadi, Uthaipon Tantipongpipat, Jamie Morgenstern, Mohit Singh, and Santosh Vempala. The price of fair PCA: One extra dimension. In Advances in Neural Information Processing Systems, volume 31, 2018. 14 [36] Bernhard Schölkopf, Dominik Janzing, Jonas Peters, Eleni Sgouritsa, Kun Zhang, and Joris M. Mooij. On causal and anticausal learning. In Proceedings of the 29th International Conference on Machine Learning, pages 1255–1262. Omnipress, 2012. 15 [37] James R. Schott. Partial common principal component subspaces. Biometrika, 86(4):899–908, 1999. doi: 10.1093/biomet/86.4.899. 8, 15, 20 [38] Xinwei Shen and Nicolai Meinshausen. Distributional principal autoencoders. arXiv preprint arXiv:2404.13649, 2024. 15, 40 [39] Uthaipon Tantipongpipat, Samira Samadi, Mohit Singh, Jamie Morgenstern, and Santosh Vempala. Multi-criteria dimensionality reduction with applications to fairness. In Advances in Neural Information Processing Systems, volume 32, 2019. 15 [40] Matthew Turk and Alex Pentland. Eigenfaces for recognition. Journal of Cognitive Neuroscience, 3(1):71–86, 1991. doi: 10.1162/jocn.1991.3.1.71. 1 [41] Alexander Vergara. Gas sensor array drift at different concentrations. UCI Machine Learning Repository, 2012. Dataset. 9, 36 [42] Alexander Vergara, Shankar Vembu, Tuba Ayhan, Margie A. Ryan, Margie L. Homer, and Ramón Huerta. Chemical gas sensor drift compensation using classifier ensembles. Sensors and Actuators B: Chemical, 166–167:320–329, 2012. doi: 10.1016/j.snb.2012.01.074. 9 12

[43] Bingkai Wang, Xi Luo, Yi Zhao, and Brian Caffo. Semiparametric partial common principal component analysis for covariance matrices. Biometrics, 77(4):1175–1186, 2021. doi: 10.1111/ biom.13369. 15 [44] Haoxiang Wang, Haozhe Si, Bo Li, and Han Zhao. Provable domain generalization via invariantfeature subspace recovery. In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 23018– 23033. PMLR, 17–23 Jul 2022. 1, 15 [45] Zhenyu Wang, Molei Liu, Jing Lei, Francis Bach, and Zijian Guo. StablePCA: Distributionally robust learning of shared representations from multi-source data. arXiv preprint arXiv:2505.00940, 2025. 1, 15 [46] Ian J. Wright, Peter B. Reich, Mark Westoby, David D. Ackerly, Zdravko Baruch, et al. The worldwide leaf economics spectrum. Nature, 428(6985):821–827, 2004. doi: 10.1038/ nature02403. 1 [47] Yi Yu, Tengyao Wang, and Richard J. Samworth. A useful variant of the Davis–Kahan theorem for statisticians. Biometrika, 102(2):315–323, 2015. doi: 10.1093/biomet/asv008. 22

Contents of the Appendix A Further related work

A.1 Multi-domain dimension reduction . . . . . . . . . . . . . . . . . . . . . . . . . .

A.2 Invariance for prediction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

B Further details on Anchor PCA

B.1 Ordered bases for direction-level statements . . . . . . . . . . . . . . . . . . . . .

B.2 Additional theoretical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

B.3 Test statistic of FindS⋆ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C Proofs

C.1 Additional intermediate results . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.2 Proof of Thm. 3.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.3 Proof of Thm. 4.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.4 Consistency of Anchor PCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

D Details on experiments

D.1 Details of the motivating example . . . . . . . . . . . . . . . . . . . . . . . . . .

D.2 Details on the perturbation-path illustration . . . . . . . . . . . . . . . . . . . . .

D.3 Details on the simulations

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

D.4 Gas sensor array drift details . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

E Nonlinear extensions

E.1 Vanilla autoencoder (AE) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

E.2 Anchor autoencoder (AnchorAE) . . . . . . . . . . . . . . . . . . . . . . . . . . .

E.3 Distributional principal autoencoder (DPA) . . . . . . . . . . . . . . . . . . . . .

E.4 Multi-domain DPA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Further related work

AnchorPCA is a dimensionality reduction over several source domains that builds on ideas of invariance. We summarize work in multi-domain dimension reduction and related ideas in prediction that use invariance. A.1

Multi-domain dimension reduction

Multi-domain dimension reduction can be separated into three groups: worst-case methods, invariance-based methods, and other methods (e.g., those which do not relate to generalization, but consider multiple groups and PCA). Worst-case methods. Worst-case methods replace pooled explained variance or reconstruction by a max-min (or min-max) objective, which seeks to optimize performance on the worst domain, group, or criterion. FairPCA [35] minimizes the maximum increase in reconstruction error incurred by using a shared projection instead of the domain-specific optimal projection. Multi-criteria dimensionality 14

reduction [39] generalizes this idea to several criteria over multiple groups. However, neither provides a robustness perspective. In contrast, wcPCA [15] includes out-of-sample guarantees for target covariances in the convex hull of the observed source covariances. Along similar lines, StablePCA [45] casts multi-source PCA as a group-DRO problem. AnchorPCAλ is also worst-case optimal, but over a different uncertainty set of perturbed covariances (Thm. 4.1). Importantly, it does not only optimize a worst-case objective: it also recovers the maximal invariant subspace S⋆ (Thm. 3.1), a target none of the methods above consider. Invariance-based methods. Invariance-based dimension reduction methods seek representations whose distributions, moments, or predictive relations are stable across domains, often with the goal of improving downstream prediction for a target Y . Transfer Component Analysis (TCA) [28] aligns source and target distributions through MMD-type criteria while preserving variance; DomainInvariant Component Analysis (DICA) [25] learns transformations that reduce dissimilarity across training domains while preserving predictive structure. In contrast, Principal Invariant Component Analysis (PICA) [26] does not require information about Y and finds linear components whose projections have the same distribution across domains. AnchorPCAλ does not require information about a target Y either and uses a geometric notion of invariance: agreement between the domainspecific top-k principal projectors. This is weaker than what PICA aims for (since directions can be invariant even when their variances or ranks differ); We trade off invariance with explained variance, which ensures our solution is always non-empty. While DICA provides generalization guarantees for downstream prediction, TCA, DICA, and PICA do not provide unsupervised PCA-style robustness guarantees for explained variance or reconstruction of the kind studied here. Other methods. Other multi-group dimensionality reduction methods model shared structure for statistical or exploratory purposes rather than for OOD robustness. Common principal component analysis assumes simultaneous diagonalizability of the covariance matrices [11, 12]; partial common principal component methods estimate or test for shared eigendirections or shared principal subspaces [37, 43]. Our FindS⋆ procedure builds on this testing tradition, but the population target is different: S⋆ is the intersection of local top-k eigenspaces, not a set of eigenvectors that are shared as individual axes or ordered in the same way across domains. Multi-block methods such as JIVE [22] and AJIVE [10] decompose data sources into joint and individual low-rank variation, whereas contrastive PCA [1] searches for directions enriched in one dataset relative to a background dataset. These methods provide useful decompositions, but they do not couple shared-structure estimation to the projector-agreement and structured-perturbation guarantees studied here. Nonlinear PCA and autoencoder variants provide another axis of generalization [38]; § E discusses how the same agreement principle could be ported to autoencoder-based representations. A.2

Invariance for prediction

In supervised domain generalization, invariance is usually formulated as stability of a conditional law, a prediction rule, or a risk across domains. This principle is rooted in causal invariance [17, 29, 30, 36] and appears in invariant models for causal transfer learning, invariant risk minimization, invariant-feature subspace recovery, invariant causal representation learning, and invariant subspace decomposition for time-varying prediction problems [3, 21, 23, 32, 44]. Several prediction methods relax strict invariance rather than enforcing it exactly: for example, anchor regression interpolates between ordinary least squares and an invariant solution through a penalty parameter [8, 33]. AnchorPCAλ imports this soft-invariance logic into unsupervised dimension reduction: finite λ penalizes projector disagreement, while AnchorPCA∞ recovers the hard agreement target and, when it exists, the maximal invariant subspace S⋆ .

Further details on Anchor PCA

B.1

Ordered bases for direction-level statements

The optimization problems defining AnchorPCAλ (Defn. 2.3) and AnchorPCA∞ (Defn. 2.2) depend on W only through the projection matrix ΠW = W W ⊤ . Hence, the canonical population objects are solution subspaces or solution projectors. However, the projection-based formulation allows us to fix an ordered representation, see § 2.3. We now state their precise construction. 15

B.1.1

AnchorPCAλ

Solutions of AnchorPCAλ can be obtained by solving rank-k PCA on Mλ = Σ + 2Eλ Π (§ 2.3). Hence, we use the standard PCA convention and choose Wk⌣,λ = (w1⌣,λ , . . . , wk⌣,λ ) as an ordered orthonormal eigenbasis of Mλ associated with its k largest eigenvalues. Thus, w1⌣,λ maximizes w⊤ Mλ w over unit vectors, w2⌣,λ maximizes the same quantity subject to being orthogonal to w1⌣,λ , and so on. If Mλ has eigenvalue ties, the ordered basis is not unique and we choose any orthonormal basis inside tied eigenspaces. If all ties occur within the selected top-k spectrum, this affects only the solution basis, not the rank-k solution subspace and projector. If there is no eigengap at rank-k, i.e. λk (Mλ ) = λk+1 (Mλ ), then the rank-k solution projector is non-unique, and any maximizing rank-k solution projector may be used. B.1.2

AnchorPCA∞

We now define an analogous ordered basis for AnchorPCA∞ in a two-stage construction. First, we determine the solution subspace: directions are ordered by agreement with the domain-specific top-k principal subspaces (equivalently, by decreasing eigenvalue of Π). Second, within each eigenspace, we choose a basis by iteratively maximizing the pooled explained variance over directions in that block. Step 1a: Ordering by agreement. The agreement constraint of AnchorPCA∞ in (1) is equivalent to maximizing k X wi⊤ Πwi Tr(W ⊤ ΠW ) = (6) i=1

over W = [w1 , . . . , wk ] ∈ Op×k (§ 2.3). We now characterize this solution. Let ρ1 > ρ2 > · · · > ρL be the distinct eigenvalues of Π, and let Eigℓ := ker(Π − ρℓ Ip ),

dℓ := dim(Eigℓ ).

Eigℓ is the eigenspace of Π corresponding to the eigenvalue ρℓ and dℓ is its dimension. Larger ρℓ implies larger average agreement with the domain-specific top-k principal subspaces. Then, let r be the unique index satisfying X X dℓ < k ≤ dℓ . ℓ<r

ℓ≤r

By Ky Fan’s maximum principle [6, Problem I.6.15], the trace in (6) is maximized precisely when the columns of W are drawn from the leading eigenspaces of Π: every maximizer spans a k-dimensional subspace of the form Eig1 ⊕ · · · ⊕ Eigr−1 ⊕ Ur , P where Ur ≤ Eigr is an arbitrary subspace of dimension q := k − ℓ<r dℓ . Step 1b: Determining Ur . The eigenspaces Eig1 , . . . , Eigr−1 are selected completely, and only the final eigenspace Ur can be determined via the secondary Σ-criterion of (2.2), which chooses Ur to maximize pooled explained variance among all q-dimensional subspaces of Eigr . Equivalently, inside Eigr , the subspace Ur is chosen by the top-q eigenspace of ΠEigr ΣΠEigr . Step 2: Choosing an ordered basis. Within each fully selected eigenspace Eig1 , . . . , Eigr−1 , we choose an orthonormal eigenbasis ordered by non-increasing eigenvalue of Σ, and concatenate these bases in decreasing order of ρℓ . We then append an ordered top-q eigenbasis of Σ within Eigr . With this convention, AnchorPCA∞ directions are ordered first by invariance, that is, by decreasing Π-eigenvalue. Inside a fixed eigenspace, they are ordered by pooled explained variance. Remark B.1 (Handling non-uniqueness.). If Σ has eigenvalue ties in one of the fully selected eigenspaces Eig1 , . . . , Eigr−1 , the basis within the tied eigenspace is not unique; we make an arbitrary choice, and this affects only direction-level quantities, not the solution subspace or projector. If instead ΠEigr ΣΠEigr has ties among its top q eigenvalues, the subspace Ur itself is not unique, and the rank-k AnchorPCA∞ solution subspace and projector are likewise non-unique; the uniqueness statements elsewhere in the paper (e.g., Thm. C.9(ii)) exclude this case by assumption. 16

B.2

Additional theoretical results

Notation. Throughout the appendix, all matrices are represented with respect to the standard basis of Rp , unless otherwise stated. For pairwise orthogonal subspaces S1 , . . . , Sj ⊆ Rp , we write S1 ⊕ · · · ⊕ Sj for their orthogonal direct sum. Suppose that a symmetric matrix M ∈ Rp×p preserves each Si , that is, M v ∈ Si whenever v ∈ Si . Then we can block-diagonalize M with respect to (Si )ji=1 as follows. Let Bi ∈ Rp×di be an orthonormal basis matrix for Si , where di := dim(Si ), and let Q := [B1 · · · Bj ]. In this aligned basis, M is block diagonal, that is   A1 0   .. Q⊤ M Q = diag(A1 , . . . , Aj ) =  Ai := Bi⊤ M Bi . , . 0 Aj We abbreviate this block-diagonal representation by writing M = A1 ⊕ · · · ⊕ Aj with respect to the decomposition S1 ⊕ · · · ⊕ Sj . B.2.1

AnchorPCA∞ as a limit of AnchorPCAλ

Intuitively, as the penalty λ grows, solutions of AnchorPCAλ approach the hard-agreement solutions of AnchorPCA∞ . The following proposition states this at the level of projection matrices. Thm. C.8 (ii) strengthens this at the level of individual directions rather than subspaces and relates the limit to the invariant subspace S⋆ . ,λ Proposition B.2 (Large-penalty limit). Let Π⌣ be any rank-k AnchorPCAλ solution projector. If k ,λ ⌣ n λn → ∞, every accumulation point of (Πk ) is a rank-k AnchorPCA∞ solution projector. ,∞ If, in addition, the population rank-k AnchorPCA∞ solution projector is unique, denote it by Π⌣ , k then ,λ ,∞ Π⌣ →0 as λ → ∞. − Π⌣ k k op

Proof. Let Pk := {X ∈ Rp×p : X = W W ⊤ , W ∈ Op×k } denote the set of rank-k orthogonal projectors. Since Pk is closed and bounded in Rp×p , it is compact. For X ∈ Pk , define H(X) := Tr(ΣX),

G(X) := Tr(ΠX),

Fλ (X) := H(X) + 2Eλ G(X).

Since both AnchorPCA objectives depend on W only through W W ⊤ , rank-k AnchorPCAλ is the problem of maximizing Fλ over Pk , while rank-k AnchorPCA∞ is the problem of first maximizing G over Pk and then maximizing H among the maximizers of G. ,λn Fix a sequence λn → ∞, and for each n, let Π⌣ ∈ Pk be a rank-k AnchorPCAλ solution projector k ,λn at penalty λn . By compactness of Pk , the sequence (Π⌣ ) has accumulation points in Pk . Let Π⋆ k ,λn be any such accumulation point; passing to a subsequence (and relabeling), Π⌣ − Π⋆ k

→ 0. op

Step 1: Π⋆ maximizes G over Pk . By contradiction, assume that Π⋆ does not maximize G and that there exist Y ∈ Pk and δ > 0 such that G(Y ) ≥ G(Π⋆ ) + 2δ. By continuity of G, ,λn G(Π⌣ ) ≤ G(Π⋆ ) + δ k

for all sufficiently large n. Further, for all X ∈ Pk , |H(X)| = | Tr(ΣX)| ≤ k Σ op (because X ,λn ) ≥ −2k Σ op . Therefore, for all is a rank-k orthogonal projector), and hence H(Y ) − H(Π⌣ k sufficiently large n, ,λn ,λn ,λn Fλn (Y ) − Fλn (Π⌣ ) = H(Y ) − H(Π⌣ ) + 2Eλn G(Y ) − G(Π⌣ ) k k k

≥ −2k Σ op + 2Eλn δ. 17

,λn ) ≥ G(Π⋆ ) + 2δ − (G(Π⋆ ) + δ) = 2δ − δ = δ. Since λn → ∞, the where G(Y ) − G(Π⌣ k ,λn right-hand side eventually is positive, contradicting the optimality of Π⌣ . Thus, Π⋆ maximizes G k over Pk .

Step 2: Π⋆ maximizes H among the maximizers of G. By contradiction, assume that there exists Y ∈ Pk and δ > 0 such that Y maximizes G and H(Y ) ≥ H(Π⋆ ) + 2δ. Since Y maximizes G over ,λn Pk , G(Y ) ≥ G(Π⌣ ) for all n. Further, by continuity of H, k ,λn H(Π⌣ ) ≤ H(Π⋆ ) + δ k

for all sufficiently large n. Hence, for all sufficiently large n, ,λn ,λn ,λn Fλn (Y ) − Fλn (Π⌣ ) = H(Y ) − H(Π⌣ ) + 2Eλn G(Y ) − G(Π⌣ ) ≥ δ, k k k ,λn which contradicts the optimality of Π⌣ . Therefore, Π⋆ is a rank-k AnchorPCA∞ solution projeck tor.

Step 3: full convergence under uniqueness. Assume now that the rank-k AnchorPCA∞ solution ,∞ projector is unique and equals Π⌣ . By contradiction, assume that there exist ε > 0 and a sequence k λn → ∞ such that ,λn ,∞ Π⌣ − Π⌣ ≥ε for all n. k k op

,λn ) has an accumulation point Π⋆ . By the first part of the By compactness of Pk , the sequence (Π⌣ k ,∞ proposition, Π⋆ is a rank-k AnchorPCA∞ solution projector, hence Π⋆ = Π⌣ , a contradiction. k

Therefore

,λ ,∞ − Π⌣ Π⌣ k k

B.2.2

→ 0. op

Recovery of S⋆ by AnchorPCAλ with finite penalty λ

The following proposition shows that AnchorPCAλ recovers the invariant subspace at two levels: if λ > λ⋆ , every rank-k solution subspace contains S⋆ , and if λ > λ⋆⋆ , the leading rank-m AnchorPCAλ directions (in the convention of § B.1) span S⋆ . Proposition B.3 (Finite-penalty recovery). Assume 1 ≤ m ≤ k < p and that Σ is block-diagonal with respect to S⋆ ⊕ S⋆⊥ . Let ΣS⋆ := ΠS⋆ ΣΠS⋆ and ΣS⋆⊥ := ΠS⋆⊥ ΣΠS⋆⊥ be the diagonal blocks ,λ of Σ with respect to S⋆ ⊕ S⋆⊥ and let R := ΠS⋆⊥ ΠΠS⋆⊥ . If Π⌣ is a rank-k AnchorPCAλ solution k projector, then the following statements hold. ,λ (i) If λ > λ⋆ := λ1 (ΣS⋆⊥ ) − λm (ΣS⋆ ) /2E 1 − λk−m+1 (R) , then Π⌣ ΠS⋆ = ΠS⋆ , equivak ⌣,λ lently S⋆ ⊆ Sk . (ii) If λ > λ⋆⋆ := λ1 (ΣS⋆⊥ ) − λm (ΣS⋆ ) /2E 1 − λ1 (R) , then the top m-dimensional eigenspace of Mλ = Σ + 2EλΠ is S⋆ . Proof. Since rank-k AnchorPCAλ is rank-k PCA on Mλ = Σ + 2Eλ Π, the proof proceeds by block-diagonalizing Mλ along S⋆ ⊕ S⋆⊥ and then showing that, once λ is large enough, the S⋆ -block eigenvalues dominate the relevant S⋆⊥ -block eigenvalues. Step 1: block-diagonalize Mλ along S⋆ ⊕S⋆⊥ . Let S ∈ Op×m and Z ∈ Op×(p−m) have orthonormal columns respectively spanning S⋆ and S⋆⊥ . Let A := S ⊤ ΣS,

B := Z ⊤ ΣZ,

R := Z ⊤ ΠZ

Then, in the orthonormal coordinates (S, Z), (S, Z)⊤ Mλ (S, Z) =

A + 2EλIm 0 18

0 B + 2EλR

because S ⊤ ΣZ = 0 (because, by assumption, Σ is block-diagonal with respect to S⋆ ⊕ S⋆⊥ ) and S ⊤ ΠS = Im and S ⊤ ΠZ = 0 (as S⋆ is the eigenspace of Π associated with eigenvalue 1). Step 2: compare the top eigenvalues across the two blocks. First, observe that ΣS⋆ = SAS ⊤ ,

ΣS⋆⊥ = ZBZ ⊤ ,

R = ZRZ ⊤ .

Hence the eigenvalues of A, B and R that we will require are as follows λm (A) = λm (ΣS⋆ ),

λ1 (B) = λ1 (ΣS⋆⊥ ),

λi (R) = λi (R), i ∈ {1, . . . , p − m}.

Furthermore, λ1 (R) < 1: to prove this, assume, by contradiction that Rv = v for some nonzero (e) v ∈ Rp−m . Then, as R = Z ⊤ ΠZ, z := Zv ∈ S⋆⊥ satisfies Πz = z, which implies Πk z = z for every e, hence z ∈ S⋆ , a contradiction. The smallest eigenvalue of the invariant block and the eigenvalues of the non-invariant block are thus am (λ) := λm (A) + 2Eλ,

bj (λ) := λj (B + 2EλR) for j ∈ {1, . . . , p − m}.

Step 3: show inclusion of S⋆ for λ > λ⋆ (part (i)). By Weyl’s inequality [19, Theorem 4.3.7] used in decreasing-order notation, applied to B, 2EλR with indices (1, k − m + 1) (recall that we are using decreasing order notation of eigenvalues, whereas Thm. 4.3.7 in [19] uses increasing order), bk−m+1 (λ) ≤ λ1 (B) + 2Eλ λk−m+1 (R). λ1 (B)−λm (A) Since λ⋆ = 2E{1−λ , the condition λ > λ⋆ is equivalent to k−m+1 (R)}

2Eλ{1 − λk−m+1 (R)} > λ1 (B) − λm (A). and therefore am (λ) = λm (A) + 2Eλ > λ1 (B) + 2Eλ λk−m+1 (R) ≥ bk−m+1 (λ). Hence at most k − m eigenvalues from the non-invariant coordinate block are at least am (λ), whereas all m eigenvalues from the invariant coordinate block are at least am (λ). Therefore every top-k ,λ ΠS⋆ = ΠS⋆ . eigenspace of Mλ contains S⋆ , that is, Π⌣ k Step 4: show that the leading m directions span S ⋆ for λ > λ⋆⋆ (part (ii)). Again by Weyl’s inequality [19, Theorem 4.3.7] used in decreasing-order notation, with indices (1, 1), b1 (λ) ≤ λ1 (B) + 2Eλ λ1 (R). 1 (B)−λm (A) Since λ⋆⋆ = λ2E{1−λ , the condition λ > λ⋆⋆ is equivalent to 1 (R)}

2Eλ{1 − λ1 (R)} > λ1 (B) − λm (A), and therefore am (λ) = λm (A) + 2Eλ > λ1 (B) + 2Eλ λ1 (R) ≥ b1 (λ). Thus every eigenvalue of the invariant block strictly exceeds every eigenvalue of the non-invariant block. Consequently, the leading m-dimensional eigenspace of Mλ is S⋆ .

Remark B.4 (Basis-level interpretation). Part (ii) has a basis-level interpretation: the first m directions of any ordered AnchorPCAλ representative basis (in the sense of § B.1) span S⋆ and diagonalize Σ within S⋆ . Under the weaker threshold in part (i), only the subspace-level statement S⋆ ⊆ Sk⌣,λ is guaranteed. If the eigenvalues of an invariant block and a non-invariant block of Mλ tie, the ordered eigenbasis is determined only up to a rotation inside the tied eigenspace. Hence, the k ordered directions need not contain m vectors spanning S⋆ . Specifically, one can always choose a basis whose first m columns span S⋆ , but this need not hold for every ordered basis. 19

B.2.3

Invariance of S⋆ to variance inflation

T (e) Proposition B.5. Assume λk (Σe ) > λk+1 (Σe ) for every e ∈ E, and let S⋆ = e∈E Im(Πk ) be the maximal invariant subspace of (Σe )e∈E . Then S⋆ is also the maximal invariant subspace of (e) (Σe + ∆e )e∈E for every choice of symmetric matrices ∆e with 0 ⪯ ∆e and Im(∆e ) ⊆ Im(Πk ), for all e ∈ E. Proof. We first show that each admissible perturbation leaves the rank-k principal subspaces unchanged. Intersecting these unchanged subspaces then shows that S⋆ is unchanged. Fix e ∈ E and (e) write Se := Im(Πk ). Since ∆e is symmetric with Im(∆e ) ⊆ Se , its kernel contains Se⊥ , so both Σe and ∆e are block-diagonal with respect to Rp = Se ⊕ Se⊥ . On Se⊥ , the perturbed covariance Σe + ∆e agrees with Σe and so has largest eigenvalue λk+1 (Σe ). On Se , since ∆e ⪰ 0, the smallest eigenvalue of Σe + ∆e is at least λk (Σe ). The gap λk (Σe ) > λk+1 (Σe ) then forces the rank-k principal subspace of Σe + ∆e to be exactly Se , so its rank-k principal projector, which we denote (e) (e) by Πk,∆ equals Πk . Intersecting over e gives that the maximal invariant subspace of the perturbed PE (e) covariances is S⋆ , and averaging gives E −1 e=1 Πk,∆ = Π.

B.3

Test statistic of FindS⋆

We now state the Wald-type test statistic used in the sequential procedure for estimating dim(S⋆ ). Since our construction is adapted from Schott’s partial common principal component subspace test [37], we first recall the standard version in our notation and then explain the modification for S⋆ . We present only the reduced representation, since this is the form used in computation. To simplify notation, throughout this appendix section we hide the dependence on k. For example, we (e) let Πe := Πk be the projector onto the top-k eigenspace in domain e. B.3.1

Standard construction.

Suppose that, for each e ∈ {1, . . . , E}, we observe an i.i.d. Gaussian sample with population covariance matrix Σe , and assume that the kth and (k +1)st eigenvalues of Σe are separated by a strict eigengap. Let Πe denote the projector onto the top-k eigenspace of Σe , and define r as the dimensionof the union of the top-k eigenspaces of (Σe )e∈E , that is, r := dim span Im(Π1 ) ∪ · · · ∪ Im(ΠE ) . Schott considers the hypotheses H0,s : r = s

HA,s : r > s.

versus

b e be the sample projector onto the top-k We now show how the test statistic is constructed. Let Π b e , and let Π b 0 be the projector onto the span of the eigenspace of the sample covariance matrix Σ b1 + · · · + Π b E . Let C be ∈ Rp×k be an orthonormal basis of Im(Π b e ), and let top s eigenvectors of Π p×(p−s) b b be an orthonormal basis of Im(Ip − Π0 ). The statistic depends on three objects: the C0 ∈ R c ⋆ , and the projector Fb onto the deviation vector v ⋆ , its estimated asymptotic covariance matrix W c ⋆ corresponding to its ns = (Ek − s)(p − s) largest eigenvalues. eigenspace of W Deviation vector v ⋆ .

Define

 b⊤ C b1 ) vec(C 0   .. v ⋆ :=  . . b⊤ C bE ) vec(C 

c⋆. Covariance matrix W

Here, we need several intermediate components. First, let s s X X + bj qbj qb⊤ , b−1 qbj qb⊤ , b ⋆ := b Π λ Π := λ j ⋆ j j j=1

j=1

b1 ≥ · · · ≥ λ bs and qb1 , . . . , qbs are the top-s eigenvalues and corresponding eigenvectors of where λ be,1 ≥ · · · ≥ λ be,p and qbe,1 , . . . , qbe,p denote the b b Π1 + · · · + ΠE . Further, for all e ∈ {1, . . . , E}, let λ 20

b e , in decreasing order. Then, for all e ∈ {1, . . . , E}, eigenvalues and corresponding eigenvectors of Σ let p k X X be,j λ be,ℓ nλ ⋆ ⊤ b b b0⊤ qbe,ℓ qbe,ℓ Ye := cj c⊤ ⊗C C0 , j b b 2 j=1 ℓ=k+1 ne (λe,j − λe,ℓ ) PE where cj is the jth standard basis vector in Rk , n := e=1 ne , and ⊗ denotes the Kronecker product. be,j − λ be,ℓ )−2 : if the gap is small, the estimated The eigengap enters here through the factors (λ eigenspace becomes unstable. Now, for h, i ∈ {1, . . . , E}, define the k(p − s) × k(p − s) matrix Vbhi :=

E X

bh⊤ Π b+ b b ⋆ b⊤ b + b b⊤ b + b b ⋆ b ⋆ b⊤ b + b (C ⋆ Cf ⊗Ip−s )Yf (Cf Π⋆ Ci ⊗Ip−s )−(Ch Π⋆ Ci ⊗Ip−s )Yi −Yh (Ch Π⋆ Ci ⊗Ip−s )

f =1

and let the Ek(p − s) × Ek(p − s) partitioned matrix whose (h, i)th block is Vbhi be denoted by Vb . Finally, let c ⋆ := diag(Yb1⋆ , . . . , YbE⋆ ) + Vb . W c ⋆ corresponding to its Projector onto the leading eigenspace. Let Fb be the eigenprojection of W ns = (Ek − s)(p − s) largest eigenvalues. Reduced Wald statistic.

Schott’s reduced Wald statistic is then c ⋆ Fb)+ v ⋆ , Ts := n (v ⋆ )⊤ (FbW

where (·)+ denotes the Moore–Penrose pseudoinverse. Under the Gaussian assumptions and the eigengap condition, Ts is asymptotically χ2ns under H0,s . B.3.2

Adaptation to S⋆

For FindS⋆ , we apply the same construction not to the top-k projectors Πe , but to their orthogonal complements Π⊥ e := Ip − Πe . A vector belongs to S⋆ if and only if it lies in the null space of every Π⊥ e . Hence S⋆ is the common ⊥ null space of Π⊥ , . . . , Π , and estimating dim(S ) is equivalent to estimating the dimension of the ⋆ 1 E span of these complementary projectors. In practice, this means that Schott’s statistic is applied exactly as above, but with the selected eigenspaces replaced by the (p − k)-dimensional complementary eigenspaces, that is, by Π⊥ e , and for t = p − k, . . . , p − mmin − 1, test H0,t : r = t

versus

HA,t : r > t,

r := p − dim(S⋆ ).

Let rb be the first t for which we do not reject. The estimated invariant dimension is then m b = p − rb. This is the same procedure as in the main text under the change of variables s = p − t. be is formed from the empirical eigenvectors qbe,k+1 , . . . , qbe,p . For the complementary-space test, C Accordingly, in the covariance formula the selected indices range over j = k + 1, . . . , p and the non-selected indices over ℓ = 1, . . . , k; equivalently, one may apply Schott’s formula after reordering the spectrum so that the complementary eigenspace is treated as the selected subspace. Under the complement parametrization, the degrees of freedom become nt = (E(p − k) − t)(p − t), so Tt is asymptotically χ2nt under H0,t . The eigengap assumption λk (Σe ) > λk+1 (Σe ) separates the top-k and bottom-(p − k) eigenspaces symmetrically, so Schott’s regularity conditions apply verbatim to the complementary eigenspaces.

Proofs

C.1

Additional intermediate results

Here we collect additional results used in the proofs of the main text results. We first introduce principal angles and the corresponding projector-distance identity, then restate the Davis–Kahan theorem, and finally record the three lemmas used in the consistency proof § C.4. 21

Definition C.1 (Principal angles). Let U, V ∈ Op×k and let σ1 ≥ · · · ≥ σd be the singular values of U ⊤ V . The principal angles between U and V are cos−1 (σ1 ), . . . , cos−1 (σd ). We write Θ(U, V ) for the diagonal matrix of principal angles and sin Θ(U, V ) for its entrywise sine. Proposition C.2 (Projector distance and principal angles). Let U, V ⊆ Rp be subspaces of equal dimension d, and let U, V ∈ Rp×d have orthonormal columns spanning U and V, respectively. Then ∥ΠU − ΠV ∥op = ∥sin Θ(U, V )∥op .

(7)

Proof. This is the standard projector-distance identity for equal-dimensional subspaces; see, for example, [9, Thm. 2.1]. Since U and V have orthonormal columns, their orthogonal projectors are ΠU = U U ⊤ and ΠV = V V ⊤ . b ∈ Rp×p be symmetric, with eigenvalues Theorem C.3 (Davis–Kahan sin Θ theorem [47]). Let Σ, Σ b b λ1 ≥ · · · ≥ λp and λ1 ≥ · · · ≥ λp and let r, s ∈ {1, . . . , p} such that r ≤ s. Let V and Vb contain br , . . . , λ bs , respectively. orthonormal eigenvectors associated with the eigenvalues λr , . . . , λs and λ Let b − λ| : λ ∈ [λs , λr ], λ b ∈ (−∞, λ bs+1 ] ∪ [λ br−1 , ∞)}, δ := inf{|λ b0 := ∞ and λ bp+1 := −∞. If δ > 0, then with λ b −Σ Σ sin Θ(Vb , V )

≤ op

(8)

b be symmetric p × p matrices, and let Lemma C.4 (Top-d eigenspace perturbation bound). Let A, A b denote the orthogonal projector onto the spans of the top d d ∈ {1, . . . , p−1}. Let Πd (A) and Πd (A) b b−A eigenvectors of A and A, respectively. If γd (A) := λd (A) − λd+1 (A) > 0 and A < γd (A) , op

then b−A 2 A b − Πd (A) Πd (A)

b−A Proof. Let ε := A

≤ op

γd (A)

. By Weyl’s inequality [19, Theorem 4.3.7] in decreasing-order notation

b − A), λd+1 (A) b ≤ λd+1 (A) + ε, so with indices (d + 1, 1) applied to (A, A b ≥ λd (A) − λd+1 (A) − ε = γd (A) − ε > 0, λd (A) − λd+1 (A) where the last inequality uses ε < γd (A)/2 which holds by assumption. Let V and Vb have b respectively. Applying Thm. C.3 orthonormal columns spanning the top-d eigenspaces of A and A, b − λ| : λ ∈ b = A, b r = 1, and s = d. The separation quantity is δ = inf{|λ with Σ = A, Σ b b b [λd (A), λ1 (A)], λ ∈ (−∞, λd+1 (A)]} = λd (A) − λd+1 (A) ≥ γd (A) − ε > γd (A)/2. Therefore, b − Πd (A) Πd (A)

= sin Θ(Vb , V )

≤ op

ε 2ε ≤ , δ γd (A)

where the first equality follows from (7). Recall the following asymptotic notation. Let (bN )N ≥1 be a deterministic positive sequence and let (RN )N ≥1 be a sequence of random scalars or matrices. Fix a norm ∥ · ∥ on the relevant space; for scalars this is the absolute value. We write RN = O(bN ) almost surely if, there exists an event of probability one on which, for each realization, there exist a finite constant C < ∞ and index N0 ∈ N+ (both allowed to depend on the realization) such that ∥RN ∥ ≤ CbN for all N ≥ N0 . We write RN = o(bN ) almost surely if ∥RN ∥/bN → 0 almost surely. When the notation is applied entrywise, it is understood to hold for each fixed entry separately. 22

Lemma C.5 (Polynomial rate for the averaged projectors). Assume that, for all e ∈ E, E[xe ] = 0, E∥xe ∥42 < ∞, and λk (Σe ) > λk+1 (Σe ). Then, for all α ∈ (0, 1/2) and all e ∈ E, b e − Σe Σ

= o(n−α e )

and

b (e) − Π(e) Π k k

= o(n−α e )

almost surely. Consequently, with n := mine∈E ne , b −Σ Σ

= o(n−α )

and

b −Π Π

= o(n−α ) op

almost surely. Proof. Step 1: convergence of coordinates in each domain. Let e ∈ E and α ∈ (0, 1/2). In this proof only, let N := ne and Xi := xe,i . Choose s ∈ (1, 2) such that α < 1 − 1/s. For all a, b ∈ {1, . . . , p} and all i ∈ {1, . . . , N }, let Yiab := (Xi )a (Xi )b (where (Xi )a denotes the a-th coordinate of Xi ). Since s < 2 and E∥xe ∥42 < ∞, we have E|Yiab |s = E[|(Xi )a |s |(Xi )b |s ] ≤ E∥Xi ∥2s 2 <∞ where the final bound uses 2s < 4. By the Marcinkiewicz–Zygmund strong law of large numbers [4, Theorem 8.4.4, statement (ii)], N 1 X ab Y − EYiab N 1/s i=1 i converges almost surely, and is therefore bounded almost surely. Dividing by N 1−1/s gives N 1 X ab Yi − EYiab = O N −(1−1/s) = o(N −α ) almost surely, N i=1 where the last equality uses α < 1 − 1/s. Applying the same strong law to the coordinate variables PN (Xi )a and using E(Xi )a = 0, gives X a := N1 i=1 (Xi )a = O(N −(1−1/s) ) almost surely, and hence X a X b = O(N −2(1−1/s) ) = o(N −α ) almost surely, since 2(1 − 1/s) > α. PN Step 2: convergence of covariances in each domain. Let AN := N −1 i=1 Xi Xi⊤ . Since b e = N AN − X X ⊤ , we have Σ N −1 ⊤ b e − Σe = (AN − Σe ) − X X ⊤ + 1 Σ AN − X X . N −1 The first two terms are entrywise o(N −α ) almost surely by the bounds of Step 1: for AN − Σe this uses E[Xi Xi⊤ ] = Σe , which holds since E[Xi ] = 0, so that the entrywise bound on AN from Step 1 is a bound on AN − Σe . The last term is entrywise O(N −1 ) = o(N −α ) almost surely, because ⊤ AN − X X = O(1) entrywise almost surely and α < 1. Since p is fixed, entrywise convergence at −α rate o(N ) implies operator-norm convergence at the same rate. Therefore b e − Σe Σ

= o(n−α e )

almost surely.

(9)

Step 3: convergence of projection matrices in each domain. Let γk (Σe ) := λk (Σe )−λk+1 (Σe ) > b e − Σe 0, by assumption. By (9), for all sufficiently large ne , Σ < γk (Σe )/2 almost surely. op

b=Σ b e , and d = k gives Applying Lemma C.4 with A = Σe , A 2 b (e) − Π(e) b e − Σe Π ≤ Σ = o(n−α e ) k k γk (Σe ) op op where the last equality follows from (9).

almost surely,

Step 4: convergence of the average covariance and average projection matrix. Since E is fixed and ne ≥ n for all e ∈ E, averaging over e ∈ E yields E

b −Σ Σ

1 X b ≤ Σe − Σ e = o(n−α ) E e=1 op op 23

and

b −Π Π

≤ op

1 X b (e) (e) Π − Πk = o(n−α ) E e=1 k op

almost surely. Lemma C.6 (Recovery of the eigenvalue blocks of Π). Let 1 ≥ ρ1 > · · · > ρL ≥ 0 be the distinct eigenvalues of Π, and let Πℓ denote the orthogonal projector onto the eigenspace of ρℓ , with multiplicity dℓ := rank(Πℓ ). Listing the eigenvalues of Π in nonincreasing order with multiplicity, ρℓ occupies a block of dℓ consecutive indices which we denote Iℓ , so that I1 , . . . , IL partition {1, . . . , p}. Let γΠ be the smallest gap between the distinct population eigenvalues, ( min{ρ − ρ′ : ρ, ρ′ ∈ {ρ1 , . . . , ρL }, ρ > ρ′ }, L ≥ 2, γΠ := +∞, L = 1. b (including duplicates). Suppose that ∥Π b − Π∥op < toln /4 Let ρb1 ≥ · · · ≥ ρbp be the eigenvalues of Π and toln < γΠ /2. Then the block-stabilized construction of AnchorPCA∞ in § 5.1 partitions the indices 1, . . . , p into empirical blocks that coincide with I1 , . . . , IL . Proof. Let cℓ :=

j<ℓ dj , so that Iℓ = {cℓ + 1, . . . , cℓ + dℓ }, and let εn := ∥Π − Π∥op . Recall that

b from largest to smallest to group them the construction in § 5.1 scans the ordered eigenvalues of Π into empirical blocks: each block starts at an eigenvalue and includes subsequent eigenvalues as long b − Π∥op and let as they are at most toln smaller than the first eigenvalue of the block Let εn := ∥Π ℓ ∈ {1, . . . , L − 1}. Step 1: no block Iℓ is split. By Weyl’s inequality [19, Theorem 4.3.1] in decreasing-order notation b − Π, we have that the eigenvalues of Π b and Π, applied to the two symmetric matrices Π and Π b b matched by rank order and counted with multiplicity, satisfy |λj (Π) − λj (Π)| ≤ ∥Π − Π∥op = εn for all j ∈ {1, . . . , p}. Since λj (Π) = ρℓ for all j ∈ Iℓ , it follows that |b ρj − ρℓ | ≤ εn for all j ∈ Iℓ . Hence, for all j, j ′ ∈ Iℓ , |b ρj − ρbj ′ | ≤ |b ρj − ρℓ | + |b ρj ′ − ρℓ | ≤ 2εn < toln /2 < toln , b − Π∥op < toln /4. In particular, where the third inequality follows by the assumption εn = ∥Π every eigenvalue in Iℓ lies within toln of the first eigenvalue of the block, ρbcℓ +1 , so the rule keeps all indices in Iℓ in the same empirical block. Step 2: adjacent blocks Iℓ and Iℓ+1 are not merged. Applying Weyl’s inequality [19, Theob b rem 4.3.1] in decreasing-order notation to Π and Π−Π, for any index j we have λj (Π)+λp (Π−Π) ≤ b b b λj (Π) ≤ λj (Π) + λ1 (Π − Π), and hence λj (Π) − εn ≤ λj (Π) ≤ λj (Π) + εn . Taking j = cℓ + 1 b ≥ ρℓ − εn . Similarly, taking j = cℓ + dℓ + 1 and and using λcℓ +1 (Π) = ρℓ gives ρbcℓ +1 = λcℓ +1 (Π) b ≤ ρℓ+1 + εn . Therefore, using λcℓ +dℓ +1 (Π) = ρℓ+1 gives ρbcℓ +dℓ +1 = λcℓ +dℓ +1 (Π) toln ρbcℓ +1 − ρbcℓ +dℓ +1 ≥ ρℓ − ρℓ+1 − 2εn ≥ γΠ − 2εn > γΠ − > toln , 2 where the second last inequality uses that εn < toln /4, and the final inequality uses toln < γΠ /2 (both hold by assumption). Thus the first eigenvalue past Iℓ lies more than toln below ρbcℓ +1 , the first eigenvalue of the current block, so the rule closes that block at index cℓ + dℓ and starts a new one. Going through the eigenvalues from largest to smallest, Steps 1 and 2 show that each block starts at index cℓ + 1 and ends at index cℓ + dℓ : Step 1 keeps a block from ending too early, and Step 2 forces it to end at the right place. The blocks the construction produces are therefore exactly I1 , . . . , IL . C.2

Proof of Thm. 3.1

Thm. 3.1 is a projector-level statement. We prove the following stronger version, which uses the ordered convention of § B.1 and additionally records what happens to the ordered directions. Part (i) concerns AnchorPCA∞ , and Part (ii) concerns AnchorPCAλ as λ → ∞; both are stated first at the projector level and then at the direction level under a simple-spectrum Assm. C.7. Together, the two parts imply Thm. 3.1. If m = 0, the statements below are vacuous, so we assume m ≥ 1. 24

Assumption C.7 (Simple spectrum on the invariant block). The restriction of ΣS⋆ to S⋆ (i.e., ΣS⋆ := ΠS⋆ ΣΠS⋆ ) has simple spectrum within S⋆ : there exist µ1 > · · · > µm and an orthonormal basis s1 , . . . , sm ∈ S⋆ such that ΣS⋆ si = µi si for all i ∈ {1, . . . , m}. Theorem C.8 ( Thm. 3.1, extended version). Assume 1 ≤ m := dim(S⋆ ) ≤ k < p, and let ,∞ ,λ ΣS⋆ = ΠS⋆ ΣΠS⋆ . Let Π⌣ and Π⌣ be any rank-k solution projectors of AnchorPCA∞ and k k AnchorPCAλ , respectively, and let Wk⌣,∞ = (u1 , . . . , uk ),

Wk⌣,λ = (w1λ , . . . , wkλ )

be corresponding ordered bases (§ B.1). (i) (AnchorPCA∞ .) The first m vectors u1 , . . . , um form an ordered orthonormal basis of S⋆ ,∞ that diagonalizes ΣS⋆ within S⋆ ; in particular, span(u1 , . . . , um ) = S⋆ and Π⌣ ΠS⋆ = k ΠS⋆ . If additionally Assm. C.7 holds, then ui ∈ {+si , −si } for all i ∈ {1, . . . , m}. λ ,λ ), that is, the denote the projector onto span(w1λ , . . . , wm (ii) (AnchorPCA, λ → ∞.) Let Π⌣ m leading rank-m spectral projector of Mλ . Then, ,λ Π⌣ − ΠS⋆ op → 0 m

as λ → ∞.

If additionally Assm. C.7 holds, then, with a suitable choice of signs, as λ → ∞, wiλ → si

λi (Mλ ) − 2Eλ → µi ,

for all i ∈ {1, . . . , m}.

Proof. Part (i) AnchorPCA∞ . Since S⋆ is the eigenspace of Π corresponding to its largest eigenvalue 1, and dim(S⋆ ) = m ≤ k, the ordering convention of § B.1 draws the first m directions u1 , . . . , um entirely from S⋆ . Hence span(u1 , . . . , um ) = S⋆ ,

,∞ Π⌣ ΠS⋆ = ΠS⋆ . k

Within this block, the convention orders directions by diagonalizing Σ, equivalently ΣS⋆ , within S⋆ . Thus u1 , . . . , um is an orthonormal eigenbasis of ΣS⋆ on S⋆ , ordered by decreasing eigenvalue. Under Assm. C.7, this eigenbasis is unique up to sign, so ui ∈ {+si , −si } for each i ∈ {1, . . . , m}. Part (ii) AnchorPCAλ . We first show the claimed subspace convergence in steps 1–2 invoking the Davis–Kahan theorem (Thm. C.3). The direction-level convergence under Assm. C.7 is proved in steps 3–6. Step 1: Uniqueness of the top-m eigenspace of (2Eλ)−1 Mλ . Let S ∈ Rp×m have orthonormal columns spanning S⋆ , and complete it to an orthogonal matrix (S, Z) with Z ∈ Rp×(p−m) spanning S⋆⊥ . Since S⋆ is the eigenvalue-1 eigenspace of Π, we have I 0 (S, Z)⊤ Π(S, Z) = m , R := Z ⊤ ΠZ. 0 R fλ := 1 Mλ = Π + 1 Σ. By Weyl’s inequality [19, Theorem 4.3.1] in decreasingLet M 2Eλ 2Eλ order notation applied to (Π, (2Eλ)−1 Σ), we have, for each index j, λj (Π) − ∥(2Eλ)−1 Σ∥op ≤ fλ ) ≤ λj (Π) + ∥(2Eλ)−1 Σ∥op . Since ∥(2Eλ)−1 Σ∥op = (2Eλ)−1 ∥Σ∥op , applying this with λ j (M fλ ) ≥ 1 − (2Eλ)−1 ∥Σ∥op . Applying the same bound with j = m and using λm (Π) = 1 gives λm (M fλ ) ≤ λ1 (R) + (2Eλ)−1 ∥Σ∥op . Hence, for j = m + 1 and using λm+1 (Π) = λ1 (R) gives λm+1 (M Σ ∥ ∥ fλ ) > λm+1 (M fλ ). Thus, the leading rank-m spectral all λ >:= E(1−λ1op(R)) =: λ0 ,, we have λm (M projector of Mλ is unique for all λ > λ0 . Step 2: Apply Davis–Kahan to obtain subspace convergence. We have that 1 Σ op > 0, (10) 2Eλ where the first inequality follows from step 1 and the second holds for all λ > λ0 as λ1 (R) < 1 (shown b := M fλ with r = 1 and in step 2 in the proof of Prop. B.3). We apply Thm. C.3 to Σ := Π and Σ fλ ) ≥ 1 − λ1 (R) − 1 − λm+1 (M

λ λ s = m. Hence, V = S, Vb = (w1λ , . . . , wm ) =: Wm . With these choices, the separation parameter b − λ| : λ ∈ [λm (Π) = 1, λs (Π) = 1], λ b ∈ (−∞, λm+1 (M fλ )]} = |1 − λm+1 (M fλ )|, so δ = inf{|λ the theorem states that

λ sin Θ(Wm , S) op ≤

fλ Π−M

fλ )| |1 − λm+1 (M

We apply (10) for all λ > λ0 and then consider λ → ∞, which gives λ sin Θ(Wm , S) op ≤

(2Eλ)−1 Σ op 1 − λ1 (R) − (2Eλ)−1 Σ op

→ 0.

,λ λ λ ⊤ Since Π⌣ = Wm (Wm ) and ΠS⋆ = SS ⊤ , by (7), as λ → ∞, m ,λ Π⌣ − ΠS⋆ op → 0. m

Step 3: Write the eigenvalue equation in S⋆ ⊕ S⋆⊥ coordinates. Assume now, in addition, that Assm. C.7 holds. For the remainder of the proof, set S := (s1 , . . . , sm ) and complete it to an orthogonal matrix (S, Z). In the basis (S, Z), 2EλIm + A C ⊤ (S, Z) Mλ (S, Z) = , C⊤ 2EλR + B where A := S ⊤ ΣS⋆ S = S ⊤ ΣS = diag(µ1 , . . . , µm ),

B := Z ⊤ ΣZ,

C := S ⊤ ΣZ.

The equality S ⊤ ΣS⋆ S = S ⊤ ΣS follows from ΠS⋆ S = S. Let i ∈ {1, . . . , m}. Recall that Mλ wiλ = λi (Mλ )wiλ , i.e. wiλ is the ith leading eigenvector of Mλ . Since the columns of (S, Z) form an orthonormal basis of Rp , the vector wiλ has unique coordinates uλi ∈ Rm and viλ ∈ Rp−m in this basis, so that wiλ = Suλi + Zviλ . We also let ηi (λ) := λi (Mλ ) − 2Eλ. The eigenvalue equation Mλ wiλ = λi (Mλ )wiλ then is equivalent to (A − ηi (λ)Im )uλi + Cviλ = 0, C ⊤ uλi − 2Eλ(Ip−m − R) + ηi (λ)Ip−m − B viλ = 0. | {z }

(11) (12)

=:Li (λ)

Step 4: Reduce to an eigenvalue equation in S⋆ . Rearranging (12) gives Li (λ)viλ = C ⊤ uλi . Because R ⪯ λ1 (R)Ip−m with λ1 (R) < 1, we have Ip−m − R ⪰ (1 − λ1 (R))Ip−m . By Weyl’s inequality [19, Theorem 4.3.1], used in decreasing-order notation and applied to (2EλΠ, Σ), for each i ≤ m we have λi (2EλΠ) + λp (Σ) ≤ λi (Mλ ) ≤ λi (2EλΠ) + λ1 (Σ). Since S⋆ is the eigenvalue-1 eigenspace of Π, λi (2EλΠ) = 2Eλ for i ≤ m. Hence λp (Σ) ≤ λi (Mλ ) − 2Eλ ≤ λ1 (Σ), and therefore |ηi (λ)| ≤ ∥Σ∥op . (13) Combining this with ∥B∥op ≤ Σ op , we have Li (λ) = 2Eλ(Ip−m − R) + ηi (λ)Ip−m − B ⪰ 2Eλ(1 − λ1 (R)) − 2 Σ op Ip−m . Thus, for all λ > λ0 , the matrix Li (λ) is positive definite and therefore invertible. In addition, for all λ > λ0 , 1 . (14) Li (λ)−1 op ≤ 2Eλ(1 − λ1 (R)) − 2 Σ op 26

Since Li (λ)viλ = C ⊤ uλi and ∥uλi ∥2 ≤ 1, it follows that ∥viλ ∥2 ≤ Li (λ)−1 op ∥C∥op ∥uλi ∥2 ≤ Li (λ)−1 op ∥C∥op . Hence, ∥viλ ∥2 → 0 as λ → ∞ by (14). Because wiλ = Suλi + Zviλ has unit norm and (S, Z) is orthogonal, we have ∥uλi ∥22 + ∥viλ ∥22 = 1. Therefore ∥uλi ∥2 → 1

(15)

as λ → ∞. Substituting viλ = Li (λ)−1 C ⊤ uλi into (11) gives (A − ηi (λ)Im )uλi + CLi (λ)−1 C ⊤ uλi = 0, or equivalently (A + CLi (λ)−1 C ⊤ )uλi = ηi (λ)uλi . (16) Step 5: Identify the limiting eigenvalues and directions along subsequences. Fix any sequence λn → ∞. For all i ∈ {1, . . . , m} , the sequence ηi (λn ) is bounded (see (13)). By (15), after discarding finitely many terms we have ∥uλi n ∥2 > 0 for all i ∈ {1, . . . , m}. For such n, define uλn

u bλi n := ∥uλin ∥ , which belongs to the compact unit sphere of Rm . By the Bolzano–Weierstrass i

theorem, each of these sequences has a convergent subsequence. Passing to a common subsequence (which we continue to index by n) along which all 2m sequences converge, it holds that for all i ∈ {1, . . . , m}, there exist ηi⋆ ∈ R and a unit vector u ei ∈ Rm such that u bλi n → u ei .

ηi (λn ) → ηi⋆ ,

Dividing the reduced eigenvalue equation (16) by ∥uλi n ∥2 gives λn A + CLi (λn )−1 C ⊤ u bi = ηi (λn )b uλi n . Passing to the limit λn → ∞ along the chosen subsequence, the term CLi (λn )−1 C ⊤ vanishes since 2 CLi (λ)−1 C ⊤ op ≤ ∥C∥op Li (λ)−1 op → 0 by (14), and we obtain Ae ui = ηi⋆ u ei . This implies that u ei is a unit eigenvector of A with eigenvalue ηi⋆ . Since A = diag(µ1 , . . . , µm ) with µ1 > · · · > µm , there exists a unique index ji ∈ {1, . . . , m} such that ηi⋆ = µji ,

u ei = ±eji .

It remains to show that ji = i, which we do using orthogonality and eigenvalue ordering. Let i, ℓ ∈ {1, . . . , m} with i ̸= ℓ. Orthogonality of wiλn and wℓλn implies 0 = ⟨wiλn , wℓλn ⟩ = ⟨uλi n , uλℓ n ⟩ + ⟨viλn , vℓλn ⟩. Since viλn , vℓλn → 0 and ∥uλi n ∥2 , ∥uλℓ n ∥2 → 1, passing to the limit yields ⟨e ui , u eℓ ⟩ = 0. Hence j1 , . . . , jm is a permutation of {1, . . . , m}. On the other hand, by the ordering convention, λ1 (Mλn ) ≥ · · · ≥ λm (Mλn )

=⇒

η1 (λn ) ≥ · · · ≥ ηm (λn ).

Passing to the limit gives µj1 ≥ · · · ≥ µjm . Because µ1 > · · · > µm are strictly decreasing, the only possible permutation (j1 , . . . , jm ) of {1, . . . , m} for which µj1 ≥ · · · ≥ µjm holds is the identity. Therefore, for all i ∈ {1, . . . , m} ηi⋆ = µi ,

u ei = ±ei .

Step 6: Full convergence of eigenvalues and eigenvectors. Since the original sequence λn → ∞ was arbitrary, every convergent subsequence of ηi (λ) has limit µi . A bounded sequence whose convergent subsequences all share the same limit converges to that limit, so for all i ∈ {1, . . . , m} λi (Mλ ) − 2Eλ = ηi (λ) → µi

λ → ∞.

Similarly, every subsequential limit of u bλi equals ±ei . We can choose the sign of wiλ so that λ ⟨wi , si ⟩ ≥ 0 rules out the −ei limit and u bλi → ei . As ∥uλi ∥2 → 1 by (15), this gives uλi → ei . Since λ vi → 0, wiλ = Suλi + Zviλ → Sei = si as λ → ∞, which is the asserted direction-level convergence. 27

C.3

Proof of Thm. 4.1

We first show (4) and then use it to prove statements (i) and (ii). Let W ∈ Op×k and ΠW := W W ⊤ . By definition of Cλ , for every perturbed covariance Σ′e , there exists ∆e ∈ Sp+ with 0 ⪯ ∆e ⪯ (e) 2Eλ Πk such that Σ′e = Σe + ∆e . Hence,  max ′ ′

(Σ1 ,...,ΣE )∈Cλ

E 1 X Tr Σe (Ip − ΠW ) + E e=1 

R(W ; Σ′1 , . . . , Σ′E ) =

 Tr ∆e (Ip − ΠW )  .

max

∆e ∈Sp +, (e)

0⪯∆e ⪯2Eλ Πk

Applying the trace-monotonicity inequality Tr(AC) ≤ Tr(BC) (valid for A ⪯ B and C ⪰ 0; see, (e) e.g., Boyd and Vandenberghe [7, Example 2.24]) with A = ∆e , B = 2Eλ Πk , and C = Ip − ΠW , we obtain (e) Tr ∆e (Ip − ΠW ) ≤ Tr 2Eλ Πk (Ip − ΠW ) , (e)

with equality for ∆⋆e = 2Eλ Πk . Therefore, max ′ ′

(Σ1 ,...,ΣE )∈Cλ

R(W ; Σ′1 , . . . , Σ′E ) = =

E E X 1 X (e) Tr Σe (Ip − ΠW ) + 2λ Tr Πk (Ip − ΠW ) E e=1 e=1 E E X 1 X (e) Tr Σe (Ip − ΠW ) + λ ∥ΠW − Πk ∥2F , E e=1 e=1 (e)

where the second line follows since ΠW and Πk are rank-k orthogonal projectors, so (e) (e) (e) (e) 2 Tr Πk (Ip − ΠW ) = 2k − 2 Tr(ΠW Πk ) = Tr (ΠW − Πk )2 = ∥ΠW − Πk ∥2F . This proves (4). PE Part (i). A direct computation gives E1 e=1 Tr Σe (Ip − ΠW ) = Tr(Σ) − Tr(W ⊤ ΣW ) so minimizing the worst-case reconstruction error is equivalent to maximizing ⊤

Tr(W ΣW ) − λ

E X

(e)

∥ΠW − Πk ∥2F ,

e=1

which is Definition 2.3. (e)

(e)

Part (ii). Using the fact that ∥ΠW − Πk ∥2F = 2k − 2 Tr(W ⊤ Πk W ) for all e ∈ E and summing over e, E X (e) ∥ΠW − Πk ∥2F = 2Ek − 2E Tr(W ⊤ ΠW ). e=1

Hence B = arg min

E X

f ∈Op×k e=1 W

(e)

2 f⊤ f ∥ΠW f − Πk ∥F = arg max Tr(W ΠW ). f ∈Op×k W

Combining this with R(W ; Σ1 , . . . , ΣE ) = Tr(Σ) − Tr(W ⊤ ΣW ), minimizing R over B amounts to maximizing Tr(W ⊤ ΣW ) over the maximizers of Tr(W ⊤ ΠW ), which is Definition 2.2. C.4

Consistency of Anchor PCA

We now state the consistency results of the empirical AnchorPCAλ estimator and of the blockstabilized empirical AnchorPCA∞ construction (§ 5.1), already stated in § 5.3 informally. Corollary C.11 then shows that, when S⋆ ̸= {0}, the first block returned by AnchorPCA∞ consistently estimates S⋆ . Throughout, for all e ∈ E, let xe,1 , xe,2 , . . . , xe,ne be i.i.d. copies of xe , and assume n := mine∈E ne → ∞. Theorem C.9 (Consistency of the empirical projector-agreement estimators). Assume that, for all e ∈ E, E[xe ] = 0, E∥xe ∥42 < ∞, and λk (Σe ) > λk+1 (Σe ). Then the following statements hold. 28



b ⌣,λ be an empirical rank-k AnchorPCAλ solution projector. If λk (Mλ ) > (i) Let λ > 0, and let Π k λk+1 (Mλ ) then b ⌣,λ − Π⌣,λ Π k k

→0

almost surely as n → ∞.

b ⌣,∞ be the block-stabilized AnchorPCA∞ estimator with tolerance toln > 0. Assume that (ii) Let Π k there exists α ∈ (0, 1/2) and c > 0 such that for all n, toln = cn−α and the population rank-k AnchorPCA∞ solution projector is unique, then b ⌣,∞ − Π⌣,∞ Π k k

→0

almost surely as n → ∞.

Remark C.10. The uniqueness assumption in part (ii) concerns the population projector, i.e. the optimal subspace, not any particular basis. For example, it holds automatically if there is a Πeigenvalue block boundary at rank-k. Otherwise, it requires an eigengap condition (21). To prove Thm. C.9, we use three auxiliary lemmas collected in § C.1: Lemma C.4 transfers operatornorm perturbation bounds to leading eigenspace projectors, Lemma C.5 gives polynomial almost-sure b and Π, b and Lemma C.6 guarantees recovery of the population spectral blocks of Π by the rates for Σ empirical block clustering. Proof of Thm. C.9. Choose any α0 ∈ (0, 1/2). By Lemma C.5, b −Σ Σ

→0 op

b −Π Π

and

→0 op

almost surely.

(17)

Statement (i): AnchorPCAλ . By (17), there is an event Ω0 with P(Ω0 ) = 1 on which b −Π cλ = Σ b + cλ − Mλ b −Σ M ≤ Σ + 2Eλ Π → 0 almost surely, where M op

b Mλ = Σ + 2Eλ Π. Let γk (Mλ ) := λk (Mλ ) − λk+1 (Mλ ); by assumption, γk (Mλ ) > 0. 2Eλ Π, cλ − Mλ Fix a realization in Ω0 . Then M < γk (Mλ )/2 for all n ≥ N (with N depending on the op

b=M cλ , d = k) gives realization), so Lemma C.4 (with A = Mλ , A 2 cλ − Mλ b ⌣,λ − Π⌣,λ M → 0. Π ≤ k k γk (Mλ ) op op As this holds on Ω0 , the convergence is almost sure. Statement (ii): AnchorPCA∞ . Recall that the block-stabilized construction of the solution to b from largest to smallest to group them into AnchorPCA∞ (§ 5.1 ) scans the ordered eigenvalues of Π empirical blocks: each empirical block starts at an eigenvalue and includes subsequent eigenvalues as long as they are at most toln smaller than the first eigenvalue of the block. Let 1 ≥ ρ1 > · · · > ρL ≥ 0 be the distinct eigenvalues of Π, and let Πℓ denote Pthe orthogonal projector onto the eigenspace corresponding to ρℓ . Let dℓ := rank(Πℓ ) and cℓ := j<ℓ dj . Listing the eigenvalues of Π in nonincreasing order with multiplicity, ρℓ occupies a block of dℓ consecutive b − Π∥op indices Iℓ := {cℓ + 1, . . . , cℓ + dℓ }, so that I1 , . . . , IL partition {1, . . . , p}. Let εn := ∥Π b and ξn := ∥Σ − Σ∥op . Step 1: Match the indices of the empirical and population blocks. By Lemma C.5, with the same α as in toln = cn−α , we have εn = o(n−α ) and ξn = o(n−α ) almost surely. Since toln = cn−α with c > 0, this implies εn = o(toln ) almost surely. We now work on an event Ω1 with P(Ω1 ) = 1 on which these convergences hold. Fix a realization in Ω1 ; every convergence in the remainder of the proof is along this realization, and each “for all sufficiently large n” refers to a realization-dependent threshold. Then, for all sufficiently large n, εn < toln /4.

(18)

Moreover, since toln decreases monotonically to 0 deterministically, for all sufficiently large n we also have toln < γΠ /2, (19) 29

where γΠ is the smallest gap between distinct eigenvalues of Π as defined in Lemma C.6. Hence Lemma C.6 applies, and the block-stabilized construction in § 5.1 partitions the indices 1, . . . , p into b ℓ denote the orthogonal empirical blocks that coincide with I1 , . . . , IL . For each ℓ ∈ {1, . . . , L}, let Π b with indices in Iℓ , i.e. onto Eig dℓ in the notation projector onto the span of the eigenvectors of Π d of § 5.1; Eigℓ estimates Im(Πℓ ), and it has dimension dℓ . b ℓ − Πℓ Step 2: Show consistency of each block projector. We show Π

→ 0 for all ℓ ∈

{1, . . . , L}. If L = 1, then Π has a single eigenvalue, so Π1 = Ip . For all sufficiently large n, the b 1 = Ip = Π1 , which proves the claim. empirical block construction also returns a single block, so Π b = Π, b r = cℓ + 1, Assume now that L ≥ 2, and let ℓ ∈ {1, . . . , L}. Apply Thm. C.3 with Σ = Π, Σ b b and s = cℓ + dℓ : then V spans Im(Πℓ ), V spans Im(Πℓ ), and [λs (Π), λr (Π)] = {ρℓ }. It remains to lower bound the separation δ between ρℓ and the empirical eigenvalues outside Iℓ (defined in Thm. C.3). By Weyl’s inequality [19, Theorem 4.3.1] in decreasing-order notation applied to Π and b − Π at indices cℓ and cℓ + dℓ + 1, Π ρbcℓ ≥ ρℓ−1 − εn

ρbcℓ +dℓ +1 ≤ ρℓ+1 + εn ,

and

where for ℓ = 1 there is no eigenvalue above the block and for ℓ = L none below, so only the relevant one-sided bound is used. In either case, the empirical eigenvalues outside Iℓ lie at distance at least γΠ − εn from ρℓ , so δ ≥ γΠ − εn > 0, where the last inequality follows from εn < toln /2 < γΠ /4 using (18) and (19). Thus, we get εn b ℓ − Πℓ Π = sin Θ(Vb , V ) → 0, ≤ γΠ − ε n op op b ℓ and Πℓ have equal rank dℓ (Step 1). where we additionally used Prop. C.2 for the first equality as Π Step 3: Assemble the rank-k projector. Let r ∈ {1, . . . , L} be the unique index such that cr < k ≤ cr + dr , and let q := k − cr , so, in the population case, blocks 1, . . . , r − 1 are taken in full and q directions are taken from block r. Case 1: q = dr . The first r blocks are fully taken for the population solution. By Steps 1a and 1b Pr ,∞ in § B.1.2, the population AnchorPCA∞ projector is Π⌣ = ℓ=1 Πℓ . For all sufficiently large n, k b ⌣,∞ = Pr Π b the empirical construction selects the first r empirical blocks, so Π ℓ=1 ℓ . By Step 2, k b ⌣,∞ − Π⌣,∞ Π k k

≤ op

r X ℓ=1

b ℓ − Πℓ Π

→ 0. op

Case 2: q < dr . Blocks 1, . . . , r − 1 are taken in full, and the secondary variance criterion selects br := Π b rΣ bΠ b r , and let Qr and Q b r denote the q directions from block r. Let Br := Πr ΣΠr and B b orthogonal projectors onto the spans of the top q eigenvectors of Br and Br , respectively. We now br → Br in operator norm. Note that show that B br − Br = Π b r (Σ b − Σ)Π b r + (Π b r − Πr )ΣΠ b r + Πr Σ(Π b r − Πr ), B br together with Π

= ∥Πr ∥op = 1 and (17) gives b r − Πr ∥op → 0. br − Br ∥op ≤ ξn + 2∥Σ∥op ∥Π ∥B

,∞ By Steps 1a and 1b in § B.1.2, the population AnchorPCA∞ projector is Π⌣ = k

(20) Pr−1

ℓ=1 Πℓ + Qr , where the first sum is a fixed projector independent of the selection within block r. Hence Qr = Pr−1 ,∞ Πk⌣,∞ − ℓ=1 Πℓ is determined by Π⌣ , so the uniqueness assumption forces the top-q eigenspace k of Br to be unique; equivalently, Br has a positive eigengap at rank q:

ηr := λq (Br ) − λq+1 (Br ) > 0.

(21)

br − Br ∥op < ηr /2. Applying Lemma C.4 with A = Br , Thus, for all sufficiently large n, ∥B b b A = Br , and d = q, and then using (20) gives br − Br ∥op 2∥B b r − Qr Q ≤ → 0. ηr op 30

For all sufficiently large n, the empirical and population AnchorPCA∞ projectors can be written as b ⌣,∞ = Π k

r−1 X

,∞ = Π⌣ k

bℓ + Q br , Π

ℓ=1

r−1 X

Πℓ + Q r .

ℓ=1

Hence b ⌣,∞ − Π⌣,∞ Π k k

≤ op

r−1 X

b ℓ − Πℓ Π

ℓ=1

b r − Qr + Q

→ 0. op

All displayed convergences hold for the fixed realization in Ω1 . Since P(Ω1 ) = 1, it follows that b ⌣,∞ − Π⌣,∞ Π → 0 almost surely. k

Corollary C.11 (Consistency of the empirical invariant block). Assume the conditions of Thm. C.9 (ii) b 1 denote the first empirical block returned by AnchorPCA∞ , i.e., and suppose that S⋆ ̸= {0}. Let Π d the projector onto Eig1 Then b 1 − ΠS Π ⋆

→0

almost surely.

d1 consistently estimates S⋆ . Equivalently, Eig Proof. The eigenvalue-one block of Π is S⋆ , so Π1 = ΠS⋆ . Step 2 in the proof of Thm. C.9 (ii) gives b ℓ − Πℓ Π → 0 almost surely for every ℓ; the case ℓ = 1 is the claim. op

Details on experiments

This appendix provides additional details for the motivating example in § 1.2, the synthetic S⋆ recovery simulations in § 6.1, a finite-penalty threshold simulation illustrating Prop. B.3, and the gassensor drift application in § 6.2. Code to reproduce all simulations and figures is publicly available at https://github.com/benesei/anchor-pca. D.1

Details of the motivating example

This appendix records the explicit construction used in the motivating example (§ 1.2) which is also used for the perturbation-path plot Fig. 2. Construction. Recalling § 1.2, let c1 , . . . , c4 be the canonical basis vectors of R4 , and set a := c1 , b := c2 , u := c3 . We further define two directions in the c3 –c4 plane, v := cos(50◦ )c3 + sin(50◦ )c4 ,

w := cos(100◦ )c3 + sin(100◦ )c4 ,

and let u⊥ , v ⊥ , w⊥ denote the unit directions in the c3 –c4 plane that are orthogonal to u, v, w, respectively. The population covariance matrices of the three domains are then defined as Σ1 = 220 uu⊤ + 140 aa⊤ + 90 bb⊤ + 25 u⊥ u⊥⊤ , Σ2 = 120 aa⊤ + 90 vv ⊤ + 70 bb⊤ + 10 v ⊥ v ⊥⊤ , Σ3 = 320 ww⊤ + 120 bb⊤ + 80 aa⊤ + 10 w⊥ w⊥⊤ ; these are shown in Table 1. Recovered subspaces. The solution of poolPCA spans span(a, c3 , c4 ) and therefore drops b. Both AnchorPCA25 and AnchorPCA∞ recover span(a, b). Their difference lies only in the third direction: AnchorPCA∞ selects v, and AnchorPCA25 selects q = 0.246c3 + 0.969c4 , a direction capturing slightly more variance on average than v (126.47 vs. 111.22), but is less invariant (q ⊤ Πq = 0.57 < 0.61 = v ⊤ Πv). 31

D.2

Details on the perturbation-path illustration

To visualize Thm. 4.1, we reuse the population motivating example and plot, for several methods (fit on the population source covariances), the worst-case average reconstruction error over Cλ as the perturbation strength ρ := 2Eλ increases. By Eq. (4), for any rank-k projector Π, this worst-case error is E E 2 1 X ρ X (e) , Π − Πk Tr (Ip − Π)Σe + E e=1 2E e=1 F so each method traces a line in ρ: the intercept is its average reconstruction error on the source covariances, and the slope is its disagreement with the local top-k projectors. The oracle minimizer at perturbation strength ρ is AnchorPCAρ/(2E) . We plot these lines for poolPCA, AnchorPCA25 , AnchorPCA∞ , and the wcPCA baselines maxRCS, maxRegret, and norm-maxRegret from [15]. Fig. 2 is an unsupervised analog of the perturbation-path plot in Anchor Regression [33]: poolPCA is optimal only near ρ = 0, a fixed finite-penalty AnchorPCAλ is optimal over a broad middle range of perturbations centered around its penalty, and AnchorPCA∞ is the right choice once robustness to top-k variance inflation matters more than reconstruction accuracy. The figure additionally shows that, in this example, the wcPCA baselines norm-maxRegret and maxRCS perform best only on short transition windows near the origin, though they target a different form of distribution shift: shifts within the convex hull of the source covariances. D.3

Details on the simulations

We now describe the details of the experiment described in § 6.1. D.3.1

Data-generating process

Notation. We write E = |E| for the number of domains, p for the number of covariates. Each domain e ∈ E is associated with a covariance Σe , whose construction we detail in the following paragraphs. Let k be the target rank for the low-dimensional representation, and m = dim(S⋆ ) be the dimension of the invariant subspace S⋆ . Let q := p − k, d := p − m. So, q is the dimension of the complement of each top-k subspace, and d is the dimension of S⋆⊥ . Throughout this section, a configuration corresponds to one choice of (E, p, k, m); for a given configuration and eigenvalue regime (see § D.3.1 below), a distribution draw is a multi-domain distributional instance, namely the full collection of covariances (Σe )e∈E with its own S⋆ (if we talk about multiple distribution draws, we assume they are chosen i.i.d.); a sample is one finite dataset drawn from a fixed distribution draw; and individual variables Xe,i are called observations or data points. Generating population covariance matrices. Each population covariance is constructed so that the top-k eigenspaces of the E covariances have a nontrivial maximal intersection of dimension m. We build this in three steps. 1. Invariant subspace. We draw a Haar orthogonal matrix U = (u1 , . . . , up ) ∈ Rp×p and split its columns into S := (u1 , . . . , um ), R := (um+1 , . . . , up ), so that S is an orthonormal basis of S⋆ := Im(S) and R is an orthonormal basis of S⋆⊥ . The remaining construction works in R-coordinates inside S⋆⊥ . 2. Bottom subspace per domain. For each domain e ∈ E, we draw a bottom subspace of dimension q inside S⋆⊥ : sample a d × q matrix with i.i.d. standard Gaussian entries, take the Q factor of its reduced QR decomposition to obtain an orthonormal matrix He ∈ Rd×q , and set Be := RHe . If the resulting collection (B1 , . . . , BE ) does not jointly span S⋆⊥ , we discard this collection and redraw all E bottom subspaces until dim span(B1 , . . . , BE ) = p − m, which is what makes the intersection of the top-k subspaces equal to S⋆ exactly, rather than something larger. 32

3. Top subspace per domain. The remaining k − m top directions in domain e ∈ E lie in the complement of Be inside S⋆⊥ . Let Ge ∈ Rd×(k−m) be an orthonormal basis of ker(He⊤ ) and set Ce := RGe . The top-k eigenspace for all e ∈ E is then

By construction,

Ue := Im(S, Ce ) = S⋆ ⊕ Im(Ce ). ⊥ T S ⊥ = (S⋆⊥ )⊥ = S⋆ . e∈E Ue = e∈E Im(Be ) = e∈E Im(Be )

For all e ∈ E, given the orthonormal basis (S, Ce , Be ), we define the covariance matrix as Σe := S diag(θe,⋆ )S ⊤ + Ce diag(θe,env )Ce⊤ + Be diag(ηe )Be⊤ , where θe,⋆ ∈ Rm are the eigenvalues on the invariant directions, θe,env ∈ Rk−m on the domainspecific top directions, and ηe ∈ Rq on the bottom directions. The eigenvalue regimes below ensure a local k versus k + 1 eigengap in every domain. Easy and hard eigenvalue regimes. We use the same subspace construction with two different choices of eigenvalues. In the easy regime, the invariant and domain-specific top eigenvalues are well separated from the bottom eigenvalues: θe,⋆,j ∼ U [5, 8],

θe,env,j ∼ U [5, 8],

ηe,j ∼ U [0.5, 3].

In the hard regime, the invariant directions still sit in each top-k subspace but carry less variance than the domain-specific top directions: θe,⋆,j ∼ U [3.2, 4.2],

θe,env,j ∼ U [6, 9],

ηe,j ∼ U [0.5, 2.5].

All distribution draws are generated with master seed 42. For each configuration (p, k, m, E) and each eigenvalue regime, we generate 100 independent distribution draws; S⋆ is redrawn for every replicate. Gaussian sampling. independently,

Given a distribution draw and sample size N , we sample N data points

Xe,i ∼ N (0, Σe ), i ∈ {1, . . . , N }, e ∈ E. For each distribution draw and sample size N , we generate 20 independent samples, each consisting of N observations per domain. The standard sample-size grid is N ∈ {50, 100, 200, 500, 1000, 2000, 5000}. For the small-E configuration described below, we extend the grid to include {10000, 30000, 100000}. Non-Gaussian covariance-preserving sampling. As a robustness check, we also use a symmetric two-component Gaussian mixture with the same mean and covariance as thepGaussian model. For each domain, let re be the leading eigenvector of Σe and set µe = 0.75 λmax (Σe ) re ; if needed, this vector is shrunk to ensure that Σe − µe µ⊤ e is positive semidefinite. We then sample for i ∈ {1, . . . , N }, e ∈ E, Xe,i = Ze,i + ξe,i µe ,

Ze,i ∼ N (0, Σe − µe µ⊤ e ),

P(ξe,i = 1) = P(ξe,i = −1) = 1/2.

Thus EXe,i = 0, Cov(Xe,i ) = Σe . The mixture preserves the population eigenspaces but is not Gaussian. D.3.2

Estimators

In the experiments, we consider AnchorPCA∞ and FindS⋆ . AnchorPCA∞ is fit with the default block tolerance of the package2 . tolN := min{0.05, 0.5 n−0.4 min }. In these experiments, each domain has the same sample size, so nmin = N . The first empirical block contains all eigenvalues ρbj satisfying b Thus the oracle block corresponding ρb1 − ρbj ≤ tolN , where ρbi is the ith largest eigenvalue of Π. to S⋆ is separated whenever ρb1 − ρbm ≤ tolN < ρb1 − ρbm+1 . FindS⋆ is used with significance level α = 0.05. 2 https://github.com/benesei/anchor-pca

D.3.3

Recovery configurations

We consider three configurations (E, p, k, m). 1. Main random-subspace configuration: (E, p, k, m) = (5, 10, 5, 2), shown in the main text. 2. Small-m configuration: (E, p, k, m) = (5, 8, 5, 1) shown in Fig. 6. 3. Small-E configuration: (E, p, k, m) = (2, 8, 5, 2) shown in Fig. 7. This configuration uses the minimal feasible invariant dimension. Indeed, since each domain has a q = p − k = 3 dimensional bottom space and the bottom spaces must span S⋆⊥ , feasibility requires p − m ≤ Eq, D.3.4

equivalently

m ≥ p − E(p − k) = 2.

Recovering the invariant subspace S⋆

We now repeat the experiments of Fig. 3 for non-Gaussian data, small-m, and small-E configurations. Performance quantities and aggregation. For each distribution draw and sample size N , we compute two quantities over the 20 independent samples: the dimension-recovery indicator, equal to 1 if the estimated first-block dimension equals the true m, and the operator-norm projector error for the estimated invariant subspace. Within each distribution draw, we average the 20 dimensionrecovery indicators and take the median of the 20 projector errors. The displayed dimension-recovery probability is the mean of these distribution-draw-level averages across distribution draws, while the displayed subspace error is the median of the distribution-draw-level median projector errors. Shaded bands, when shown, indicate the 10%–90% quantiles across distribution draws. Agreement-separation gap of the configurations. separation gap as

For a distribution draw, we define the agreement-

γΠ := ρm (Π) − ρm+1 (Π) = ρ1 (Π) − ρm+1 (Π) = 1 − ρm+1 (Π), where ρi (Π) denotes the ith largest eigenvalue of Π and the second and third equality use ρ1 (Π) = · · · = ρm (Π) = 1 on S⋆ . The gap γΠ measures how cleanly S⋆ separates from the rest of the agreement spectrum and reflects one aspect of difficulty of the configuration. Pooling the 100 easy and 100 hard distribution draws for each configuration, the median γΠ values, with 10–90% ranges, are median γΠ (10%–90%) configuration main random-subspace 0.277 (0.212–0.343) small-m 0.119 (0.072–0.184) small-E 9.8 · 10−3 (1.8 · 10−4 –5.5 · 10−2 ) The gap decreases across the three configurations, illustrating increased hardness. For the small-E configuration, the agreement gap is often very small, so this configuration is best viewed as a stress test rather than a typical regime. b − ρm+1 (Π). b This is the quantity relevant At the empirical level, we analogously track γ bΠ := ρ1 (Π) for the block-stabilized estimator: the first empirical block has correct dimension m if b − ρm (Π) b ≤ tolN < ρ1 (Π) b − ρm+1 (Π). b ρ1 (Π) Non-Gaussianity. Fig. 5 repeats the main-text experiment from Fig. 3 (E = 5, p = 10, k = 5, m = 2), but draws data from the covariance-preserving Gaussian mixtures introduced in § D.3.1 instead of Gaussians. The recovery curves are close to the Gaussian ones in Fig. 3. This is in line with § 5.3, which does not require Gaussianity for AnchorPCA∞ under its moment assumptions. The dashed curves show performance of FindS⋆ and are included to test its robustness to violations of its Gaussian data assumption. Small-m configuration. Fig. 6 shows Gaussian results for the small-m configuration (E = 5, p = 8, k = 5, m = 1). This configuration has a smaller population agreement-separation gap than the main configuration, but the finite-sample recovery curves are comparable on the standard grid. In both eigenvalue regimes, AnchorPCA∞ reaches correct first-block dimension recovery by the largest standard-grid sample sizes. 34

Dimension recovery

1.0 c = m) Pr(m

0.6 0.4

0.8

AnchorPCA ∞

0.6

FindS

color = DGP easy hard

0.4

0.2 0.0

texture = method

kΠ S − Π S k op

0.8

Subspace error

1.0

0.2

102

103

0.0

102

N per environment

103

N per environment

Figure 5: Main random-subspace configuration with Gaussian mixtures. Same distribution draws and aggregation as in Fig. 3, but observations are drawn from a symmetric two-component Gaussian mixture preserving each covariance Σe . The panels show the fraction of samples for which m was estimated correctly and median projector error. Solid lines use the first grouped eigenspace from AnchorPCA∞ , while dashed lines show FindS⋆ using significance level 0.05. Shaded bands are 10–90% quantiles across distribution draws and are shown for AnchorPCA∞ . As suggested by the consistency results (§ 5.3), for AnchorPCA∞ the estimated recovery probability converges to one and the subspace error to zero and, for FindS⋆ , the estimated recovery probability approaches 0.95. The close agreement with the Gaussian curves in Fig. 3 suggests that, in this configuration, the AnchorPCA∞ block-recovery mechanism and FindS⋆ are robust to this covariance-preserving non-Gaussian sampling perturbation.

Dimension recovery

1.0 c = m) Pr(m

0.6 0.4

0.8

AnchorPCA ∞

0.6

FindS

color = DGP easy hard

0.4

0.2 0.0

texture = method

kΠ S − Π S k op

0.8

Subspace error

1.0

0.2

102

103

0.0

102

N per environment

103

N per environment

Figure 6: Small-m configuration. Plots show the same as Fig. 5 for the small-m configuration. Again, the results are as suggested by § 5.3: e.g., in the left plot, the errors of AnchorPCA∞ and FindS⋆ approach zero and 0.05, respectively.

Small-E configuration. For the small-E configuration (E = 2, p = 8, k = 5, m = 2), we use the extended sample-size grid up to N = 100000. This configuration has a much smaller population agreement-separation gap than the other configurations, so we expect it to be a particularly difficult case. Fig. 7 shows the same quantities as Fig. 5 but also shows the tolerance diagnostic. In the top row, the first-block dimension is correct at N = 100000 in about 0.68 of easy-regime runs and 0.57 of hard-regime runs. FindS⋆ improves dimension recovery in the easy regime, reaching about 0.96 correct recovery, but in the hard regime it tends to overestimate the dimension, with median m b = 3. Intuitively, this is indeed a difficult case for FindS⋆ because it estimates m indirectly, by (k) (k) testing r = dim(span{Ip − Π1 , Ip − Π2 }) = dim(S⋆⊥ ) = p − m for all p − k to p − mmin ; (k) see § B.3. Thus, since p − m = 6 and rank(Ip − Πe ) = 3, this configuration has the smallest possible m = mmin = 2; equivalently, r = p − m is maximal, so the correct dimension m = 2 is reached only if the test rejects all smaller candidate dimensions r = 3, 4, 5, leading to a small empirical power of the test. If sampling noise makes it fail to reject r = 5, it stops there and returns m b = p − 5 = 3. AnchorPCA∞ uses a different criterion: it thresholds the empirical eigenvalues directly using tolN , so it can still return m = 2 whenever the first two eigenvalues form a separated 35

Dimension recovery

1.0

Subspace error

c = m) Pr(m

0.6

0.8

AnchorPCA ∞

0.6

FindS

0.4

color = DGP easy hard

0.4 0.2 0.0

texture = method

kΠ S − Π S k op

0.8

1.0

0.2

102

103

104

N per environment

105

0.0

102

agreement-separation gap

10 1

10 2

10 3

10 3 103

N per environment

104

105

Small-E configuration

Main random-subspace configuration

102

103

N per environment

block tolerance tol N

102

103

104

N per environment

105

γbΠ (easy regime) γbΠ (hard regime)

Figure 7: Small-E configuration and agreement-separation gap. The top row shows the same as Fig. 5 for the small-E configuration. As expected (see the discussion around γΠ ), this is a particularly hard case; indeed, the error of AnchorPCA∞ decreases but this time more slowly. Even for large N , some of the distribution draws are difficult in that for none of the 20 samples the estimated m is correct, which results in wide shaded areas. The bottom row (log scale) compares the automatic block tolerance tolN (dark blue) with the median γ bΠ in the easy and hard regimes for the main random-subspace configuration (left) and the small-E configuration (right). The plot shows that the small-E setting is difficult because γ bΠ is much smaller than in the main configuration and, for many of the considered N , is larger than tolN .

block, even when sampling noise makes S⋆⊥ appear five dimensional. The bottom row of Fig. 7 illustrates why the small-E configuration is also hard for AnchorPCA∞ : γ bΠ is much smaller than in the main configuration and is smaller than the automatic block tolerance over much of the extended grid. We hypothesize that the performance of AnchorPCA∞ could be improved by a more involved choice of the block tolerance. D.4

Gas sensor array drift details

Data and domains. The UCI gas-sensor drift archive contains 13,910 measurements from 16 metal-oxide chemical sensors exposed to six gases at multiple concentration levels [13, 41], see Fig. 8. Each measurement comes from a 16-channel metal-oxide gas sensor array and is summarized by eight features per channel, giving p = 128 features in total. Two are steady-state features: the maximal resistance change DR with respect to the baseline, and a normalized version of DR, obtained by dividing DR by the sensor resistance measured under gas exposure. The remaining six are transient features computed from rising and decaying exponential moving averages of the signal. The archive is organized into ten temporal batches: B1 contains months 1–2, B2 months 3,4,8,9,10, B3 months 11–13, B4 months 14–15, B5 month 16, B6 months 17–20, B7 month 21, B8 months 22–23, B9 months 24 and 30, and B10 month 36. We use these batches as domains. Data are available from the UCI Machine Learning Repository [41] at https://doi.org/10.24432/C5MK6M and are licensed under CC BY 4.0. 36

Gas composition (%)

100 80

Ethanol Ethylene Ammonia Acetaldehyde Acetone Toluene

60 40 20 0 nBi

B1 445

B2 1244

B3 1586

B4 161

B5 197

B6 2300

B7 3613

B8 294

B9 470

B10 3600

Figure 8: Gas composition by temporal batch. Stacked bars show the percentage of each gas class in batches B1–B10; numbers below the bars give the batch sizes nBi . The batches are strongly imbalanced and their gas mixtures vary over time, with some gases absent in several early batches and B10 balanced across the six gases. Baselines and hyperparameters. We consider the following baselines. First, we fit PCA on the averaged source covariances and term this poolPCA. We solve AnchorPCAλ for λ = 1 and AnchorPCA∞ (using the block-stabilized plug-in estimator from § 5.1 with default tolerance) Lastly, we consider a worst-case PCA variance, norm-maxRegret [15] which is fit with 10 restarts, 2000 iterations, learning rate 0.01, and deterministic seed 12020 + 10000s + 101k. Protocol. For each last source batch s ∈ {3, . . . , 8}, and for each k ∈ {10, 20, 30}, we fit the rank-k methods on B1–Bs and evaluate them on both the source batches and the held-out target batches B(s + 1)–B10. For each split, feature standardization is applied once using the pooled source observations from B1–Bs only, and the resulting transformation is applied unchanged to all source b b is computed from that batch’s and target batches. After this preprocessing, each batch covariance Σ observations centered by their own empirical batch mean. For a representation W ∈ Op×k and batch b b , the reported percentage of explained variance (EV) is covariance Σ EVb (W ) := 100

b bW ) Tr(W ⊤ Σ . b b) Tr(Σ

Rolling temporal splits. Fig. 9 repeats the source–target evaluation of Fig. 4 across all source– target splits s ∈ {3, . . . , 8}. The top row shows mean source EV, the bottom row mean target EV; shaded bands are min–max ranges across batches (shown only for poolPCA and AnchorPCA∞ ). On source batches, poolPCA is best as it directly maximizes source average variance. On target batches, AnchorPCA∞ usually outperforms poolPCA and AnchorPCAλ=1 , with the largest gains at k = 20 and k = 30 and competitive performance at k = 10. Only for s = 8, poolPCA seems to perform better but this split is noisier because its target set contains only B9 and B10. The overall pattern matches the main B1–B6/B7–B10 result in the main text. Explained variance of estimated S⋆ . The AnchorPCA∞ solution has rank k, but its first grouped eigenspace is a finite-sample estimate of the invariant subspace S⋆ with dimension m b (Thm. C.9). For s ∈ {3, . . . , 8} and k ∈ {10, 20, 30}, Fig. 10 compares the percentage of explained variance contributed by this estimated block with source-only poolPCA using the same dimension m b on B9 and B10. The comparison is evaluated separately on target batches B9 and B10. Across s = 3, . . . , 8, the estimated dimensions are m b = (3, 3, 3, 2, 2, 2) for k = 10, m b = (5, 5, 5, 5, 5, 5) for k = 20, and m b = (9, 8, 7, 7, 7, 6) for k = 30. On the target batch B9, the estimated block explains more variance than the same-dimensional poolPCA baseline for k = 20 and k = 30. On B10, both subspaces explain high variance and are often close. 37

k=10

% explained variance (source batches)

100

k=20

k=30

90 80 70 60

% explained variance (target batches)

100 80 60 40 20

Last source batch s

AnchorPCA λ = 1

poolPCA

Last source batch s

AnchorPCA ∞

Figure 9: Rolling source–target splits on gas-sensor drift data. Each column fixes k; the horizontal axis is the last source batch s. Lines show mean explained variance over the source batches (top row) or held-out target batches (bottom row). Shaded bands show min–max ranges over batches for poolPCA and AnchorPCA∞ .

k=10

B9 % explained variance

k=20

k=30

30 20 10 0

B10 % explained variance

100 90 80 70 60

m̂ = 3

m̂ = 2

m̂ = 5

m̂ = 9

m̂ = 8

m̂ = 7

m̂ = 6

Last source batch s

poolPCA top m̂

Last source batch s

estimated S

Figure 10: Estimated invariant block vs. same-dimensional poolPCA. For each (s, k), the red curve evaluates the first grouped eigenspace (estimated top eigenspace of Π, equivalently an estimate of S⋆ ) selected by AnchorPCA∞ ; the blue curve evaluates the top m b source-only poolPCA directions, where m b is the estimated dimension of S⋆ . Values are shown separately for target batches B9 and B10. 38

Nonlinear extensions

In this appendix, we discuss nonlinear extensions of our approach to finding rank-k robust or invariant subspaces via Anchor PCA translated to autoencoder-based approaches. Autoencoders (AEs) are a popular approach for nonlinear dimension reduction [18, 34]. As mentioned in § 1.1, the projection matrix ΠW = W W ⊤ can be interpreted as a linear autoencoder which first projects a given x to a k-dimensional embedding via W ⊤ ∈ Rk×p and then reconstructs x by mapping back to the original p-dimensional space via W ∈ Rp×k . An autoencoder functions analogously, but with possibly nonlinear encoder and decoder, typically parameterized as (deep) neural networks and optimized via stochastic gradient descent [31]. E.1

Vanilla autoencoder (AE)

An autoencoder (f, g) consists of an encoder f : Rp → Rk , which produces k-dimensional embeddings z = f (x), and a decoder g : Rk → Rp , which produces reconstructions x̂ = g(z) = g(f (x)). A standard autoencoder is trained to minimize the mean squared reconstruction error h i 2 LAE (f, g; P) = Ex∼P ∥x − g(f (x))∥2 , where, in practice, P is replaced by its empirical distribution. Given data from multiple domains, one approach to learning a shared nonlinear k-dimensional embedding is to train an autoencoder on the pooled data, similar to poolPCA in the linear setting. We refer to this approach as poolAE. However, like poolPCA, this approach is susceptible to focusing on spurious directions that exhibit high variation in only a few domains. E.2

Anchor autoencoder (AnchorAE)

We now discuss a method we call AnchorAE as a more robust alternative. Let {(fe , ge )}E e=1 be domain-specific AEs, i.e., (fe , ge ) ∈ arg min LAE (f, g; Pe ) . f,g

In analogy to Eq. (2), we define AnchorAEλ as the following unconstrained optimization problem E h i X 2 min Ex∼P̄ ∥x − g(f (x))∥2 + λ R (g ◦ f, ge ◦ fe ) f,g

(22)

e=1

PE where P̄ = E1 e=1 Pe denotes the (empirical) pooled distribution and R is a suitable invariance penalty between the shared and domain-specific AEs (discussed in more detail below). The first term in Eq. (22) is the pooled reconstruction error, which can be viewed as the AE-equivalent or proxy of explained variance in PCA. The second term is meant to capture disagreement among the nonlinear projections produced by the shared g ◦ f and domain-specific ge ◦ fe AEs, similar to the disagreement term (e)

ΠW − Π k

between projection matrices in AnchorPCAλ in Eq. (2). On the choice of invariance penalty. Whereas in the linear case, it is possible to compare the projection matrices directly, it is less clear how to assess agreement between nonlinear functions, as these are infinite-dimensional objects whose global behavior is not determined by their local behavior. In principle, we can consider measuring invariance w.r.t. domain-specific nonlinear projections through suitable norms ∥(g ◦ f ) − (ge ◦ fe )∥ in function space. However, each AE is unconstrained outside the support of the (empirical) distribution used to train it, i.e., Pe for (fe , ge ) and P̄ for (f, g). Hence, we should evaluate disagreement only on the respective shared supports. We therefore propose the following invariance penalty h i 2 R (g ◦ f, ge ◦ fe ) = Ex∼Pe ∥g(f (x)) − ge (fe (x))∥2 , 39

which can be viewed as the regret incurred by using the shared instead of the domain-specific AE for reconstructing data from domain e [2, 14]. Limitations. AnchorAE requires training E + 1 separate AEs, since to the best of our knowledge — and unlike in the linear case—the problem cannot be reduced to a modified single learning objective for (f, g). This approach may thus be computationally expensive, especially when the number of domains is large. In the following subsections, we discuss an approach for training a single, distributional AE from multi-domain data. E.3

Distributional principal autoencoder (DPA)

The distributional principal autoencoder (DPA) [38] is a nonlinear dimension reduction approach, which can be interpreted as a nonlinear variant of PCA that, in contrast to standard autoencoders, emphasizes distributional rather than mean reconstruction. We first review DPA in the single-domain setting and then discuss a possible extension for our multi-domain setting. Given a number k ≤ p of components, DPA consists of an encoder f : Rp → Rk , which produces embeddings z = f (x), and a stochastic decoder or generator g : Rk ×Rp−k → Rp , which, in addition to z, takes (p − k)-dimensional noise ϵ ∼ N (0, I) as input and allows for sampling reconstructions x̂ = g(z, ϵ). The DPA objective is constructed such that, for a fixed encoder f , the optimal DPA decoder g ∗ maps a given embedding z to the distribution of X given f (X) = z, d

g ∗ (z, ϵ; f ) = (X|f (X) = z) ,

(23)

where = denotes equality in distribution. In other words, the optimal decoder evaluated at z matches the distribution of realizations of X that are mapped by the encoder to z. At the same time, the DPA encoder minimizes the variability in the distributions in (23) by encoding the first k ‘principal’ components. Both of these goals are achieved by minimizing the following DPA objective LβDPA (f, g; P) 1 β β ∥g(f (X), ϵ) − g(f (X), ϵ′ )∥2 , = EX∼P,ϵ∼N (0,I) ∥X − g(f (X), ϵ)∥2 − E iid 2 X∼P,ϵ,ϵ′ ∼ N (0,I) where β ∈ (0, 2). This objective corresponds to the negative expected energy score [16] between X and the corresponding stochastic decoder output, conditional on the encoding f (X) of X. Due to this conditioning, the DPA objective differs from an energy distance by a normalization constant which depends on the encoder and encourages capturing principal (i.e., variation-minimizing) components. E.4

Multi-domain DPA

One could extend the DPA approach to our task of learning invariant subspaces from multi-domain data as follows. To map data from different distributions {Pe }e∈E to an invariant subspace, we consider a single shared encoder f : Rp → Rk , whereas to allow for accurate distributional reconstruction across domains, we consider domain-specific decoders ge : Rk × Rp−k → Rp e∈E . We could now train E domain-specific DPAs. However, in contrast to AnchorPCAλ and AnchorAE, it is unclear how to then compare or penalize disagreement between a shared DPA and domain-specific ones since, by construction, each DPA can fully reconstruct any distribution if given sufficient capacity in the form of additional p − k-dimensional noise as input to the stochastic decoder. Moreover, as discussed above, training a separate model for each domain is computationally demanding. This motivates considering the following objective for training a single multi-domain DPA: X β LβMD - DPA f, {ge }e∈E ; {Pe }e∈E := LDPA (f, ge ; Pe ).

(24)

e∈E

Since domain-specific decoders are used, this is different from a distributional, nonlinear version of poolPCA. Intuitively, the additional reconstruction capacity from domain-specific decoders and noise may enable the encoder to focus on shared directions, rather than explaining pooled variance.

Record · ID 259500 · SHA-256 6d3f2357d5142dcf
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.