arXiv CSopen access

Online Supervised Dimension Reduction with Random Features: Diagnostics and Computational Trade-offs

Online Supervised Dimension Reduction with Random Features:

arXiv:2609.20454v1 [stat.ML] 17 Sep 2026

Diagnostics and Computational Trade-offs Zhenlin Yao, Wei Xiong School of Statistics University of International Business and Economics Beijing, China

Abstract Accurate optimization of a supervised spectral objective need not produce an accurate population subspace or a better predictive representation. We investigate these distinctions for Online Kernel Supervised Principal Component Analysis (OKSPCA), which combines a centered cross-moment in finite random-feature coordinates with an Adam-style orthonormal basis update for an established objective. Fixed-map consistency, concentration and perturbation results describe the estimator and its exact subspace; same-target comparisons then assess the practical iterate separately. Across six predictive benchmarks, performance depends on the declared pipeline: replacing the tracker with the exact empirical target leaves the two regression deficits largely unchanged. Direct classification-rank models capture nearly all terminal objective energy on average, but a saved intermediate state exhibits substantial geometric deviation; a controlled sample-size study further separates empirical accuracy from population recovery. In distinct numerical-service workloads, exact on-request computation is faster in the tested classification settings, whereas Adam saves time relative to the tested full thin-SVD service for some dense wider-regression requests, alongside persistent geometric error. These diagnostics limit explanations based solely on terminal optimization accuracy and distinguish numerical cost from quality, rank coverage and freshness; they establish neither practical-tracker convergence nor predictive or deployment benefits from basis availability.

Keywords: supervised dimension reduction; random features; cross-covariance; subspace tracking; numerical basis access

1

Introduction

1.1

Motivation

Supervised dimension reduction seeks predictor coordinates that preserve response-associated information. Supervised PCA, HSIC-based supervised embeddings and covariance-based methods offer several ways to express this association [1–4]. In a stream, explicit features permit the associated empirical moments to be updated without retaining a growing sample-level Gram matrix. An 1

orthonormal basis can then be maintained approximately, computed when requested, or cached between factorizations. These choices raise two separate questions: how accurately the basis solves its empirical objective, and whether improving that solution changes the prediction or numerical service for which it is used. This paper studies those questions through a common empirical target. A same-target comparison holds the accumulated moment fixed and replaces the maintained subspace with its exact spectral solution. It can therefore test whether a particular discrepancy is explained by terminal tracking accuracy, without treating different supervised criteria as interchangeable. Prediction remains a property of the complete representation and evaluation pipeline. Likewise, a fast basis response can be stale, rank-limited, or geometrically inaccurate. The purpose is to determine which explanations the observed comparisons support, rather than to infer usefulness from objective capture alone.

1.2

Related work

The operator is established. In finite coordinates, HSIC-supervised PCA has a centered cross-product form [2, 5]; SRP and KSRP also use centered explicit-feature products for supervised embeddings [6]. The positive-spectrum input subspace agrees with that of orthonormal covariance-maximizing two-block PLS, although the respective objectives sum squared singular values and singular values. CCA additionally whitens marginal covariances [7–9]. Random Fourier features provide the nonlinear coordinates, not a new supervised objective [10]. Sequential spectral computation also has substantial precedents: recursive and incremental PLS [11–14], stochastic PLS [15–17], paired-stream singular-subspace tracking [18, 19], online supervised reduction [20], and online SIR and kernel SIR [21, 22]. Incremental SVD, online PCA and streaming kernel PCA supply related numerical mechanisms [23–27]. We use covariance-free incremental PLS (CIPLS) as a predictive comparator [28]; it is not an exhaustive representation of these alternatives. Existing KSPCA population-objective bounds likewise differ from a guarantee for the practical tracker considered here [29].

1.3

Study objectives and contributions

OKSPCA combines recursively centered moments with a Stiefel-constrained Adam-style update. We first define this actual computation and state its fixed-map statistical foundation (Section 2). Three comparisons then organize the evidence. First, can exact solution of the same terminal moment explain the predictive gaps in the declared pipelines? The two regression deficits largely survive this replacement. Second, do high terminal objectives describe the maintained geometry? Direct classification-rank reconstructions include a marked intermediate deviation, while a controlled sample-size study separates empirical and population subspaces (Sections 3 and 4). Third, what is the cost of supplying a basis under common request schedules? Exact on-request, projected gradient, periodic caches and applicable structured exact services expose conditional trade-offs that differ between the numerical workloads (Section 5).

2

The contribution is this set of controlled empirical-target diagnostics and its computational consequences, not a new HSIC/PLS objective or a general superiority claim. The predictive comparisons retain their adapter, coordinatewise evaluator and tuning differences. The service study varies several workload factors jointly and compares specified implementations, including full thin SVD; it does not isolate target width or establish competitiveness against all leading-k spectral solvers. These distinctions determine how far the findings can be generalized.

2

Methodology and theoretical properties

2.1

Supervised cross-covariance objective

The diagnostic comparisons require an explicit distinction between the operator being estimated, an exact basis for that operator, and the basis returned by an implementation. We first define their common objective, then give the causal update actually used in the experiments. The final two subsections explain what objective capture can say about geometry and which statistical guarantees require a different, fixed-map setting. Prediction is evaluated separately because an objective on cross-moments does not specify a predictive loss or its fitted evaluator. Let ϕ(x) ∈ RDx and ψ(y) ∈ RDy be finite feature coordinates. With fixed maps and preprocessing, define C∗ = E[(ϕ − µϕ )(ψ − µψ )⊤ ],

M∗ = C∗ C∗⊤ .

(1)

For 1 ≤ k ≤ Dx , the established supervised spectral objective is max

U ∈St(Dx ,k)

∥C∗⊤ U ∥2F ,

St(Dx , k) = {U : U ⊤ U = Ik }.

(2)

An exact maximizer spans leading left singular directions, with possible nonuniqueness at a tied boundary. For any moment C and 1 ≤ q ≤ rank(C), orthonormal covariance-maximizing two-block PLS selects the same input-side maximizing subspaces: its criterion sums the leading singular values, whereas Eq. (2) sums their squares [15, 16, 30, 31]. This relation concerns the optimization problems, not arbitrary PLS/CIPLS implementations or trajectories; CCA additionally whitens marginal covariances. For feature-column matrices Z, W , H = In − 11⊤ /n and L = W ⊤ W , the Bessel-normalized empirical moment satisfies Cn = (n − 1)−1 ZHW ⊤ ,

Cn Cn⊤ = (n − 1)−2 ZHLHZ ⊤

(n ≥ 2).

(3)

This is the finite-coordinate HSIC supervised-PCA operator [2]. Random-feature supervised embeddings, including SRP/KSRP centered cross-products, precede the present construction [6]. Our comparisons concern this common finite-feature target; neither the objective nor the centering identity is claimed as new. The analysis supplies no full-kernel approximation guarantee.

3

The cross-product used by a direct supervised embedding need not itself be an orthonormal leading-subspace basis. Conversely, agreement of optimal input-side subspaces does not require two methods to return identical coordinate systems. Scaling, rotation, feature construction and the downstream fitting rule can therefore remain relevant even when two methods are related through the same spectral operator. In OKSPCA, “kernel” denotes the stationary kernels approximated by the sampled random Fourier features (RFF), whose finite coordinates remain the evaluation target. Classification uses unweighted one-hot labels, Dy = C. For positive class probabilities pc , feature means µc within classes and overall mean µ, C∗ [:, c] = pc (µc − µ),

M∗ =

X

p2c (µc − µ)(µc − µ)⊤ .

(4)

c

Thus conditional feature-mean contrasts receive squared-probability weights, not the weights of ordinary between-class scatter. Zero-probability classes contribute zero columns. This criterion guarantees neither sufficient prediction nor macro-F1 optimality. Since centered labels sum to zero, C∗ 1C = 0,

rank(C∗ ), rank(C̄t ) ≤ C − 1.

C̄t 1C = 0,

(5)

The sample identity is pathwise even with evolving input preprocessing. The cap C − 1 is an upper bound, not an asserted rank or positive gap. If k > r = rank(C∗ ), exact maximizers include a generally nonunique (k − r)-dimensional null-space completion; at k = Dx the whole space is fixed. Rotated basis columns cannot individually be designated as supervised or completion coordinates.

2.2

Online moment estimation and basis updates

x , Sy At observation t, input and regression-target standardizers St−1 t−1 use only prior observations.

The input map is ϕt =

q

x 2/Dx cos(Wx St−1 (xt ) + bx ).

(6)

The entries of Wx ∈ RDx ×dx are sampled from N (0, σx−2 ), and phases from Unif(0, 2π), once per run. Regression uses the analogous map of its causally standardized scalar response; classification uses its declared one-hot vector. Starting from zero means and raw moment, update µϕ,t = (1 − t−1 )µϕ,t−1 + t−1 ϕt , µψ,t = (1 − t−1 )µψ,t−1 + t−1 ψt ,

(7)

Rt = (1 − t−1 )Rt−1 + t−1 ϕt ψt⊤ , C̄t = Rt − µϕ,t µ⊤ ψ,t . Induction and expansion give C̄t = t−1

t X

(ϕi − µϕ,t )(ψi − µψ,t )⊤ .

i=1

4

(8)

This exact identity retains the history of realized features. It is not the covariance obtained by replaying every observation through the final standardizer. Set C1 = 0 and Ct = tC̄t /(t − 1) thereafter. The positive Bessel rescaling preserves singular subspaces and normalized objective ratios, but is an unbiasedness correction only under the fixed-map assumptions below. The random weights are fixed, but raw-to-feature coordinates evolve with the standardizer. Prior-only preprocessing is therefore part of the operator’s definition: exact factorization removes subspace-optimization error, not the effects of its feature history. For the current Ct , the Frobenius gradient is 2Ct (Ct⊤ U ). With sym(A) = (A + A⊤ )/2, its tangent projection and subsequent global clipping are Gt = 2Ct (Ct⊤ Ut−1 ), ⊤ GT,t = Gt − Ut−1 sym(Ut−1 Gt ),

e T,t = G

 GT,t ,

(9) ∥GT,t ∥F ≤ c,

T,t /(∥GT,t ∥F + 10

cG

−12 ),

∥GT,t ∥F > c.

The optimizer-call counter advances to s = st before bias correction. With elementwise operations, the actual update is e T,t , ms = β1 ms−1 + (1 − β1 )G e ⊙2 , vs = β2 vs−1 + (1 − β2 )G T,t b s = ms /(1 − β1s ), m

(10)

vbs = vs /(1 − β2s ),

p

b s ⊘ ( vbs + ϵ), ∆t = α m

Ut = qf + (Ut−1 + ∆t ).

Signed reduced QR, denoted qf + , makes the diagonal of its triangular factor nonnegative. The arrays remain in ambient coordinates without transport. Elementwise preconditioning need not preserve tangency or right-rotation equivariance, and a finite step need not increase the current objective. We therefore call this a Stiefel-constrained Adam-style update, not a Riemannian-Adam convergence construction [32–34]. The associative gradient multiplication uses the current plug-in moment; it is not a fresh unbiased observation of the population gradient. In particular, the quadratic operation C 7→ CC ⊤ does not preserve finite-sample unbiasedness of C. For Algorithm 1, write ϕx (z) =

p

2/Dx cos(Wx z + bx ) and ψy (z) =

q

2/Dy cos(Wy z + by ) for

the fixed RFF maps, and let κ(y) be the declared class index. The operator Clipc is the piecewise map in Eq. (9). In the Adam steps, squares, square roots and ⊘ act elementwise. The operation WF is a Welford state update with the small-count, variance-divisor and threshold conventions in Section S1; it does not transform the current observation again.

5

Algorithm 1 One-observation OKSPCA update: normal finite-data path Require: (xt , yt ), index t; fixed map parameters and class encoding; α, c, β1 , β2 , ϵ. x ; also S y State: Ut−1 , µϕ,t−1 , µψ,t−1 , Rt−1 , mst−1 , vst−1 , st−1 , St−1 t−1 for regression. D ×k x Initialization: U0 = qf + (Z0 ) for seeded Gaussian Z0 ∈ R ; feature means, raw moment, Adam arrays and counters zero; raw standardizers empty. A. Form features using prior standardizers x (x )) 1: ϕt ← ϕx (St−1 t ( eκ(yt ) , classification, 2: ψt ← y ψy (St−1 (yt )), regression. B. Accumulate the realized-feature cross-moment 3: µϕ,t ← (1 − t−1 )µϕ,t−1 + t−1 ϕt , µψ,t ← (1 − t−1 )µψ,t−1 + t−1 ψt 4: Rt ← (1 − t−1 )Rt−1 + t−1 ϕt ψt⊤ ( 0, t = 1, ⊤ 5: C̄t ← Rt − µϕ,t µψ,t , Ct ← tC̄t /(t − 1), t ≥ 2. C. Maintain the basis with ambient Adam 6: Gt ← 2Ct (Ct⊤ Ut−1 ) ⊤ G ) 7: GT,t ← Gt − Ut−1 sym(Ut−1 t e 8: GT,t ← Clipc (GT,t ), st ← st−1 + 1 e T,t 9: mst ← β1 mst−1 + (1 − β1 )G ⊙2 e 10: vst ← β2 vst−1 + (1 − β2 )GT,t mst vst b st ← 11: m vbst ← st , 1 − β1 p 1 − β2st b st ⊘ ( vbst + ϵ), Ut ← qf + (Ut−1 + ∆t ) 12: ∆t ← αm D. Advance raw standardizers last x ,x ) 13: Stx ← WF(St−1 t y y 14: St ← WF(St−1 , yt ) for regression only. Ensure: Updated carried state; Ut is the maintained basis, not necessarily a later query’s returned block.

6

On this normal path, st = t; the optimizer counter is nevertheless a distinct state variable. The first zero-gradient observation advances it and executes signed QR, whose zero diagonal receives sign +1. Settings and numerical-failure branches remain in Section S1; query-time extraction and checks are separate service operations. Evolving scalers, plug-in moments, clipping, fixed learning rate and ambient adaptive arrays fall outside the standard stochastic-retraction convergence assumptions; convergence of this coupled practical recursion is not established [35–37].

2.3

Objective values and subspace geometry

For U ∈ St(Dx , q), define Jq (U ; C) = ∥C ⊤ U ∥2F ,

Jq∗ (C) =

q X

σj (C)2 ,

ρq = Jq /Jq∗ .

(11)

j=1

Singular spectra are extended by zeros; the ratio requires Jq∗ > 0. An exact reference is a leading subspace of the same moment. To extract the best q ≤ k block within a maintained U ∈ St(Dx , k), put BU = U ⊤ C̄t . The Ky Fan principle gives max Jq (U Q; C̄t ) =

Q∈St(k,q)

q X

σj (BU )2 .

(12)

j=1

Leading left singular vectors of BU give the coefficient block, followed by signed QR of U Q. The value uses the projected spectrum, not the unrestricted moment spectrum; uniqueness can fail even at a zero tie. Extracting this block does not refit the tracker at rank q. The prediction and service comparisons retain their different training, comparison-rank and returned-rank conventions in Sections 3 and 5. At equal rank, extraction may still rotate coordinates without changing the span, its objective or projector. The coordinate-standardized penalized probes in Section 3 can depend on that rotation. Lower-rank and cached query outputs are instead evaluated under the rank and freshness conventions in Section 5. (q)

For a fixed positive-semidefinite M , any rank-q projector P , and a leading rank-q projector P∗ , assume 1 ≤ q < Dx , Jq∗ =

j≤q λj (M ) > 0 and δq (M ) = λq (M ) − λq+1 (M ) > 0. Then

P

δq (M ) (q) ∥P − P∗ ∥2F , 2 s 2Jq∗ (1 − ρq ) (q) ∥P − P∗ ∥F ≤ . δq (M )

Jq∗ − tr(M P ) ≥

(13)

Here ρq = tr(M P )/Jq∗ . This empirical gap is distinct from a population gap or numerical-rank threshold; at a tied boundary the division bound is unavailable. A small relative gap permits high objective capture with substantial projector discrepancy. In the nonzero binary rank-one case,

7

(1)

ρ1 = cos2 θ and ∥P − P∗ ∥2F = 2(1 − ρ1 ), linking the mean squared distance, not the mean distance, to mean objective capture. Proofs are in Section S1.

2.4

Fixed-map theory and computational complexity

The following statement collects standard estimator, concentration and perturbation arguments with their separate conditions. It concerns any exact leading-k spectral projector Pbt of Mt = C̄t C̄t⊤ , and the leading-k population projector P∗ of M∗ . Finite second moments suffice for consistency. Almost-sure bounds on whole feature vectors supply the stronger concentration statement. Identification of a particular leading population subspace additionally needs a positive separator after the selected rank. These conditions serve different purposes: neither bounded RFF outputs nor the centered-label rank cap alone supplies the population gap. Likewise, the pathwise recursion above remains correct for dependent or evolving-map features even when the probabilistic conclusions do not apply. Proposition 1 (Conditional estimator and exact-subspace guarantees). Condition on sampled map parameters and fixed preprocessing constants. If the feature pairs are i.i.d. with finite second moments, then EC̄t = (t − 1)C∗ /t,

C̄t −→ C∗

almost surely.

(14)

If also ∥ϕt ∥2 ≤ Bx and ∥ψt ∥2 ≤ By almost surely, then for each fixed integer t ≥ 1 and δ ∈ (0, 1), with probability at least 1 − δ, C̄t − C∗

op

q √ ≤ bt (δ) = 6 2Bx By log(6/δ)/t.

(15)

Separately, if 1 ≤ k < Dx and the population boundary gap ∆k = λk (M∗ ) − λk+1 (M∗ ) is positive, then deterministically √

2 2k ∥Pbt − P∗ ∥F ≤ ∥Mt − M∗ ∥op , ∆k

∥Mt − M∗ ∥op ≤ C̄t − C∗

op

{2 ∥C∗ ∥op + C̄t − C∗

(16) op

}.

The full proof in Section S1 retains the centering term, fixed-time event and population-gap conditions [38–40]. Combining the bounds gives a leading O(t−1/2 ) high-probability order for the exact projector at fixed dimensions, confidence and positive gap. RFF maps satisfy ∥ϕ∥22 ≤ 2, and similarly for regression targets; one-hot targets have norm one. The concentration bound has no explicit feature-width factor, although embedding bounds, selected rank and eigengap can depend on dimension. Bounded features neither make evolving preprocessing i.i.d. nor control full-kernel approximation. The bounded fixed-feature study in Section 4 provides a finite-grid diagnostic, not an estimated convergence rate. The three objects remain distinct throughout the empirical sections: the population projector belongs to a known fixed-feature distribution, the exact empirical projector belongs to its accumulated 8

sample moment, and the practical projector belongs to the implemented recursion. Distances between these objects compare different subspace pairs; they do not form an additive attribution of observed error. The fixed-time probability statement likewise does not certify every requested basis along an adaptive or arbitrarily long stream. With dimensions fixed, the core named arrays occupy O(Dx Dy + 3Dx k + Dx + Dy ) entries, plus dense maps, scalers and implementation buffers specified in Section S1. Per-observation work has order O(Dx dx + Ireg Dy dy + Dx Dy k + Dx k 2 + k 3 ),

(17)

where Ireg indicates regression. Lower-order scaling, one-hot, mean and adaptive-array operations are included in this bound for positive dimensions and k ≥ 1. Associative multiplication avoids a Dx × Dx matrix. These arithmetic and named-array bounds are not process-memory or measured runtime laws; Section 5 charges the work needed to deliver a requested basis.

3

Predictive performance and same-target diagnostics

3.1

Experimental setup

The predictive comparison asks whether a more accurate solution of the empirical target explains the observed differences between representation pipelines. Four classification tasks use held-out macro-F1, and two regression tasks use held-out R2 . All reduced methods are trained at qcap = C − 1 for classification and at k = 4 for regression (Table 1). The cap is an algebraic upper bound, not an estimated population rank or an optimal budget for every method. Ten final seeds, 100, . . . , 109, share data, splits, and stream permutations across methods. Synthetic seeds change the generated data; real-data seeds preserve held-out membership. Their results are not pooled into a common task ranking. HAR preserves its official subject-disjoint split; Sensorless uses a fixed stratified split, and Kin8nm preserves its existing held-out partition [41–43]. These are benchmarks rather than application Table 1: Predictive tasks and representation dimensions. Counts are final training/test sizes. The last column is the classification training cap qcap = C − 1 or regression rank four. Original budgets for the unreduced-control comparison are stated separately in Table S4.

Task

ntrain

ntest

dx

Dx

Dy

k

XOR-Gauss Spirals-3 Friedman-1 HAR Sensorless Kin8nm

4,000 4,000 4,000 7,352 46,807 6,553

1,000 1,000 1,000 2,947 11,702 1,639

20 20 10 561 48 8

64 64 128 1,024 256 128

2 3 128 6 11 128

1 2 4 5 10 4

9

case studies. Synthetic generators, development partitions, bandwidths, optimizer settings, and uncertainty conventions are specified in Section S2. After fitting the representation, prediction uses zU (x) = U ⊤ ϕ(STx (x)),

b U (SU (zU (x))) . ybU (x) = h

(18)

b U is logistic regression with inverse Here SU is fitted only to replayed training embeddings, and h

regularization one or Ridge with penalty one. Every representation receives its own standardizer and probe. Thus scores depend on coordinates, coordinatewise scaling, and penalized prediction, and are generally not rotation invariant. Terminal replay is retrospective: its fixed STx differs from the preceding scaler states that formed the causal moment. Test observations update neither representation nor probe. OKSPCA and RFF-CIPLS share input RFF coordinates. For multiclass tasks, the CIPLS comparators concatenate one-versus-rest reducers and Gaussian-compress their outputs to the requested dimension; the displayed OVR–RP names identify these adapters. Regression instead compares OKSPCA’s standardized RFF target with a scalar-target PLS recursion. CIPLS uses raw inputs, while CCIPCA and RFF-CCIPCA are unsupervised trackers. OKSPCA received development selection, whereas comparator settings were inherited without separate grids; the primary rank-cap fits were not retuned. These are complete-pipeline comparisons, not isolated interventions on optimization, target mapping, or supervised rank. Exact replacement holds the accumulated moment and its feature history fixed; it leaves possible mismatch with the final replay coordinates unresolved. It is a numerical reference, not a population oracle.

3.2

Predictive performance and exact-target comparisons

Figure 1 retains all four classification tasks and five reduced pipelines. OKSPCA has the highest mean macro-F1 on HAR and Sensorless. Its mean differences from RFF–OVR–RP–CIPLS are +0.0504 and +0.1513, respectively. On Spirals-3 its corresponding difference is +0.0170, but OVR– RP–CIPLS has the higher mean. On binary XOR-Gauss the difference from RFF-CIPLS is −0.0002, with a fixed-seed interval spanning zero. The paired intervals in Table S3 describe the declared seeds and do not establish a general multiclass advantage. On both regression tasks, OKSPCA scores below RFF-CIPLS, although above the two unsupervised trackers; linear CIPLS has the highest mean on Friedman-1 (Table 2). Solving the same terminal empirical target exactly leaves these deficits largely unchanged. The paired exact-minus-full mean R2 changes are +0.000141 for Friedman-1 and +0.000027 for Kin8nm. These are descriptive contrasts, not equivalence tests. They weaken the explanation that terminal tracking error alone causes the observed disadvantage, without identifying a unique alternative among feature mapping, target construction, output budget, and the evaluator.

10

Seed score

Mean

XOR-Gauss (q = 1)

Spirals-3 (q = 2) 0.40 Held-out macro-F1

Held-out macro-F1

0.64 0.56 0.48 0.40

0.36 0.32 0.28 0.24

0.32

CA

SP OK

S IPL

RF

F-C

S

PL

CI

A

PC

I CC

RF

F-

I CC

OKSPCA RFFOVRRPCIPLS

A PC

Sensorless (q = 10)

Held-out macro-F1

Held-out macro-F1

HAR (q = 5)

OVR- CCIPCA RFFRPCCIPCA CIPLS

0.84 0.80 0.76

0.80 0.72 0.64 0.56

0.72 OKSPCA RFFOVRRPCIPLS

OVR- CCIPCA RFFRPCCIPCA CIPLS

OKSPCA RFFOVRRPCIPLS

OVR- CCIPCA RFFRPCCIPCA CIPLS

Figure 1: Classification at the centered-label rank cap. Small points show all ten final-seed macro-F1 scores; diamonds show means. All five pipelines were trained from the beginning at qcap = C − 1 without retuning. Multiclass OVR–RP labels distinguish the adapters from binary scalar-target CIPLS.

11

Table 2: Regression at k = 4. The reduced-method benchmark reports ten-seed R2 means and sample standard deviations. The same-target comparison distinguishes the full maintained representation, its extracted block, and the exact empirical target; their absolute means accompany the descriptive paired exact-minus-full changes. Extracted-minus-exact intervals retain their separate identity in Table S3.

Reduced pipeline OKSPCA RFF-CIPLS CIPLS CCIPCA RFF-CCIPCA Task Friedman-1 Kin8nm

Friedman-1

Kin8nm

0.6630 ± 0.0417 0.7194 ± 0.0242 0.7224 ± 0.0161 0.3012 ± 0.1298 0.2893 ± 0.1150

0.4130 ± 0.0247 0.4546 ± 0.0236 0.4043 ± 0.0002 0.1828 ± 0.0979 0.1747 ± 0.0957

Full

Extracted

Exact

Exact − full

0.662962 0.413004

0.662963 0.413005

0.663103 0.413032

+0.000141 +0.000027

These regression endpoints originate from the larger-span terminal study, not from classificationrank reconstructions. Full, extracted, and exact representations have k = qT = 4, but extraction rotates coordinates before separate standardization and probe fitting. Their small observed differences do not make the evaluator rotation invariant. Section S3 defines this comparison and reports its geometry. Unreduced Raw-X and RFF-X controls also retain an important limitation: RFF-X has higher original-budget mean scores than OKSPCA on five tasks, except Spirals-3 (Table S4). Their different dimensions preclude reduced-method rankings. Exact-minus-full measures replacement of the benchmark representation; extracted-minus-exact compares its within-span maximizing block with the unrestricted empirical solution. Only the latter has the intervals in Table S3. Neither is a rotation-free subspace score or an equivalence test.

3.3

Primary-rank classification diagnostics

Forty primary-rank OKSPCA fits were reconstructed with the original inputs and settings, reproducing their terminal scores over all ten final seeds. Five retained checkpoints per trajectory supply 200 same-moment spectral comparisons. At termination, task-mean objective ratios range from 0.998824 to 0.999951; original and exact representations have similar mean scores under the stated probe (Table 3). This compares directly trained ranks, not extraction from an earlier larger-rank classification fit. Terminal agreement is not uniform along the observed path. Sensorless seed 104 reaches a maximum angle of 78.099◦ at observation 23,404, with ρ = 0.970687. Its five-checkpoint record is retained in Table S5; the finite checkpoints establish neither cause nor duration of the deviation. The reconstructed states are distinct from recovered historical states. Together, these diagnostics separate target capture, geometric agreement, and predictive usefulness while retaining the limitations of the original pipelines.

12

Table 3: Direct primary-rank classification fits and their own exact terminal targets. Entries are means over ten seeds. Original/exact columns are macro-F1; ρ, maximum principal angle, and unnormalized Frobenius projector distance compare the retained state with the same empirical moment. These are terminal summaries, not uniform pathwise guarantees.

Task

k

Original F1

Exact F1

ρ

θmax (deg)

dP

XOR-Gauss Spirals-3 HAR Sensorless

1 2 5 10

0.5960 0.3663 0.8528 0.7625

0.5962 0.3648 0.8527 0.7625

0.999745 0.998824 0.999539 0.999951

0.882 2.373 1.462 0.459

0.0218 0.0683 0.0448 0.0151

4

Empirical and population subspace accuracy

The preceding comparisons have an empirical reference but no known population projector. A bounded fixed-feature experiment distinguishes these objects under the conditions of Proposition 1. Its known cross-covariance C∗ has four positive squared singular values (0.49, 0.36, 0.25, 0.15); the selected subspace has rank three in dimensions (Dx , Dy ) = (32, 8). Twenty seeds use nested sample sizes 256, 512, . . . , 32,768, with fixed population orientations. Section S4 gives the generator and complete conventions. This setting does not reproduce the evolving preprocessing or random-feature approximation of the predictive benchmarks. Let P∗ , Pexact,n , and Ponline,n denote the population, exact empirical, and practical projectors. The three distances are dstat = ∥Pexact,n − P∗ ∥F ,

dopt = ∥Ponline,n − Pexact,n ∥F ,

dtotal = ∥Ponline,n − P∗ ∥F .

(19)

They compare different subspace pairs and are not an additive error decomposition; dtotal − dstat can have either sign. Objective ratios likewise depend on whether the reference operator is empirical or population. At n = 256, median sample-online objective capture is 0.995834, while median dstat is 1.416004 and the population objective ratios are about 0.731. Near-optimal empirical objective values therefore coexist with appreciable population error. Across the sample-size grid, median dstat decreases to 0.1281 and median dtotal decreases from 1.395 to 0.1292. At the final checkpoint, median dopt = 0.01128. The sample-online ratio is normalized by the empirical optimum; both population ratios use the known population optimum. Thus near-one sample capture and population error concern different references, not fractions of predictive information retained. The finite-grid contraction is compatible with the fixed-map analysis, but estimates no convergence rate and establishes no convergence of the practical update. The bands condition on fixed population orientations while varying streams and tracker initialization. This experiment cannot attribute the predictive deficits in Section 3 to population estimation. 13

Sample online / exact

dstat

Population online / truth

dopt

Population exact / truth

dtotal Frobenius distance (log scale)

Normalized objective ratio

A Objective capture 1.0 0.9 0.8 0.7 0.6 28

211 Sample size, n

215

B Three projector distances 100

10−1

10−2 28

211 Sample size, n

215

Figure 2: Objective capture and projector distances over eight nested sample sizes in the bounded generator. The objective panel distinguishes sample-online, population-online, and population-exact ratios; the distance panel retains dstat , dopt , and dtotal . Thin traces show all 20 seeds, emphasized lines and bands their medians and interquartile ranges. The distances are not additive.

14

5

Computational trade-offs in basis access

5.1

Service definitions and experimental setup

The service experiment asks what computation is required to deliver a basis for the current finitefeature moment. It has no predictive endpoint. Five synthetic workloads combine raw and feature dimensions, rank caps, and stream lengths as specified in Table 4; R128 and R256 name targetfeature widths, not raw response dimensions or predictive datasets. Five final seeds and five request schedules give 825 applicable method–regime–seed–schedule cells. All methods within a regime and seed process the same raw stream and causal frontend; dimensions vary jointly across regimes, so between-regime contrasts do not isolate target width. Table 4: Numerical-service workloads. Raw responses are scalar labels or regression values. All applicable services use final seeds 6200–6204 and request periods 1, 10, 100, 1000, or Bernoulli probability 0.01, with termination included. Periodic caches retain separate update periods and receive no forced terminal refresh.

Regime

dx

dy

Dx

Dy

k

n

C2 C4 C11 R128 R256

20 64 48 10 64

1 1 1 1 1

64 256 256 128 256

2 4 11 128 256

1 3 10 4 8

4096 4096 4096 4096 2048

The offline full-moment reference defines rref (t) = #{j : σj (C̄t )2 > 10−10 σ1 (C̄t )2 },

qref (t) = min{k, rref (t)}.

(20)

A zero leading singular value gives rref = qref = 0. An approximate service applies the same relative squared-spectrum rule to Ut⊤ C̄t , obtains rlocal , and returns qout = min(k, rlocal ) columns. It receives no offline full-spectrum rank; zero local rank does not establish a zero full moment. A request factors this small matrix and constructs its leading within-span block, including signed QR, copying, and checks. This paid query work leaves the maintained basis, moments, scalers, and optimizer unchanged; same-timestamp factors may be reused. Six services apply throughout: Adam-style maintenance, normalized projected gradient (PG), full thin SVD after each observation, full thin SVD on request, and exact caches updated every ten or one hundred observations. Binary direct extraction additionally applies to C2; contrast QR followed by a small SVD applies to C4 and C11. Generic exact services compute the full thin factorization before truncation. Structured services exploit centered-class identities, with floating-point structural residuals measured separately. A cache returns its last available basis and timestamp τ , of age t − τ ; before initialization it returns no basis.

15

For binary labels, C̄t = [ct , −ct ] gives C̄t C̄t⊤ = 2ct c⊤ t and exact direction ct /∥ct ∥2 when ct ̸= 0. More generally, an orthonormal Helmert contrast HC spans 1⊥ C . Centered labels imply C̄t HC HC⊤ C̄t⊤ = C̄t C̄t⊤ . Factoring the smaller contrast matrix therefore preserves the target. A plain thin QR alone identifies the maximizing span only with sufficient input dimension, full contrast rank, and the full requested contrast span. The evaluated contrast service instead takes an SVD of its triangular factor to detect rank and select the required block (Section S5). PG shares Adam’s seeded initial basis but uses Ht = 2Ct (Ct⊤ Ut−1 ), Ut = qf +

⊤ Ht ), HT,t = Ht − Ut−1 sym(Ut−1

γHT,t Ut−1 + 2∥C̄t ∥2F + 10−12

(21)

!

.

The gradient uses Bessel-normalized Ct , whereas the denominator uses C̄t . PG has no Adam arrays, clipping, or line search. An exactly zero tangent gives zero displacement followed by signed QR. Its four-candidate, three-seed development selection differs from Adam’s fixed settings (Section S5). Quality is defined only for valid outputs with qout = qref > 0, against the current moment. Zero-target and unequal-rank geometry remain undefined. A target boundary is flagged tied when its squared-spectrum gap is at most 10−12 max{1, σ12 }. This numerical flag includes near ties; an exactly tied boundary has an objective value but does not identify a unique projector. Shortfall counts positive-reference requests without valid same-rank output. Staleness overlaps these categories. These conventions differ from the direct-rank reconstruction conventions in Section S3. Freshness certifies when a block was constructed, not its geometric accuracy. An old cache can satisfy the current-rank requirement at one request and fail it at another; denominators and stale counts therefore accompany conditional geometry.

5.2

Request schedules and service time

The primary interval A0 includes causal scaling, feature mapping, moments, maintenance, and requested output construction after initialization. Initialization is reported separately and explicitly added for an inclusive total. Stream preparation, offline references, untimed quality checks, hashing, and serialization are excluded. Separate pass B measures inclusive batch means and synchronous query durations; pass C checks quality. They neither decompose A0 additively nor provide extra input seeds, and the synchronous experiment has no external-arrival queue. Table 5 shows every regime–schedule combination. Entries are medians of five seed-paired Adam/exact-request A0 wall-time ratios, not ratios of medians. Exact-request is faster in all 75 classification seed pairs. Adam is faster for every-observation R128 requests and the two densest periodic R256 schedules. R128 period ten favors Adam in only three of five pairs; sparser regression requests favor exact-request at the group median. This grid does not establish a universal crossover frequency.

16

Table 5: All 25 request-schedule comparisons. Values below one favor Adam; each entry is the median of five paired A0 wall-time ratios against exact-request. R128 period ten favors Adam in only 3/5 pairs. The sparse R128 anomaly remains included; ranges, faster counts, and initialization-inclusive comparisons appear in Table S8.

Regime

Every 1

Every 10

Every 100

Every 1000

Bernoulli

C2 C4 C11 R128 R256

1.726 1.806 1.513 0.277 0.085

2.210 2.259 2.628 0.923 0.567

2.242 2.488 2.830 1.841 1.625

2.044 2.452 2.886 2.097 1.622

2.199 2.541 2.899 1.719 1.435

For R128 period 1000, seed 6204 has exact-request A0 wall time 3.225 seconds and process-CPU time 0.625 seconds. Its cause is unresolved; it remains in the ranges and faster counts. Prescribed identical-workload repeats and accumulation-only controls are retained in Table S9; they do not establish grid-wide timing stability.

5.3

Subspace quality, coverage, and staleness

The joint every-observation comparison covers 33 applicable method–regime combinations, a different grid from the 25 schedule contrasts. Table 6 shows both regression workloads; Table S10 and Table S11 contain the complete classification and regression results. Conditional quality averages valid matched-rank requests within each seed, then weights the five seed means equally. Coverage and stale counts pool requests; they are not independent replications. Absolute reference times are directly summarized measured times, not products of ratio medians.

17

Table 6: Every-observation service in R128 and R256. Time is the median paired A0 ratio against exactrequest; headings give its absolute median. Conditional geometry averages within seed before equal seed weighting. Matched/positive-reference and shortfall count coverage; stale overlaps other counts. Both cache periods and PG remain visible beside Adam.

Mean angle (deg)

Matched / positive

Shortfall

Stale

R128: Exact-request A0 wall median = 5,642 ms Exactrequest 1.000 1.0000 0 Adam 0.277 0.9905 0.658 PG 0.280 0.9971 0.447 Cache-10 0.193 0.9979 0.0951 Cache-100 0.107 0.9930 0.299

1.8e-06 28.83 17.11 3.48 9.92

20475/20475 20475/20475 20475/20475 20435/20475 19985/20475

0 0 0 40 490

0 0 0 18390 19785

R256: Exact-request A0 wall median = 18,840 ms Exactrequest 1.000 1.0000 0 Adam 0.085 0.9719 2.31 PG 0.113 0.9201 2.98 Cache-10 0.134 0.9931 0.252 Cache-100 0.047 0.9612 0.815

2.2e-06 87.67 87.94 8.67 23.28

10235/10235 10210/10235 10226/10235 10187/10235 9745/10235

0 25 9 48 490

0 0 0 9175 9645

Service

Paired A0 ratio

Mean ρ

Mean dP

In R128, Cache-10 combines lower cost and higher conditional quality than Adam with 40 shortfalls and 18,390 stale outputs. In R256, Adam’s time ratio is 0.085 but mean dP = 2.31 and mean maximum angle is 87.67◦ . PG has higher equally seed-weighted full-stream mean objective capture than Adam in C2, C4, C11, and R128, but lower capture in R256; R128 endpoint ordering differs from stream-average ordering. The prescribed shortfall-first rule selected γ = 0.1 for R256 despite another candidate’s smaller objective-loss average. Contrast extraction is slower than generic exact-request in C4 and slightly faster in C11, while retaining near-unit objective capture and very small geometric residuals. Cache-10 and Cache-100 cost less than Adam in all three classification regimes, but include stale returns and early shortfalls. These orderings concern the measured implementations and full service boundaries, not complexity alone. Across schedules, 689,028 requests include 688,828 positive references and 685,784 valid matchedrank outputs. The 3,044 shortfalls comprise 2,915 missing valid outputs and 129 insufficient-rank returns; 181,473 stale responses overlap these categories. Undefined geometry is not replaced by zero. Thus low service time may describe a missing, stale, or rank-limited response.

18

R256 - every-observation requests

Objective ratio ρ

Adam and PG

Exact and cached

1.00 0.75 0.50 0.25

80 60 40 20

Projector distance dP

Maximum angle (deg)

0.00

0 4 3 2 1 0 500

1000

1500

2000

500

Observation / request index

1000

1500

2000

Observation / request index

Adam PG

Exact-each

Exact-request

Cache-10

Cache-100

Input seed (line style) 6200

6201

6202

6203

6204

Exact-each and Exact-request coincide; distinct markers identify both.

Figure 3: Complete R256 every-observation trajectories for all five seeds. Columns separate Adam/PG from exact/cached services on identical axes; rows show objective ratio, maximum principal angle, and projector distance. Undefined matched-rank metrics remain gaps. Line styles identify seeds (6202 is dotted) and colours/markers identify services. R128’s distinct full-stream behavior is shown in Fig. S1.

19

R256’s geometric discrepancy persists beyond its early low-objective period (Fig. 3). All 5,120 recorded Adam angles at observations 1025–2048 exceed 80◦ , with valid matched-rank, fresh outputs and no flagged target tie. This window and threshold are post hoc descriptions. The separate saved-state correspondence covers five terminal responses, not all scalar requests (Table S12). Maximum angle concerns the least-aligned direction; projector distances also indicate aggregate error. Small relative boundary gaps make such geometry compatible with high objective capture through Eq. (13), without identifying an optimizer failure or establishing superiority over independent leading-k services.

6

Discussion

6.1

Interpretation of the findings

The same-target intervention separates numerical accuracy from the predictive deficit: improving the terminal solution leaves the regression disadvantage relative to RFF-CIPLS. It narrows an optimization-based explanation without selecting a unique cause among target encoding, feature mapping, coordinate scaling and the separately fitted probe. The intermediate-state and population diagnostics identify what this intervention leaves open. Terminal summaries miss an observed intermediate deviation, while empirical agreement does not establish population recovery. Neither the finite sample-size grid nor the saved reconstruction checkpoints establish a limiting practical discrepancy or convergence rate.

6.2

Scope and limitations

The cost advantage of maintenance is conditional on the requested output. R256 illustrates a fresh, matched-rank service with substantial geometric error; R128’s fast, conditionally accurate Cache-10 instead incurs stale returns and rank shortfalls. Neither low time nor nominal coverage can stand in for current-subspace quality. Relative gaps explain compatibility between energy and geometry, not the cause of an optimizer’s behavior. An important computational alternative remains untested. Avoiding repeated full-spectrum computation need not require the Adam update. Projected gradient already supplies a warm-started small-block method, but does not replace a comparison against a distinct leading-k solver with a specified work subspace and budget or an incremental low-rank update. The present measurements do not establish either alternative’s accuracy or speed. Instrumented batch means are not individualupdate tail latency, and synchronous responses are not deployment queueing measurements. The prediction studies retain further limitations. Multiclass CIPLS uses OVR concatenation and random projection; changing the rank cap changes both operations. Shared input features do not remove target-construction differences, coordinatewise evaluation or asymmetric tuning. The numerical-service comparisons do not resolve those predictive confounders, and PG’s developmentselected step differs from Adam’s fixed setting. The fixed-map theory also excludes the coupled effects of evolving standardizers and the actual adaptive update. Existing stochastic and Riemannian 20

analyses use different assumptions or constructions [35–37]. The six predictive tasks and bounded numerical workloads are not application deployments or concept-drift studies [44]. The resulting diagnosis is conditional: the target, its practical solution and the requested service require separate evaluation. Their agreement cannot be presumed from any single predictive, objective or timing summary, nor do these comparisons select a universally preferable method. Reproducibility. The arXiv source package contains the manuscript source and publication assets required to build this paper. A separate local package contains frozen display inputs and rendering scripts, but is not publicly archived at the time of this version. These materials do not constitute a complete experimental rerun package.

References [1] Eric Bair, Trevor Hastie, Debashis Paul, and Robert Tibshirani. Prediction by supervised principal components. Journal of the American Statistical Association, 101(473):119–137, 2006. doi: 10.1198/ 016214505000000628. [2] Elnaz Barshan, Ali Ghodsi, Zohreh Azimifar, and Mansoor Zolghadri Jahromi. Supervised principal component analysis: Visualization, classification and regression on subspaces and submanifolds. Pattern Recognition, 44(7):1357–1371, 2011. doi: 10.1016/j.patcog.2010.12.015. [3] Theodosios Papazoglou and Guosheng Yin. Covariance supervised principal component analysis. arXiv:2506.19247, 2025. Preprint. [4] Tomojit Ghosh and Michael Kirby. Supervised dimensionality reduction and visualization using centroidencoder. Journal of Machine Learning Research, 23(20):1–34, 2022. URL https://www.jmlr.org/ papers/v23/20-188.html. [5] Arthur Gretton, Olivier Bousquet, Alex Smola, and Bernhard Schölkopf. Measuring statistical dependence with Hilbert–Schmidt norms. In Algorithmic Learning Theory, volume 3734 of Lecture Notes in Computer Science, pages 63–77. Springer, 2005. doi: 10.1007/11564089_7. [6] Amir-Hossein Karimi, Alexander Wong, and Ali Ghodsi. SRP: Efficient class-aware embedding learning for large-scale data via supervised random projections, 2018. URL https://arxiv.org/abs/1811.03166. Preprint. [7] Harold Hotelling. Relations between two sets of variates. Biometrika, 28(3-4):321–377, 1936. doi: 10.1093/biomet/28.3-4.321. [8] Ildiko E. Frank and Jerome H. Friedman. A statistical view of some chemometrics regression tools. Technometrics, 35(2):109–135, 1993. doi: 10.1080/00401706.1993.10485033. [9] Roman Rosipal and Leonard J Trejo. Kernel partial least squares regression in reproducing kernel Hilbert space. Journal of Machine Learning Research, 2(Dec):97–123, 2001. URL https://www.jmlr. org/papers/v2/rosipal01a.html. [10] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems, volume 20, pages 1177–1184, 2007. [11] Bhupinder S. Dayal and John F. MacGregor. Recursive exponentially weighted PLS and its applications to adaptive control and prediction. Journal of Process Control, 7(3):169–179, 1997. doi: 10.1016/ S0959-1524(97)80001-7. [12] S. Joe Qin. Recursive PLS algorithms for adaptive data modeling. Computers & Chemical Engineering, 22(4-5):503–514, 1998. doi: 10.1016/S0098-1354(97)00262-7.

21

[13] Xun Wang, Uwe Kruger, and Barry Lennox. Recursive partial least squares algorithms for monitoring complex industrial processes. Control Engineering Practice, 11(6):613–632, 2003. doi: 10.1016/S0967-0661(02)00096-5. [14] Xue-Qiang Zeng and Guo-Zheng Li. Incremental partial least squares analysis of big streaming data. Pattern Recognition, 47(11):3726–3735, 2014. doi: 10.1016/j.patcog.2014.05.022. [15] Raman Arora, Andrew Cotter, Karen Livescu, and Nathan Srebro. Stochastic optimization for PCA and PLS. In 2012 50th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pages 861–868. IEEE, 2012. doi: 10.1109/Allerton.2012.6483308. [16] Raman Arora, Poorya Mianjy, and Teodor Marinov. Stochastic optimization for multiview representation learning using partial least squares. In Proceedings of the 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1786–1794. PMLR, 2016. URL https://proceedings.mlr.press/v48/aroraa16.html. [17] Zhehui Chen, Lin F. Yang, Chris Junchi Li, and Tuo Zhao. Online partial least square optimization: Dropping convexity for better efficiency and scalability. In Proceedings of the 34th International Conference on Machine Learning, volume 70 of Proceedings of Machine Learning Research, pages 777–786. PMLR, 2017. URL https://proceedings.mlr.press/v70/chen17h.html. [18] Alexander Kaiser, Wolfram Schenck, and Ralf Möller. Coupled singular value decomposition of a cross-covariance matrix. International Journal of Neural Systems, 20(4):293–318, 2010. doi: 10.1142/ S0129065710002437. [19] Xiaowei Feng, Xiangyu Kong, Donghui Xu, and Jianqiang Qin. A fast and effective principal singular subspace tracking algorithm. Neurocomputing, 267:201–209, 2017. doi: 10.1016/j.neucom.2017.06.006. [20] Yao Xie, Ruiyang Song, Hanjun Dai, Qingbin Li, and Le Song. Online supervised subspace tracking. arXiv:1509.00137, 2015. Preprint. [21] Zhanrui Cai, Runze Li, and Liping Zhu. Online sufficient dimension reduction through sliced inverse regression. Journal of Machine Learning Research, 21(10):1–25, 2020. URL https://www.jmlr.org/ papers/v21/18-567.html. [22] Jianjun Xu, Yue Zhao, and Haoyang Cheng. Online kernel sliced inverse regression. Computational Statistics & Data Analysis, 203:108071, 2025. doi: 10.1016/j.csda.2024.108071. [23] Matthew Brand. Fast low-rank modifications of the thin singular value decomposition. Linear Algebra and its Applications, 415(1):20–30, 2006. doi: 10.1016/j.laa.2005.07.021. [24] Erkki Oja. Simplified neuron model as a principal component analyzer. Journal of Mathematical Biology, 15(3):267–273, 1982. doi: 10.1007/BF00275687. [25] Tat-Jun Chin and David Suter. Incremental kernel principal component analysis. IEEE Transactions on Image Processing, 16(6):1662–1674, 2007. doi: 10.1109/TIP.2007.896668. √ [26] Enayat Ullah, Poorya Mianjy, Teodor V. Marinov, and Raman Arora. Streaming kernel PCA with Õ( n) random features. In Advances in Neural Information Processing Systems, volume 31, pages 7311–7321. Curran Associates, Inc., 2018. URL https://proceedings.neurips.cc/paper_files/paper/2018/ hash/7ae11af20803185120e83d3ce4fb4ed7-Abstract.html. [27] Yichuan Deng, Jiangxuan Long, Zhao Song, Zifan Wang, and Han Zhang. Streaming kernel PCA algorithm with small space. In Conference on Parsimony and Learning, volume 280 of Proceedings of Machine Learning Research, pages 1216–1254. PMLR, 2025. URL https://proceedings.mlr.press/ v280/deng25a.html.

22

[28] Artur Jordão, Maiko Lie, Victor Hugo Cunha de Melo, and William Robson Schwartz. Covariancefree partial least squares: An incremental dimensionality reduction method. In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision (WACV), pages 1420–1428, 2021. doi: 10.1109/WACV48630.2021.00146. [29] Hassan Ashtiani and Ali Ghodsi. A dimension-independent generalization bound for kernel supervised principal component analysis. In Proceedings of the 1st International Workshop on Feature Extraction: Modern Questions and Challenges at NIPS 2015, volume 44 of Proceedings of Machine Learning Research, pages 19–29. PMLR, 2015. URL https://proceedings.mlr.press/v44/Ashtiani2015.html. [30] Ky Fan. On a theorem of Weyl concerning eigenvalues of linear transformations I. Proceedings of the National Academy of Sciences of the United States of America, 35(11):652–655, 1949. doi: 10.1073/pnas.35.11.652. [31] Michael L. Overton and Robert S. Womersley. On the sum of the largest eigenvalues of a symmetric matrix. SIAM Journal on Matrix Analysis and Applications, 13(1):41–45, 1992. doi: 10.1137/0613006. [32] P.-A. Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2008. [33] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference on Learning Representations, 2015. URL https://openreview.net/forum?id=8gmWwjFyLj. [34] Gary Bécigneul and Octavian-Eugen Ganea. Riemannian adaptive optimization methods. In International Conference on Learning Representations, 2019. URL https://openreview.net/forum?id=r1eiqi09K7. [35] Silvère Bonnabel. Stochastic gradient descent on Riemannian manifolds. IEEE Transactions on Automatic Control, 58(9):2217–2229, 2013. doi: 10.1109/TAC.2013.2254619. [36] Simon Vary, Pierre Ablin, Bin Gao, and Pierre-Antoine Absil. Optimization without retraction on the random generalized Stiefel manifold. In Proceedings of the 41st International Conference on Machine Learning, volume 235 of Proceedings of Machine Learning Research, pages 49226–49248, 2024. URL https://proceedings.mlr.press/v235/vary24a.html. [37] Hiroyuki Sakai and Hideaki Iiduka. A general framework of Riemannian adaptive optimization methods with a convergence analysis. Transactions on Machine Learning Research, 2025. ISSN 2835-8856. URL https://openreview.net/forum?id=knv4lQFVoE. [38] Iosif Pinelis. Optimum bounds for the distributions of martingales in Banach spaces. The Annals of Probability, 22(4):1679–1706, 1994. doi: 10.1214/aop/1176988477. [39] Chandler Davis and William M. Kahan. The rotation of eigenvectors by a perturbation. III. SIAM Journal on Numerical Analysis, 7(1):1–46, 1970. doi: 10.1137/0707001. [40] 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. [41] Jorge Reyes-Ortiz, Davide Anguita, Alessandro Ghio, Luca Oneto, and Xavier Parra. Human activity recognition using smartphones, 2013. URL https://doi.org/10.24432/C54S4K. [42] Martyna Bator. Dataset for sensorless drive diagnosis, 2013. URL https://doi.org/10.24432/C5VP5F. [43] OpenML. kin8nm. OpenML data ID 189, version 1, 2014. URL https://www.openml.org/d/189. Uploaded 2014-04-23; accessed 2026-07-23. [44] Daniel Bienstock, Minchan Jeong, Apurv Shukla, and Se-Young Yun. Robust streaming PCA. In Advances in Neural Information Processing Systems, volume 35, pages 4231–4243. Curran Associates, Inc., 2022. doi: 10.52202/068431-0306. URL https://proceedings.neurips.cc/paper_files/paper/ 2022/hash/1b11d918b08f781a6c194c6c522edfd6-Abstract-Conference.html.

23

Supplementary material S1

Proofs and implementation details

S1.1

Update rules and numerical conventions

In the predictive configuration, the update in Eq. (9)–Eq. (10) uses fixed α = 0.03, c = 0.5, (β1 , β2 ) = (0.9, 0.999) and ϵ = 10−8 outside the square root. The moment step is 1/t, forgetting is absent, and the optional displacement clip and maximum-basis-change restriction are disabled. The first zero-gradient observation advances the optimizer counter to one. During initial zero moments the initialized span remains unchanged up to QR roundoff; later zero gradients can still produce displacement through nonzero Adam history. The arrays are not transported. The formula box describes the normal finite-data path, not exception handling: the predictive implementation increments the optimizer counter before clipping, whereas the service implementation does so afterward; both use the new counter in bias correction. Signed QR flips columns only for negative diagonal entries of R; zero receives sign +1. Rankdeficient inputs are not rejected, and the condition number of R is diagnostic, not an acceptance threshold. Completion at rank deficiency is not a uniquely defined smooth local retraction. Factorization exceptions propagate without rollback. The predictive core rejects nonfinite raw input before incrementing its sample counter. Nonfinite mapped features instead return an error after advancing that counter but before updating moments or the optimizer, without rollback. A nonfinite post-QR basis triggers Gaussian reinitialization seeded by the current sample counter while retaining advanced moments and optimizer state. These branches do not define accepted-sample skipping. The plain wrapper discards status, but the terminal and rank-aligned studies check every core-update status and reject runs unless all expected updates succeed. Stored clipping counters and norm diagnostics do not supply missing historical branch counts. Each standardizer transforms before its Welford update. An empty state returns the raw observation; with fewer than two prior observations it centers without scale division. Thereafter variance uses the count as divisor, is clamped below at 10−8 , and a coordinate is divided only when its variance exceeds 10−7 . The input RFF convention in Eq. (6) corresponds to exp{−∥x − x′ ∥22 /(2σx2 )}. Regression uses analogous scalar-target features; Dy = Dx in the predictive protocol. The finitefeature bound follows from ∥ϕ∥22 = (2/Dx )

2 j cos (·) ≤ 2. It also holds with changing standardizers,

P

without extending fixed-map probabilistic assumptions. Gradient derivation and convergence assumptions. For Ft (U ) = tr(U ⊤ Ct Ct⊤ U ) and a variation E, symmetry gives dFt (U )[E] = 2 tr{(Ct Ct⊤ U )⊤ E}. Differentiating U ⊤ U = Ik gives U ⊤ Ξ + Ξ⊤ U = 0; subtracting U sym(U ⊤ G) is its orthogonal tangent projection under the Frobenius metric. This proves the direction used before clipping. It does not justify transporting the adaptive

24

arrays, which the implementation does not do, or treating their elementwise-scaled displacement as tangent. The standard stochastic-retraction result requires smoothness, conditionally unbiased uniformly bounded tangent gradients, Robbins–Monro steps, a twice continuously differentiable retraction and compact containment [35]. Under these conditions objective values converge almost surely and gradient norms vanish. The implemented fixed steps, plug-in moments, ambient adaptation, clipping and evolving standardizers are outside that argument; a separate coupled-recursion analysis would be needed. Neither finite-grid behavior nor exact moment identities supply the missing assumptions.

S1.2

Moment identities, bias, and consistency

At t = 1, the recursion gives µϕ,1 = ϕ1 , µψ,1 = ψ1 , R1 = ϕ1 ψ1⊤ and C̄1 = 0. If a running average is correct at t − 1, then (1 − t−1 )

1 X 1X ai + t−1 at = ai . t − 1 i<t t i≤t

Apply this to both features and their outer product. Expanding the centered products subtracts the mean outer product twice and adds it once, proving Eq. (8) pathwise. For the probabilistic claims, condition on fixed maps and preprocessing and assume i.i.d. pairs with finite second moments. Write M0 = E[ϕψ ⊤ ]. Diagonal and off-diagonal terms give −2 ⊤ E[µϕ,t µ⊤ ψ,t ] = t {tM0 + t(t − 1)µϕ µψ }.

Independence is across observations, not within pairs. Subtraction from ERt = M0 gives EC̄t = (t − 1)C∗ /t, so Ct = tC̄t /(t − 1) is unbiased for t ≥ 2 under these conditions only. Cauchy–Schwarz gives E|ϕa ψb | ≤ (Eϕ2a Eψb2 )1/2 < ∞. The scalar strong law therefore applies to every raw-moment entry and feature mean. Their limits hold jointly in finite dimensions, and continuity of the mean product proves C̄t → C∗ almost surely. Second moments alone supply no concentration rate.

S1.3

Concentration under bounded features

Fix t ≥ 1 and δ ∈ (0, 1) under the additional bounded-feature conditions. If Bx By = 0, one feature is zero almost surely and the result is immediate. Otherwise the centered raw outer products and feature increments obey ∥ϕi ψi⊤ − M0 ∥F ≤ 2Bx By , ∥ϕi − µϕ ∥2 ≤ 2Bx ,

∥ψi − µψ ∥2 ≤ 2By .

These follow from the outer-product norm identity, Jensen’s inequality and the whole-vector bounds. Matrices with the Frobenius inner product form a Hilbert space, as do the feature spaces. For any of these centered sequences Zi with bound L, use the observation filtration and the stopped martingale Sj =

P

i≤j∧t Zi .

Independence across pairs makes its increments conditionally mean

25

zero, with squared essential-supremum bounds summing to at most tL2 . Hilbert-space smoothness constant one in Pinelis [38, Theorem 3.5] yields Pr(∥St /t∥ ≥ u) ≤ Pr(max ∥Sj ∥ ≥ tu) ≤ 2 exp{−tu2 /(2L2 )}. j≤t

Set at =

p

2 log(6/δ)/t. Taking u = Lat makes each failure probability at most δ/3. The union

bound, requiring no independence between the three events, gives probability at least 1 − δ for ∥Rt − M0 ∥F ≤ 2Bx By at , ∥µϕ,t − µϕ ∥2 ≤ 2Bx at , ∥µψ,t − µψ ∥2 ≤ 2By at . The exact centering decomposition is ⊤ ⊤ µϕ,t µ⊤ ψ,t − µϕ µψ = (µϕ,t − µϕ )µψ,t

+ µϕ (µψ,t − µψ )⊤ . The sample target mean in the first term retains the product of mean errors. Since ∥µψ,t ∥2 ≤ By pathwise and ∥µϕ ∥2 ≤ Bx , on the same event C̄t − C∗

op

≤ ∥C̄t − C∗ ∥F ≤ ∥Rt − M0 ∥F + By ∥µϕ,t − µϕ ∥2 + Bx ∥µψ,t − µψ ∥2 q √ ≤ 6Bx By at = 6 2Bx By log(6/δ)/t.

(S1)

For t ≥ 2, expansion of Ct − C∗ = t(C̄t − C∗ )/(t − 1) + C∗ /(t − 1) also gives ∥Ct − C∗ ∥op ≤

∥C∗ ∥op t bt (δ) + . t−1 t−1

Stopping at the chosen t proves a fixed-time statement, not simultaneous control over requests or adaptive stopping times.

S1.4

Exact-subspace perturbation and objective-gap bounds

b ⊤ U∗ . For two orthonormal k-column bases, principal-angle cosines are the singular values of U t

Idempotence and trace cyclicity give b ⊤ U∗ ∥2 = 2 ∥Pbt − P∗ ∥2F = 2k − 2∥U t F

k X j=1

26

sin2 θj .

(S2)

Under the positive population separator ∆k , Yu et al. [40, Theorem 2] bounds the projector distance by

√ √ √ 2 2 2 2k min{ k ∥Mt − M∗ ∥op , ∥Mt − M∗ ∥F } ≤ ∥Mt − M∗ ∥op . ∆k ∆k

No positive empirical gap is required for this population-gap bound. With Et = C̄t − C∗ , expand Mt − M∗ = Et C∗⊤ + C∗ Et⊤ + Et Et⊤ . Triangle inequality and submultiplicativity prove the second deterministic bound in Eq. (16). Under all the probabilistic and gap conditions, substitution on the concentration event yields √

2 2k ∥Pbt − P∗ ∥F ≤ bt (δ){2 ∥C∗ ∥op + bt (δ)}. ∆k

This controls the exact empirical basis, not the implemented iterate. For Eq. (13), express P in an eigenbasis of M and set s =

j≤q (1 − Pjj ) =

P

P

j>q Pjj .

The

diagonal entries lie in [0, 1] and sum to q. Its objective loss is therefore at least (λq − λq+1 )s, while (q)

∥P − P∗ ∥2F = 2s. Rearrangement proves the displayed bound. For binary rank one, projection onto the sole positive direction gives ρ1 = cos2 θ and the stated distance identity. For Eq. (12), every contained orthonormal block is U Q, Q = U ⊤ X ∈ St(k, q); substitution gives tr(Q⊤ BU BU⊤ Q) and Ky Fan yields the leading squared-singular-value sum, including zero or tied cutoffs. For completeness, the classification column identity follows by expanding E[(ϕ − µ)(1{y = c} − pc )] = pc (µc − µ). Multiplying the moment by its transpose gives the squared-probability weights. Centered one-hot vectors have zero coordinate sum; multiplying either the population expectation or realized sample-centered sum by 1C proves its null vector and the rank upper bound. This argument requires no class-conditional mean for a zero-probability class.

S1.5

Storage and arithmetic complexity

The core arrays are R of shape Dx × Dy , feature means of lengths Dx , Dy , and U, m, v of shape Dx × k, plus scalar counters. Dense maps retain input weights/phases of shapes Dx × dx and Dx and, for regression, target weights/phases of shapes Dy × dy and Dy . Scalers retain means, variances and counts. The predictive implementation additionally allocates feature second-moment buffers of total length Dx + Dy with whitening disabled, internal input-scaler buffers of total length 2dx with internal centering disabled, and a Dx × k previous-basis diagnostic copy. These unused or diagnostic arrays remain part of that implementation’s fixed-dimensional storage; lean numerical harnesses need not allocate them. Scalar/class metadata, monitor lists and random-generator objects are distinct from these named arrays. The bound excludes raw-data/generator buffers, wrapper objects, retained responses, transient factorization workspaces and process RSS. Dense maps, moment updates, associative gradient

27

products and thin QR supply the main arithmetic terms; conditioning the small QR factor adds O(k 3 ). None of these bounds establishes observed scaling or latency.

S2

Prediction protocols and paired comparisons

S2.1

Data, feature maps, and parameter selection

Each synthetic final seed s = 100, . . . , 109 initializes a fresh RandomState(s) and generates 5,000 observations: the first 4,000 train and the last 1,000 test. Seed 10000 + s permutes training observations only. XOR-Gauss samples equally from centers (−1, −1), (−1, 1), (1, −1), (1, 1) with labels 0, 1, 1, 0, adds independent Gaussian noise of SD 0.3, and appends 18 standard-Gaussian nuisance coordinates. Spirals-3 draws independent u ∼ Unif[0, 1], c ∼ Unif{0, 1, 2}, and ϵ ∼ N (0, 0.22 ); its two signal coordinates are (r cos θ, r sin θ),

r = u + ϵ,

θ = 4πu + 2πc/3.

The label is c, negative radii are retained, and 18 independent Gaussian nuisance coordinates are appended. Friedman-1 uses ten independent uniform [0, 1] coordinates and independent ϵ ∼ N (0, 1): y = 10 sin(πx1 x2 ) + 20(x3 − 0.5)2 + 10x4 + 5x5 + ϵ. HAR preserves the official subject-disjoint test set; development uses GroupShuffleSplit, fraction 0.2 and seed 20260715, yielding 5,588 training and 1,764 validation observations. Sensorless uses the fixed stratified 60/20/20 development construction, then combines development training/validation for the final training partition in Table 1. Kin8nm preserves its seed-0 80/20 split and splits only training again at fraction 0.2 with seed 20260715. Real-data seed s changes stream order, input RFF, initialization, and probe random state through 10000 + s, 20000 + s, 30000 + s, and 40000 + s, respectively, rather than held-out membership. Regression target-RFF seed is input-RFF seed plus one.

S2.2

Comparators and prediction probes

For C-class OVR–RP–CIPLS, each scalar-target reducer receives +1/−1 and produces k coordinates. Its compressed output is r(x) = Π⊤ col(z1 (x), . . . , zC (x)),

Π ∈ RCk×k ,

Πij ∼ N (0, 1/(Ck)).

The projection is sampled once with RandomState(20000+s). RFF class reducers use identical input maps but store separate copies. Changing k changes both reducer and compression dimensions. Binary and regression CIPLS use scalar targets without this adapter. All CIPLS constructions retain tolerance 10−12 and no norm clip; CCIPCA’s inherited amnesic parameter is two. Table S2 distinguishes these constructions and their unequal tuning opportunities. 28

Table S1: Predictive settings and development-selected bandwidths. The explicit OKSPCA rate/clip grid and selection rule are distinguished from inherited comparator settings. Bandwidth calibration uses a development-training subset of at most 512 observations; final-test data do not select parameters. Feature budgets appear in Table 1.

Item

Setting / selection

Optimizer

Adam-style QR; rate 0.03; clip 0.5; β1 = 0.9, β2 = 0.999; ϵ = 10−8 outside square root; no decay. Rates {0.003, 0.01, 0.03}; clips {0.5, 1.0}; XOR-Gauss, HAR, Kin8nm; seeds 0–4; 90 runs. Minimum average rank; ties by smaller rate, then smaller clip. Both clips tied at rate 0.03. √ Median positive pairwise distance divided by 2 after causal standardization on a fixed development-training subset of at most 512 observations; then frozen. Regression uses the same bandwidth rule on the scalar standardized response; target RFF seed is input RFF seed plus one. Forgetting factor 1 and empirical step 1/t; fixed protocol choices.

Development Selection Bandwidth rule

Target mapping Moment Task XOR-Gauss Spirals-3 Friedman-1 HAR Sensorless Kin8nm

σx

σy (regression)

4.49429 4.49160 3.20342 23.0553 6.40547 2.81574

— — 0.709610 — — 0.710298

Table S2: Method, target, adapter, and selection identities. Shared input features do not imply shared response construction or symmetric tuning. Raw-X and RFF-X are unreduced controls with different output dimensions from the reduced pipelines.

Pipeline

Input

Target / construction

Selection opportunity

OKSPCA RFF-CIPLS / RFF–OVR– RP–CIPLS CIPLS / OVR–RP–CIPLS CCIPCA RFF-CCIPCA Raw-X RFF-X

Shared RFF Shared RFF

One-hot / target RFF Scalar / OVR with fixed Ck → k RP

Rate and clip selected Inherited; no method grid

Raw linear

Scalar / OVR with fixed Ck → k RP No target; covariance No target; covariance No reduction; dx outputs No reduction; Dx outputs

Inherited; no method grid

Raw linear Shared RFF Raw identity Shared RFF

29

Inherited; amnesic 2 Inherited; no method grid Untuned control Untuned control

Terminal probes use scikit-learn 1.7.1. Logistic regression uses L2 regularization, lbfgs, tolerance 10−4 , at most 2,000 iterations, no class weights, and an unpenalized intercept.

With n training

embeddings its objective is mean negative log-likelihood plus ∥W ∥2F /(2n); Ridge uses the dense direct-solver path and minimizes squared-error sum plus ∥w∥22 , also with an unpenalized intercept. Each representation receives its own training-embedding standardizer. The replay definition is Eq. (18).

S2.3

Paired effects and unreduced controls

Table S3 separates three contrasts. The original regression primary contrast is OKSPCA minus RFF-CIPLS; classification rank-cap contrasts use the appropriately named scalar or OVR–RP comparator; the same-target regression interval concerns extracted minus exact. All use 10,000 resamples of ten paired seed differences with replacement, averaging within resample and reporting 2.5th/97.5th percentiles. For the rank-cap and extracted-block extensions, each estimand starts a fresh RandomState(50000) and uses linear percentile interpolation. The original regression inputs and reported intervals are retained, but the exact historical reporting-code and RNG-state binding has not been recovered; the later extensions’ seed cannot be assigned to that analysis. Its stored intervals are reproduced as values, not claimed as a bitwise-reproduced bootstrap. Later extensions are conditional descriptions of reused seeds, not independent confirmation. No interval here tests predictive equivalence; exact-minus-full changes in Table 2 remain separate descriptive quantities. The four primary classification caps are one, two, five, and ten for XOR-Gauss, Spirals-3, HAR, and Sensorless. Every reduced pipeline is fitted anew at that cap with its selected or inherited settings unchanged. The classification intervals apply to paired method differences, not to the individual score columns. Regression benchmark dispersion is the sample standard deviation across ten seeds; it is not a confidence interval for a population mean. Reusing a fixed test partition in real-data repetitions also differs from generating a new synthetic sample. Accordingly, the intervals have task-specific conditional interpretations and should not be treated as repeated independent test-set draws.

30

Table S3: Retained paired predictive contrasts and their original 95% percentile-bootstrap intervals. Classification contrasts are rank-cap OKSPCA minus the named comparator; primary regression contrasts are OKSPCA minus RFF-CIPLS; same-target regression contrasts are extracted minus exact. The last are not exact-minus-full intervals.

OKSPCA − RFF-CIPLS-based (95% CI)

Task

Metric

XOR-Gauss Spirals-3 HAR Sensorless Friedman-1 Kin8nm

Macro-F1 Macro-F1 Macro-F1 Macro-F1 R2 R2

-0.0002 [-0.0010, +0.0008] +0.0170 [+0.0029, +0.0343] +0.0504 [+0.0334, +0.0667] +0.1513 [+0.1186, +0.1860] -0.0564 [-0.0706, -0.0416] -0.0416 [-0.0489, -0.0352]

Task

Metric

Extracted − exact (95% CI)

Friedman-1 Kin8nm

R2 R2

-0.00014 [-0.00153, +0.00099] -0.00003 [-0.00021, +0.00013]

Table S4: No-reduction controls and their original-budget OKSPCA reference. Entries are ten-final-seed means and sample standard deviations. Raw-X has dx coordinates and RFF-X has Dx ; OKSPCA uses the displayed original kold . Classification entries are not rank-cap control experiments. Different output dimensions preclude including these controls in reduced-method ranks.

Task

Metric

kold

Old OKSPCA

Raw-X

RFF-X

XOR-Gauss Spirals-3 Friedman-1 HAR Sensorless Kin8nm

Macro-F1 Macro-F1 R2 Macro-F1 Macro-F1 R2

4 5 4 20 12 4

0.605 ± 0.035 0.369 ± 0.013 0.663 ± 0.042 0.883 ± 0.009 0.773 ± 0.025 0.413 ± 0.025

0.504 ± 0.016 0.378 ± 0.009 0.723 ± 0.016 0.954 ± 0.001 0.913 ± 0.000 0.404 ± 0.000

0.660 ± 0.032 0.360 ± 0.010 0.853 ± 0.015 0.939 ± 0.005 0.945 ± 0.003 0.657 ± 0.033

31

S3

Same-target diagnostic protocols

S3.1

Primary-rank classification reconstruction

The 40 fits in Table 3 start at the primary classification ranks with the original inputs and settings, reproduce the terminal scores, and retain five checkpoints each. They reconstruct states rather than recover historical states. Numerical rank counts squared singular values strictly above 10−10 σ12 in the full saved moment; comparison rank is the smaller of this rank and the trained cap. At equal rank, geometry uses the maintained span; at lower positive rank it uses within-span extraction without changing the state. The original-rank terminal probe is retained. A smaller exact comparison dimension remains a mixed-dimension comparison, without adding a third extracted-block probe. At rank zero, basis and normalized geometric diagnostics are undefined and the exact objective is zero. For positive rank with 0 < J ∗ ≤ 10−15 , raw objectives remain recorded but normalized diagnostics are undefined. An exact tie has zero boundary gap; the numerical-tie flag uses gap at most 10−10 λ1 . These are reconstruction conventions, distinct from the service rule in Section 5. Table S5 retains all five Sensorless seed-104 checkpoints, including the adverse midpoint. They do not establish the deviation’s duration or cause. For two orthonormal equal-rank bases, principal angles are determined by the singular values of their cross-product; the largest acute angle reports the least-aligned direction. Projector distance is ∥U U ⊤ − V V ⊤ ∥F , without division by rank. Objective ratios use the same empirical moment and comparison dimension. These statistics cannot be interchanged with prediction scores, and a tied boundary does not select a unique projector even when the leading objective value remains well defined. Table S5: The five retained Sensorless seed-104 checkpoints of the directly trained primary-rank reconstruction. Rank, objective capture, projector distance, and maximum angle retain the same-state definitions. Geometric entries are rounded to six decimals; displayed ρ = 1.000000 does not assert exact equality. The observation-23,404 deviation is a measured checkpoint, not a uniform pathwise statement.

Observation 4681 11702 23404 35106 46807

ktrain

rτ

q

ρ

θmax (deg)

dP

10 10 10 10 10

10 10 10 10 10

10 10 10 10 10

0.999978 0.999922 0.970687 0.999905 1.000000

2.622743 0.689740 78.099077 0.655705 0.163814

0.087174 0.029430 1.898569 0.023416 0.006250

32

S3.2

Regression representation comparisons

The regression endpoints use final seeds 100–109 from the larger-span terminal study. Full is the maintained k = 4 representation; extracted is its maximizing within-span block as in Eq. (12); exact is the leading left singular block of the same terminal C̄T . Positive rescaling CT = T C̄T /(T − 1) does not change its singular subspaces. The full-moment threshold is 10−10 σ12 , with qT = min(k, rτ ); both regressions have qT = k = 4. Each representation has its own standardizer and penalized probe. Extraction can rotate coordinates even at equal rank. Mean full-minus-extracted R2 is −3.897 × 10−7 for Kin8nm and −5.315 × 10−7 for Friedman-1; maximum absolute seed differences are 1.358 × 10−6 and 1.956 × 10−6 . Their small size does not establish rotation invariance. Table S6 reports extracted-to-exact geometry. The extracted-minus-exact intervals in Table S3 and descriptive exact-minus-full changes in Table 2 compare different paired quantities. Table S6: Regression terminal extracted-to-exact objective and geometry, averaged over ten final seeds at qT = 4. Projector distance is unnormalized Frobenius distance; angle is the maximum principal angle in degrees. These are the larger-span study’s regression endpoints, not direct classification-rank reconstructions.

Task

qT

Objective ratio ρ

dP

θmax

Friedman-1 Kin8nm

4 4

0.999969 0.999980

0.1426 0.0300

5.347◦ 1.016◦

S4

Bounded population-recovery protocol

This experiment uses identity features ϕ(x) = x and ψ(y) = y, with no RFF or raw-data standardizer. Its independently specified tracker follows the moment and maintenance stages B–C of Algorithm 1: empirical step 1/t, Bessel-corrected gradient moment for t > 1 (zero at t = 1), ambient Adam with fixed α = 0.03, (β1 , β2 ) = (0.9, 0.999) and ϵ = 10−8 outside the square root. The tangent gradient is clipped globally at Frobenius norm 0.5, using the over-threshold denominator stabilizer 10−12 . There is no learning-rate decay, adaptive scheduling or moment transport. Adam arrays and its counter start at zero; each sample advances the counter and applies signed reduced QR. Initialization is specified below. The RFF/scaler stages A and D of the predictive algorithm are absent here. Draw independent standard Gaussian gx ∈ R32 and ge ∈ R8 , and define x=

√

32 gx /∥gx ∥2 ,

e=

√

8 ge /∥ge ∥2 ,

y = C∗⊤ x + Be,

where B = (I − C∗⊤ C∗ )1/2 is the symmetric positive-semidefinite square root. The four positive squared singular values of C∗ are (0.49, 0.36, 0.25, 0.15), all others zero. Thus ∥C∗ ∥op = 0.7 < 1, the root exists, and (Dx , Dy , r, k) = (32, 8, 4, 3); the fourth direction defines the selected boundary gap.

33

Spherical symmetry and independence give centered coordinates, Cov(x) = I32 , Cov(e) = I8 , and Cov(x, y) = C∗ , Cov(y) = C∗⊤ C∗ + BB ⊤ = I8 . √ √ √ √ Moreover, ∥x∥2 = 32, ∥e∥2 = 8, and ∥y∥2 ≤ ∥C∗ ∥op 32 + ∥B∥op 8 almost surely. These fixed-feature, bounded coordinates satisfy the stated population construction without time-varying standardization. Signed reduced QR of counter-generated Gaussian matrices with orientation seeds 81001 and 81002 fixes the left/right singular vectors across all runs. Seeds 0–19 vary the spherical streams and separately namespaced Gaussian three-column tracker initialization, not population orientations. Sample sizes 256, 512, . . . , 32,768 are nested within seed. The objective curves in Fig. 2 are ρ3 (Uonline ; C̄n ), ρ3 (Uonline ; C∗ ), and ρ3 (Uexact ; C∗ ); the empirical-exact ratio is one by construction. Projector pairs are defined in Eq. (19). Pointwise seed medians and interquartile ranges describe this finite grid, not independent checkpoints or an estimated convergence rate.

S5

Numerical-service protocols and comparisons

S5.1

Streams, query schedules, and parameter selection

For the workloads in Table 4, draw xi ∼ N (0, Idx ). Classification uses C ∈ {2, 4, 11} and selects the label maximizing 1.25{xi1 cos(2πc/C) + xi2 sin(2πc/C)} + gic ,

c = 0, . . . , C − 1,

with independent standard Gumbel gic and the lowest class index breaking exact ties. Regression uses independent ϵi ∼ N (0, 1) and yi = sin(xi1 ) + 0.5(x2i2 − 1) + 0.75xi3 xi4 + 0.25xi5 + 0.25ϵi . Remaining raw coordinates are Gaussian nuisances. All methods within regime and seed share √ stream, maps, order, and the signed-QR Gaussian initial basis. Input RFF bandwidth is dx ; regression target bandwidth is one in causally standardized units. Each service computes the frontend and moment itself. These computational generators have no predictive endpoint. The Adam service uses fixed base rate 0.03, β1 = 0.9, β2 = 0.999, global Frobenius gradient clip 0.5, and ϵ = 10−8 outside the square root. The clipping denominator uses 10−12 . These settings belong to the service implementation; its final runs do not retune them against PG’s development-selected step. p

For the structured multiclass service, Helmert column j = 1, . . . , C − 1 equals 1/ j(j + 1) in rows 1, . . . , j, −j/ j(j + 1) in row j + 1, and zero elsewhere. Hence HC HC⊤ = I − 11⊤ /C. Reduced p

QR of C̄t HC followed by SVD of the triangular factor determines the singular values and input directions; rank is not inferred from QR diagonal entries. Plain QR suffices for an exact full contrast

34

span only if Dx ≥ C − 1, rank(C̄t ) = C − 1, and the requested rank is C − 1. Deficient rank, smaller requested blocks, or canonical singular coordinates require the stated spectral work. Floating-point residual checks against the original moment remain outside timing. Final seeds are 6200–6204; PG development seeds are 6100–6102. Periodic requests use positive multiples of 1, 10, 100, or 1000 and termination, with duplicates removed. Bernoulli requests use probability 0.01 with the fixed seed-index hash construction and include termination. Cache updates retain periods ten or one hundred, with no extra terminal refresh. Same-timestamp extraction can reuse factors; cache age is measured against the request index. Returned bases carry their actual dimension and update timestamp. Approximate services use the spectrum of their inspected small matrix to determine output rank; caches retain the rank at their last scheduled factorization. A no-target response describes the inspected matrix, whereas a zero full-reference rank describes the offline target. Thus these two events need not coincide for approximate or stale services, and request status cannot be reconstructed from a conditional objective mean alone. For PG, each γ ∈ {0.1, 0.3, 1, 3} is tested on the three development seeds. Let Ts contain positive-reference observations and define

ℓs,t (γ) =

 1 − ρs,t (γ),

valid same-rank output,

1,

missing or unequal-rank output,

3 1X L(γ) = 3 s=1

P

t∈Ts ℓs,t (γ)

|Ts |

.

(S3)

Selection first minimizes total shortfalls, then chooses the smallest γ within 10−12 of the minimum loss among those candidates. Zero-reference observations are excluded from loss; numerical failures stop the trajectory. The unit penalty is a development criterion, not an imputation for final undefined geometry. Table S7 retains all 20 candidates, including R256’s selected γ = 0.1. Table S7: All PG development candidates. Cells give total shortfalls and equally seed-weighted prescribed loss; the selected candidate follows shortfall, loss, and step-size priority. The development grid does not constitute a matched tuning budget for fixed-setting Adam.

Regime C2 C4 C11 R128 R256

γ = 0.1

γ = 0.3

γ=1

γ=3

0 / 0.00955 0 / 0.03017 0 / 0.09772 0 / 0.03005 10 / 0.07625

0 / 0.00358 0 / 0.01075 0 / 0.03602 0 / 0.01138 20 / 0.03354

0 / 0.00076 0 / 0.00239 0 / 0.00981 0 / 0.00304 18 / 0.01298

0 / 0.16788 0 / 0.00081 0 / 0.00292 0 / 0.09884 16 / 0.07763

35

S5.2

Timing protocol and repetitions

The recorded execution environment was Windows on an AMD Ryzen 5 7500F six-core processor, using Python 3.13.6, NumPy 2.3.2 and SciPy 1.16.1. The loaded NumPy and SciPy OpenBLAS libraries were versions 0.3.30 and 0.3.28, respectively, each reporting one thread; the process affinity mask was one. Scientific arrays used float64. Generic exact and small-matrix factorizations called numpy.linalg.svd with full_matrices=False, then selected the required columns. No explicit LAPACK-driver option was supplied. The retained environment record does not identify the Windows edition/build or record a per-call LAPACK-driver trace; these are not inferred from the present machine. A0 is measured in a fresh process after initialization. It includes causal preprocessing, features, moments, maintenance, and required output extraction, signed QR, copying, and arithmetic/shape/finiteness checks. Exact-each pays factorization during maintenance. Libraries use one thread and the process one logical processor. Stream preparation, offline-reference construction, quality comparisons, hashing, and serialization are excluded. Table S8 retains the 125 paired base timings through their 25 group summaries, including initialization-inclusive results and the anomalous R128 observation. Pass B measures inclusive 64-observation batch means and synchronous query durations nested within batches. Its quantiles retain all completed calls, including unavailable/no-target responses, without subtracting overhead or dropping resolution-limited timings. These are not individualupdate tail latencies, queueing delays, or additive components of A0. Untimed pass C measures quality and state correspondence. Retaining responses may grow harness memory even though named numerical arrays have fixed dimensions.

S5.3

Quality and request coverage

Tables S10 and S11 cover all 33 applicable every-observation combinations: six services in every regime plus binary C2 and contrast C4/C11. Time ratios first pair A0 by seed against exact-request, then take the median. Absolute medians come directly from measured durations. Quality first averages valid positive-reference, matched-rank requests within seed, then averages five seed means equally. Counts retain request denominators; stale overlaps validity and rank categories. This grid does not replace the distinct schedule comparison. Structured services have near-unit objective ratios with measured floating-point residuals; exact-each and exact-request differ in when they pay for the same current target. Across all schedules, valid-basis, uninitialized-cache, and no-identified-target responses number 685,913, 2,975, and 140, respectively. Their sum is 689,028; these status counts differ from the positive-reference coverage partition in Section 5. Current-rank increases can make an initialized cache rank insufficient. Neither missing nor unequal-rank geometry is filled with zero.

36

Table S8: All 25 Adam/exact-request schedule groups. Paired A0 medians, ranges, and Adam-faster counts retain all five base seed pairs per group. Initialization-inclusive ratios are separate from primary service time. No anomalous timing is removed.

Regime

Spacing

C2 C2 C2 C2 C2 C4 C4 C4 C4 C4 C11 C11 C11 C11 C11 R128 R128 R128 R128 R128 R256 R256 R256 R256 R256

1 10 100 1000 Bernoulli 0.01 1 10 100 1000 Bernoulli 0.01 1 10 100 1000 Bernoulli 0.01 1 10 100 1000 Bernoulli 0.01 1 10 100 1000 Bernoulli 0.01

Median

Min

Max

Faster seeds / 5

Init.-included median ratio

1.726 2.210 2.242 2.044 2.199 1.806 2.259 2.488 2.452 2.541 1.513 2.628 2.830 2.886 2.899 0.277 0.923 1.841 2.097 1.719 0.085 0.567 1.625 1.622 1.435

1.617 2.013 2.090 1.846 1.804 1.716 2.232 2.337 2.117 2.246 1.251 2.317 2.615 2.750 2.654 0.261 0.892 1.744 0.287 1.637 0.081 0.497 1.332 1.455 1.130

1.841 2.427 2.360 2.424 2.289 1.911 2.759 2.584 2.648 2.758 2.013 2.887 2.873 3.239 3.118 0.286 1.144 1.964 2.613 1.775 0.101 0.932 2.283 2.021 2.321

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 3 0 1 0 5 5 0 0 0

1.726 2.209 2.241 2.044 2.198 1.806 2.258 2.486 2.452 2.540 1.512 2.627 2.829 2.884 2.890 0.277 0.923 1.841 2.096 1.719 0.085 0.567 1.624 1.622 1.434

37

R128 - every-observation requests Exact and cached

1.0 0.8 0.6 0.4 0.2

80 60 40 20

Projector distance dP

Maximum angle (deg)

Objective ratio ρ

Adam and PG

0

2

1

0 1000

2000

3000

4000

1000

Observation / request index

2000

3000

4000

Observation / request index

Adam PG

Exact-each

Exact-request

Cache-10

Cache-100

Input seed (line style) 6200

6201

6202

6203

6204

Exact-each and Exact-request coincide; distinct markers identify both.

Figure S1: Complete every-observation R128 trajectories for all five seeds, including early behavior and undefined segments. Objective, maximum-angle, and projector-distance rows use comparable axes across Adam/PG and exact/cached columns. Conditional cache quality does not remove early shortfalls or stale returns. Line styles identify seeds (6202 is dotted); colours and sparse markers identify services.

38

Table S9: All prescribed C4 seed-6200, period-ten repeats: three Adam, three exact-request, and six accumulation-only measurements. Wall and CPU times are in milliseconds, with initialization-inclusive values separate. Accumulation A/B identify repeated launches, not different algorithms. These 12 measurements add no input seeds and establish no grid-wide timing stability.

Service / control Adam Adam Adam Exact-request Exact-request Exact-request Accum.-A Accum.-A Accum.-A Accum.-B Accum.-B Accum.-B

S5.4

Repeat

Service wall (ms)

Service CPU (ms)

Total wall (ms)

Total CPU (ms)

0 1 2 0 1 2 0 1 2 0 1 2

683.60 692.76 661.37 302.66 295.91 299.18 247.13 260.19 266.83 249.48 233.59 258.21

609.38 609.38 656.25 296.88 265.62 296.88 203.12 250.00 218.75 218.75 234.38 218.75

683.79 692.97 661.61 302.85 296.09 299.54 247.31 260.37 267.03 249.70 233.80 258.42

609.38 609.38 656.25 296.88 265.62 296.88 203.12 250.00 218.75 218.75 234.38 218.75

R256 delivered-basis correspondence

Table S12 compares the five terminal delivered query-cache bases with their saved exact references, not with raw maintained bases. Their maximum angles range from 84.109◦ to 89.407◦ . The separate 5,120-request, latter-half description uses recorded scalar diagnostics and a post hoc threshold; it is not a recomputation from all corresponding saved states. For the five endpoints, d2P /2 =

P

j sin

2θ

j

is approximately 1.68, 1.94, 3.00, 2.95, 2.21 in seed order. Every sum exceeds one, excluding a single misaligned direction with all others exactly aligned, but not implying that all directions are orthogonal.

39

Table S10: Complete classification every-observation comparison: all 21 applicable method–regime combinations. Cost, conditional geometry, matched/positive-reference counts, shortfalls, and overlapping stale counts are read jointly. Conditional means do not measure whole-grid success rates.

Mean dP

Mean angle (deg)

Matched / positive

Shortfall

Stale

C2: Exact-request A0 wall median = 492.9 ms Exactrequest 1.000 1.0000 0 Adam 1.726 0.9946 0.0383 PG 1.798 0.9994 0.0013 Cache-10 0.714 0.9977 0.0271 Cache-100 0.572 0.9920 0.0731 Exact-each 0.984 1.0000 0 Binary 0.878 1.0000 1.6e-14

6.6e-07 1.59 0.07 1.11 2.99 6.6e-07 4.5e-07

20473/20473 20473/20473 20473/20473 20435/20473 19985/20473 20473/20473 20473/20473

0 0 0 38 488 0 0

0 0 0 18390 19785 0 0

C4: Exact-request A0 wall median = 594.2 ms Exactrequest 1.000 1.0000 0 Adam 1.806 0.9898 0.199 PG 1.752 0.9990 0.0472 Cache-10 0.624 0.9962 0.102 Cache-100 0.569 0.9797 0.308 Exact-each 1.019 1.0000 0 Contrast 1.129 1.0000 1.5e-15

1.4e-06 7.33 1.91 3.56 10.36 1.4e-06 1.3e-06

20474/20474 20474/20474 20474/20474 20427/20474 19985/20474 20474/20474 20474/20474

0 0 0 47 489 0 0

0 0 0 18390 19785 0 0

C11: Exact-request A0 wall median = 880.8 ms Exactrequest 1.000 1.0000 0 Adam 1.513 0.9812 0.653 PG 1.458 0.9971 0.243 Cache-10 0.516 0.9966 0.229 Cache-100 0.423 0.9756 0.727 Exact-each 1.022 1.0000 0 Contrast 0.948 1.0000 3.8e-15

2.2e-06 16.49 6.89 6.26 15.04 2.2e-06 1.7e-06

20473/20473 20473/20473 20473/20473 20363/20473 19985/20473 20473/20473 20473/20473

0 0 0 110 488 0 0

0 0 0 18390 19785 0 0

Service

Paired A0 ratio

Mean ρ

40

Table S11: Complete regression every-observation comparison: all 12 applicable combinations, including exact-each. Statistics have the same seed weighting and request conditioning as Table S10. R128 and R256 are distinct workloads, not replicates of a width intervention.

Mean angle (deg)

Matched / positive

Shortfall

Stale

R128: Exact-request A0 wall median = 5,642 ms Exactrequest 1.000 1.0000 0 Adam 0.277 0.9905 0.658 PG 0.280 0.9971 0.447 Cache-10 0.193 0.9979 0.0951 Cache-100 0.107 0.9930 0.299 Exact-each 0.979 1.0000 0

1.8e-06 28.83 17.11 3.48 9.92 1.8e-06

20475/20475 20475/20475 20475/20475 20435/20475 19985/20475 20475/20475

0 0 0 40 490 0

0 0 0 18390 19785 0

R256: Exact-request A0 wall median = 18,840 ms Exactrequest 1.000 1.0000 0 Adam 0.085 0.9719 2.31 PG 0.113 0.9201 2.98 Cache-10 0.134 0.9931 0.252 Cache-100 0.047 0.9612 0.815 Exact-each 1.016 1.0000 0

2.2e-06 87.67 87.94 8.67 23.28 2.2e-06

10235/10235 10210/10235 10226/10235 10187/10235 9745/10235 10235/10235

0 25 9 48 490 0

0 0 0 9175 9645 0

Service

Paired A0 ratio

Mean ρ

Mean dP

Table S12: All five R256 Adam terminal responses. Full-reference rank, capped comparison rank, and actual output rank are distinct. The relative boundary gap is the stored gap divided by the capped exact objective; a positive absolute gap need not be large relatively. State correspondence applies to these five delivered responses.

Seed

rref

qref

qout

ρ

θmax (deg)

dP

Gap/J ∗

6200 6201 6202 6203 6204

16 15 13 14 14

8 8 8 8 8

8 8 8 8 8

0.996425 0.995430 0.994618 0.992189 0.995622

84.109 88.699 86.733 87.106 89.407

1.8307 1.9698 2.4510 2.4293 2.1019

0.000446 0.000497 0.000192 0.000382 0.000342

41