Disentangling Interpretable Factors with Supervised Independent Subspace Principal Component Analysis

Authors: jiayu  su
david a knowles
raul  rabadan
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 20:43:17

Abstract: The success of machine learning models relies heavily on effectively representing high-dimensional data. However, ensuring data representations capture human-understandable concepts remains difficult, often requiring the incorporation of prior knowledge and decomposition of data into multiple subspaces. Traditional linear methods fall short in modeling more than one space, while more expressive deep learning approaches lack interpretability. Here, we introduce Supervised Independent Subspace Principal Component Analysis ($\texttt{sisPCA}$), a PCA extension designed for multi-subspace learning. Leveraging the Hilbert-Schmidt Independence Criterion (HSIC), $\texttt{sisPCA}$ incorporates supervision and simultaneously ensures subspace disentanglement. We demonstrate $\texttt{sisPCA}$'s connections with autoencoders and regularized linear regression and showcase its ability to identify and separate hidden data structures through extensive applications, including breast cancer diagnosis from image features, learning aging-associated DNA methylation changes, and single-cell analysis of malaria infection. Our results reveal distinct functional pathways associated with malaria colonization, underscoring the essentiality of explainable representation in high-dimensional data analysis.

Main paper: 
 introduction: 
 High-dimensional data generated by complex biological mechanisms encapsulate an ensemble of patterns. A prime example is single-cell RNA sequencing (scRNA-seq) data (Fig. 1 ). These datasets measure the expression of tens of thousands of genes across potentially millions of cells, creating rich tapestries woven from interacting cellular pathways, gene dynamics, cell states, and inherent measurement noise. To unravel these patterns and reveal hidden relationships, it is necessary to decompose the data into meaningful, lower-dimensional subspaces.Linear representation learning methods, such as Principal Component Analysis (PCA) [Hotelling, 1933] and Independent Component Analysis (ICA) [Comon, 1994] , extract latent spaces from data using explainable linear transformations. These widely employed unsupervised tools learn a single latent space or a union of one-dimensional subspaces. Independent Subspace Analysis (ISA) extends ICA by extracting multidimensional components as independent subspaces [Cardoso, 1998 ]. Yet, the unsupervised nature of these methods precludes knowledge integration, restricting their utility and sometimes even identifiability [Theis, 2006] . In the context of Fig. 1 , the representation learned without supervision fails to separate temporal variability from technical batch effects.Conversely, recent advancements in deep generative models, especially semi-supervised approaches [Kingma et al., 2014] , have shown promise in disentangling diverse latent spaces and retaining relevant information under supervision. However, challenges remain in ensuring subspace independence within a variational autoencoder (VAE). While models like β-VAE [Higgins et al., 2016] , HCV [ Lopez et al., 2018] and biolord [Piran et al., 2024] attempt to address this, inference for deep generative models remains challenging and the learned representations are not interpretable.To bridge the gap, we propose Supervised Independent Subspace Principal Component Analysis (sisPCA)1 , an innovative method extending PCA to multiple subspaces. By incorporating the Hilbert-Schmidt Independence Criterion (HSIC), sisPCA effectively decomposes data into explainable independent subspaces that align with target supervision. It thus reconciles the simplicity and clarity of linear methods with the nuanced multi-space modeling capabilities of advanced generative models. In summary, our contributions include:• A multi-subspace extension of PCA for disentangling linear latent subspaces in highdimensional data. We additionally show that supervising subspaces with a linear target kernel can be conceptualized as linear regression regularized akin to Zellner's g-prior.• An efficient eigendecomposition-based alternating optimization algorithm to compute sisPCA subspaces. The learning process with linear kernels resembles matrix factorization, which may potentially benefit from desirable local geometric properties (Conjecture 3.1). • Demonstrated effectiveness and interpretability of sisPCA in various applications. This includes identifying diagnostic image features for breast cancer, dissecting aging signatures in human DNA methylation data, and unraveling time-independent transcriptomic changes in mouse liver cells upon malaria infection.Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? Answer: [Yes] Justification: In the paper, we proposed a multi-subspace extension of PCA to disentangle interpretable subspaces of variations using self-supervision. We believe our work is highly novel and significant to practitioners using PCA as well as general audience in computational biology.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper. 
background: 
 2.1 Hilbert-Schmidt Independence Criterion (HSIC)The Hilbert-Schmidt Independence Criterion (HSIC) serves as a methodology for testing the independence of two random variables X and Y [Gretton et al., 2005a] . It operates by embedding probability distributions into reproducing kernel Hilbert spaces (RKHS) and quantifying independence by distance between the joint distribution and the product of its marginals. Specifically, HSIC builds on the cross-covariance operator C XY : G → F, which is a generalization of the cross-covariance matrix to infinite-dimensional RKHS (i.e., the feature spaces of X and Y after transformation),C XY := E XY [ϕ(X) ⊗ ψ(Y )] -E X [ϕ(X)] ⊗ E Y [ψ(Y )].Here F := span({ϕ(X)}) and G := span({ψ(Y )}) are the RKHSs with feature maps ϕ and ψ for X and Y respectively. The tensor product operator f ⊗ g maps G to F, such that (f ⊗ g)h := f ⟨g, h⟩ G for all h ∈ G. The HSIC is then defined as, where v i and u j are the orthogonal bases of G and F. X and Y are independent if and only if HSIC(X, Y ; F, G) = 0. In practice, given a finite sample {(x i , y i )} n i=1 from the joint distribution P X,Y , the empirical HSIC can be computed as,HSIC n (X, Y ) = 1 (n -1) 2 tr(K X HL Y H)where K X and L Y are matrices of kernel evaluations over the samples {x i } n i=1 in F and {y i } n i=1 in G, and H is the centering matrix defined as H = I n -1 n 1 n 1 ⊤ n , with I n being the identity matrix.The HSIC, as a non-parametric criterion for independence, has been effectively integrated into numerous representation learning models, especially for disentangling complex data structures. One of the earliest applications of HSIC is ICA [Comon, 1994] , where the goal is to recover unmixed and statistically independent sources. Gretton et al. [2005a] showed that minimizing HSIC via gradient descent outperforms specialized linear ICA algorithms. Alternatively, HSIC can also be maximized to encode specific information in learning tasks. In supervised PCA (sPCA), HSIC is deployed to guide the identification of principal subspaces with maximum dependence on target variables [Barshan et al., 2011] . More recently, Ma et al. [2020] used HSIC as a supervised learning objective for deep neural networks to bypass back-propagation, and Li et al. [2021] subsequently extend it to the context of self-supervised learning.Our primary interest lies in extending these models to identify multiple subspaces, each representing independent, meaningful signatures. In this direction, Cao et al. [2015] suggested the inclusion of HSIC as a diversity term in multi-view subspace clustering, encouraging representations from different views to capture complementary information. Building on a similar idea, Lopez et al. [2018] incorporated HSIC-based regularization into VAE architectures to promote subspace separation.However, there has yet to be a linear multi-space model for more interpretable data decomposition.The presented work is also related to contrastive representation learning, as introduced in Abid et al.[2018] and Abid and Zou [2019] . Along this line of research, recent studies have also applied HSIC to regularize contrastive VAE subspaces [Tu et al., 2024 , Qiu et al., 2023] . See Appendix A for more discussions on connections and differences.Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes] Justification: See Section 5. Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
the sispca model: 
 We introduce Supervised Independent Subspace Principal Component Analysis (sisPCA), a linear model for disentangling independent data variation (Fig. 2 ). The model's linearity ensures explicit interpretability and enables regression-based extensions such as sparse feature selection. We formally discuss the connection between sisPCA and regularized linear regression in Section 3.2.As motivated in Fig. 1 , given a dataset {X i ∈ X p } n i=1 with n observations and p features and m associated target variables {Y i ∈ Y m } n i=1 , we aim to find m separate subspace representations of the data {{Z j i ⊂ R dj } m j=1 } n i=1 with latent dimensions {d j } m j=1 . Each subspace should maximize dependence with one target variable while minimizing dependence with other subspaces. For simplicity, we assume Euclidean space, {X i } := X ∈ R n×p and {Y i } := Y ∈ R n×m . Other types of data can be easily handled using appropriate kernels, which will be discussed later. The linear projection to the j-th subspace U j : R p → R dj is represented by the matrix U j ∈ R p×dj , with Z j := XU j ∈ R n×dj being the data representation in this j-th subspace. Similar to the concept of PCA loading, the projection U j depicts linear combinations of original features in X and thus can be directly interpreted as feature importance scores. (1)subject to U T j U j = I, ∀j ∈ {1, ..., m},where K Zj and K Yj are kernels defined on the j-th subspace Z j and the j-th target variable Y j , respect. This formulation differs from kernel PCA, where in the first term of (1) the kernel is defined over Z rather than X (the kernel extension of sisPCA will be discussed later).In the special case where K Zj := Z j Z T j is linear for all subspaces, the optimization becomesEQUATIONsubject to U T j U j = I, ∀j ∈ {1, ..., m}. The first term is the supervised PCA objective [Barshan et al., 2011] . We can thus view the above formulation, termed sisPCA-linear, as an extension of supervised PCA to multiple subspaces with additional regularization for subspace independence (Fig. 2 ). It comes with an appealing property: Remark 3.1. Maximizing the sisPCA-linear objective (2) is equivalent to minimizing the reconstruction error of a linear autoencoder plus regularization. See Appendix B.We now examine the second HSIC term in (2). Consider two subspaces Z u := XU ∈ R n×du andZ v := XV ∈ R n×dv with centered X. The HSIC regularization is tr(XU U T X T XV V T X T ) = tr(Z u Z T u Z v Z T v ) = ||Z T u Z v || 2 F ≥ 0.This term equals zero if and only if Z u and Z v are orthogonal. While not convex in U and V jointly, the coupling of U and V solely through the matrix product Z T u Z v indicates a well-behaved local geometry, as explored by Sun and Luo [2016] . Indeed, it allows sisPCA-linear to benefit from theoretical insights on matrix factorization. For example, Ge et al. [2017] showed that for a quadratic function f over the matrixU T V -in our context f = ||Z T u Z v || 2F -all local minima are also globally optimal under mild conditions achievable through proper regularization.Returning to (2), we see that spurious local optima may only emerge from subspace imbalance in the first symmetry-breaking supervision term, where some subspaces may contribute more to the overall objective. This leads to the following conjecture on the optimization landscape of (2): Conjecture 3.1 (informal): The sisPCA-linear objective (2) has no spurious local optima under balanced supervision, which is achievable through proper regularization; With unbalanced supervision, the global optima can still be recovered using local search algorithms following simple initialization based on the relative supervision strength of each subspace.Intuitively, the balanced supervision condition is to ensure that the local optima induced by subspace symmetry and interchangeability are also global optima. For more discussions on the optimization landscape, the balance condition, and Conjecture 3.1, see Appendix C.We now consider an iterative optimization approach to solve (2).Remark 3.2. The optimization problem in (2) can be solved via alternating optimization, where each iteration has an analytical update for each subspace.Using basic matrix algebra, the objective of (2) simplifies to,m j=1 tr(U T j X T HK Yj HXU j ) - λ 2 m j=1 m i̸ =j tr(K Zi HK Zj H) = m j=1 tr(U T j X T Kj XU j )whereKj := H(K Yj -λ 2 m i=1,i̸ =j K Zi )H and K Zi := Z i Z T i = XU i U T i X T .Given the set {U i̸ =j } i , U j can be updated by maximizing tr(U T j X T Kj XU j ), leading to the update U (t+1) j ← Q dj , where Q dj are the columns of Q corresponding to the d j largest eigenvalues from the eigendecomposition X T Kj X := QΛQ T .The full optimization process is outlined in Algorithm 1, Appendix B. This procedure guarantees convergence to an optimum as the objective is bounded and non-decreasing in every iteration. We further implement an initialization step to find the path (subspace update order) towards the global optimum as proposed in Conjecture 3.1. Briefly, we compare subspace contributions to the supervision loss and prioritize updates for subspaces under stronger supervision.While convenient, a zero HSIC regularization loss with a linear kernel in (2) does not guarantee independent subspaces. For strict independence, the subspace kernel K Z in (1) needs to be universal (e.g., Gaussian). We refer to this as sisPCA-general and solve it using gradient descent (Algorithm 2, Appendix D). The naive implementation with a Gaussian kernel has complexity O(n 3 ) in contrast to O(n 2 ) with a linear kernel. See Appendix D for performance difference discussions.Both sisPCA-linear and sisPCA-general are linear methods where {U m } measures direct contribution of original features to each subspace. Nevertheless, the sisPCA framework can be easily extended to incorporate nonlinear feature interactions, analogous to kernel PCA. See Appendix E.The use of kernel independence measure in (1) allows sisPCA to accommodate various data types through flexible kernel choices for K X (data), K Z (latent subspace), and K Y (target).For categorical variables Y = Y j n i=1 (e.g., cancer types), we use the Dirac delta kernel,K Y (i, j) = 1 Yi=Yj .It is also possible to use other general graph kernels for categorical variables with an intrinsic hierarchical structure, e.g., subtypes and stages [Smola and Kondor, 2003] .For continuous variables Y ∈ R n , we use the linear kernelK Y = Y Y T .When Y ∈ R n×d is multivariate, the kernel is the sum of per-dimension kernelsK Y = d i=1 Y :i Y :i T .Remark 3.3. Maximizing the sisPCA-linear objective (2) with linear kernels on the target space is equivalent to performing regularized regression against the target. See Appendix B.In a nutshell, we show that sisPCA-linear can be viewed as approximating the target Y with Z = Xu. The particular regularization on u corresponds to a zero-mean multivariate Gaussian prior, which is related to Zellner's g-prior in Bayesian regression.In practical applications, a common goal would be to recover both subspaces linked to known attributes (supervised) and to unknown attributes (unsupervised) simultaneously. In sisPCA, this is achieved by setting the target kernel for the unknown subspace to the identity matrix (K Y = I), corresponding to unsupervised PCA (Fig. 2 ). Section 4.1 provides an example of this process.However, the absence of external supervision introduces a potential identifiability issue. As indicated in Appendix B eq. 4 each supervised subspace is driven by two forces to (1) align with the target and (2) capture major variations in the data. This dual objective can lead to scenarios where unknown attributes, ideally retained in the residual unsupervised subspace, being inadvertently presented in supervised subspaces. Fig. 8 in Appendix D gives an example where sisPCA-general fails. Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. Moreover, our baseline models (PCA and sPCA) are deterministic, eliminating the need for error bars or confidence intervals that typically account for randomness.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors).• It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
applications: 
 Unless otherwise specified, in the following sections sisPCA refers to sisPCA-linear, where the objective (2) is solved using Algorithm 1. For baseline comparisons, we consider linear models including PCA and sPCA. While non-linear VAE counterparts such as HCV [Lopez et al., 2018] are included for quantitative performance benchmarking in Section 4.4, they are not considered for interpretability analyses due to their inherent complexity. The rationale and details for baseline selection are provided in Appendix F.We first consider learning latent subspaces associated with known and unknown attributes using simulated data. The dataset reflects a ground truth 6-dimensional latent space, comprising three distinct 2D subspaces (Fig. 3a ): S1 with two Gaussian distributions, S2 with a noisy 2D grid and S3 with a ring structure. The defining manifold characteristics ϕ of S3 remain unknown to the model, representing the unsupervised component. These subspaces were concatenated and linearly projected to 20-dimensions using a 6 × 20 matrix with entries uniformly distributed on [0, 1].Both unsupervised PCA (targeting S3) and supervised PCA (S1 and S2) subspaces capture structures heavily influenced by S2 (Fig. 3b and Fig. 11 in Appendix H), due to S2's pronounced variations. In contrast, sisPCA markedly improves the disentanglement of these subspaces, especially the unsupervised S3 (Fig. 3c ). Despite S2's dominant influence, sisPCA isolates the effects of each subspace, resulting in clearer separation of the three independent signals. The two supervised subspaces S1 and S2 are distinctly characterized by patterns exclusively associated with the supervision attributes. In the unsupervised subspace S3, although Principal Component 2 (PC2) picks up some categorical information from S1, sisPCA successfully uncovers the underlying circular structure.We apply sisPCA to the Kaggle Breast Cancer Wisconsin Data2 to demonstrate its utility in data compression and feature extraction. The dataset contains 569 samples with 30 summary features from breast mass imaging. Our goals are to (1) learn compressed subspaces for predicting disease status ('Malignant' or 'Benign', agnostic during training), and (2) understand the relationship between original features and how they contribute to the learned representation and diagnosis potential. The diagnosis label, invisible to all models, is used to measure subspace quality via the mean silhouette score1 n i s(i) = 1 n i (b(i) -a(i))/ max{a(i), b(i)}, where a(i) and b(i) are mean intra-cluster and nearest-cluster distances for sample i. Higher scores indicate larger diagnostic potential. In addition, we use the Geodesic Grassmann distanced(Z i , Z j ) = k n=1 θ 2 n 1/2to measure subspace separateness, where k = min{dim Z i , dim Z j } and {θ n } the principal angles [Miao and Ben-Israel, 1992] . Higher scores indicate better disentanglement. Inputs for all model are zero-centered and variance-standardized.In the PCA space, samples are well-separated based on diagnosis along PC1 (Fig. 4a ). 'symme-try_mean' and 'radius_mean' are the top two features negatively contributing to PC1, motivating us to construct separate subspaces to reflect nuclei size (using 'radius_mean' and 'radius_sd' as targets) and shape (using 'symmetry_mean' and 'symmetry_sd' as targets). The remaining 26 features are projected onto these subspaces using sPCA (Fig. 4b ) and sisPCA (Fig. 4c ). In sPCA, both subspaces better explain diagnosis status than PCA but remain highly entangled (Grassmann distance: 1.493, Pearson correlation of PC2 loading: 0.850). However, with sisPCA's explicit disentanglement, the symmetry subspace loses its predictive power as the two spaces separate further (Grassmann distance: 2.710). sisPCA subspaces are constructed from distinct feature sets (PC2 loading correlation -0.203), with 'area' and 'perimeter' contributing more to the radius subspace and 'compactness' and 'smoothness' to the symmetry one. The radius subspace also gains additional predictive power by repulsing further from the symmetry space (Silhouette scores: 0.516 in sisPCA, 0.470 in sPCA).Our sisPCA results suggest that cell nuclear size is more informative for breast cancer diagnosis than nuclear shape. We confirm this by measuring directly the predictive potential of target variables (Silhouette scores: 0.457 for 'radius_mean ' and 'radius_sd', 0.092 for 'symmetry_mean' and 'sym-metry_sd') . Our conclusion also aligns with previous clinical observations [ Kashyap et al., 2018] . In contrast, PCA and sPCA, while capable to extract new diagnostic features (PC1), cannot faithfully capture feature relationships without disentanglement and potentially overestimate symmetry-related features' relevance in malignancy.Tumorigenesis and aging are two intricately linked biological processes, resulting in cancer omics data that often display patterns of both. DNA methylation (DNAm) exemplifies this complexity, undergoing genome-wide alterations during aging while also exhibiting cancer-specific changes in particular regions, presumably silencing tumor suppressor genes or activating oncogenes.The Cancer Genome Atlas (TCGA)3 offers a comprehensive collection of DNAm datasets from patients with various cancer types. In TCGA DNAm data, methylation status is probed across genomic locations using the Illumina Infinium 450K array and quantified as beta values ranging from 0 (unmethylated) to 1 (methylated). The resulting data matrix X presents challenges due to the use of three different probe types and highly correlated features, typically requiring careful preprocessing. For illustration purpose, we use a downsampled dataset comprising the first 5,000 non-constant and non-NA CpG sites from 9,725 TCGA tumor samples across 33 cancer types.Our goal is to disentangle tumorigenic signatures from age-dependent methylation dynamics in this pan-cancer DNAm data. Traditional methylation analyses often employ regression-based methods to learn site-specific statistics, later aggregated by gene or high-level genomic annotations [Bock, 2012] . However, these approaches may suffer from high dimensionality and multi-collinearity due to potential redundancy in CpG methylation activity. We propose using sisPCA to address these limitations, which allows us to (1) learn compressed, low-dimensional representations that retain biological information, and (2) minimize confounding factors not controlled in simple regression models by enforcing disentanglement. Specifically, we aim to learn two subspaces: one aligning with chronological age (CA, aging subspace, supervised with a linear kernel) and another with TCGA cancer categories (cancer subspace, supervised with a delta kernel). We evaluate the quality of learned representations using information density measured by the Silhouette score and subspace separateness by the Grassmann distance. For the rank-one aging subspace (rank(K Y ) = rank(Y Y T ) = 1), we also measure information density using the maximum absolute Spearman correlation, max d∈[1,dj ] {|ρ(Z (d) j , CA)|}, between CA and each axis of the subspace Z j ∈ R n×dj .We first validate the disentanglement effect of formulation (2). Unlike models such as Lopez et al. [2018] that use a Gaussian kernel, sisPCA minimizes the HSIC regularization with a linear kernel. This approach trades strict statistical guarantees for improved computational efficiency (detailed in Appendix D). Our experiments demonstrate that as λ increases, the two subspaces show increasing divergence (Table 1 ). Notably, while HSIC-Gaussian is not explicitly optimized, it decreases in tandem with HSIC-Linear. The generally larger values of HSIC-linear potentially offer advantages in optimization and help mitigate numerical rounding errors.Table 1 : Separateness of the aging and cancer subspaces inferred by sisPCA.sisPCAPCA λ = 0 (sPCA) λ = 1 λ = 10HSIC-Linear (in the objective of sisPCA) 481.6 189.7 2.0e-4 2.4e-05 HSIC-Gaussian 1.1e-2 6.5e-3 7.0e-4 7.0e-4 Grassmann distance 0 3.09 4.97 4.97Furthermore, the separation of aging and cancer subspaces leads to a moderate increase in target information density and a decrease in confounding information (Table 2 ). However, stronger regularization does not always equate to better representations. This is partly due to the inherent coupling between aging and tumorigenesis. Efforts to remove aging signals inevitably result in information loss on cancer type. We discuss the tuning of λ more generally in Appendix G. In an optimal pair of subspaces, each property (infection status or time) should be more distinguishable in its corresponding subspace while showing less separation in the other.We analyze scRNA-seq data of mouse hepatocytes from Afriat et al. [2022] 4 . Our goal is to distinguish genes directly involved in parasite harboring from those associated with broader temporal changes post-infection. The processed dataset comprises gene expression profiles of 19,053 cells collected at five post-infection time points (2, 12, 24, 30, and 36h) and from control mice. Infection status was determined based on GFP expression linked to malaria. We keep the top 2,000 highly variable genes and use normalized expression as model inputs. Time points are treated as discrete categories to account for potential nonlinear dynamics. See Appendix F for full experiment details.UMAP visualizations of PCA, sPCA, and sisPCA subspaces show that while PCA primarily captures temporal variations (Fig. 5a ), both sPCA and sisPCA successfully differentiate between infected and uninfected hepatocytes in their infection subspaces. However, sPCA's infection space still exhibits significant temporal effects, suggesting uncontrolled confounding effects (Fig. 12b in Appendix H). sisPCA effectively eliminates this intermingling, yielding cleaner representations where relevant biological information is further enriched in the corresponding spaces (Fig. 5b ).Comparisons with non-linear VAE counterparts (Fig. 5c and Fig. 12 , full model description in Appendix F) and quantitative evaluations (Table 3 ) demonstrate that our HSIC-based supervision formulation (1) achieves performance comparable to neural network predictors. Notably, sisPCA outperforms HSIC-constrained supervised VAE (hsVAE) [Lopez et al., 2018] in separating infected and uninfected cells. Indeed, hsVAE's performance in the infection subspace is so poor that even under supervision, the representation contains near-random information (Fig. 5c ). To address this gap, we developed hsVAE-sc by incorporating additional domain-specific knowledge (See Appendix F, and Fig. 12d in Appendix H). This model learns a much improved infection space (Silhouette score: 0.233) while maintaining high distinguishability in the temporal space (Silhouette score: 0.634). The improvement is likely due to the fact that parasite encapsulation can induce a slight increase in total RNA counts; by using unnormalized count-level data and explicitly modeling the library size, the model captures this additional information and thus produces enhanced results.While visualizations and quantitative metrics provide useful estimates, they may not fully capture the biological relevance of the representations. For instance, a subspace with two point masses perfectly separating cells by infection status would have the highest information density, yet offer little new biological insight. Advantageously, linear models like sisPCA are inherently interpretable. Using the learned sisPCA projection U as feature importance scores, we rank genes based on their PC contributions. Chemokine genes, including Cxcl10, emerge as top positive contributors to infection-PC1, which is elevated in infected cells. This aligns with their established role as acute- Figure 6 : GO biological process enrichment results of top genes contributing to the sisPCA-infection subspace. Genes are ranked by their PC1 loading and are grouped by effect direction. 
discussion: 
 This study presents sisPCA, a novel extension of PCA for disentangling multiple subspaces. We showcase its capability in interpretable analyses for learning complex biological patterns, such as aging dynamics in methylation and transcriptomic changes during Plasmodium infection. To enhance usability, we have implemented an automatic hyperparameter tuning pipeline for λ using grid search and spectral clustering, similar to contrastive PCA [Abid et al., 2018] (Appendix G).Still, sisPCA has several limitations: Linearity constraints: The linear nature of sisPCA may miss non-linear feature interactions, potentially underperforming on more complicated datasets. While nonlinear extensions are possible (Appendix E), they come at the cost of reduced computational efficiency and interpretability. Linear kernel HSIC: The HSIC-linear regularization, while computationally convenient, does not guarantee complete subspace independence. However, minimizing HSIC-linear tends to reduce HSIC with a Gaussian kernel (Table 1 ), which suggests that the issue is less of a concern in practice. Subspace identifiability: Our formulation (1) relies on external supervision to differentiate subspaces, which could lead to identifiability issues if the supervisions are too similar, or when one subspace is unsupervised.Despite these limitations, sisPCA's ability to provide interpretable, disentangled representations of complex biological data makes it a valuable addition to the toolkit of biologists working with highdimensional datasets. We envision future applications on larger-scale omics datasets and potential new biomedical discoveries.Contrastive representation learning, as introduced in Abid et al. [2018] , aims to find representations that capture the difference between explicit positive and negative samples. For example, contrastive PCA (cPCA) takes the pair of a target dataset of interest and a background dataset as inputs and finds the principal components that account for large variances in the target but small variances in the background [Abid et al., 2018] . cVAE [Abid and Zou, 2019] and contrastiveVI [Weinberger et al., 2023] apply the idea to VAE architectures, where the target and background data are assumed to have different generative processes and the corresponding representations lie in separate subspaces.The learning of the target and background subspaces are achieved by explicitly forcing target and background samples to use different spaces. More recently, Tu et al. [2024] and Qiu et al. [2023] adopt an HSIC-based regularization approach to encourage the disentanglement of the two subspaces.Despite the similarity, our proposed sisPCA model is conceptually different from these contrastive models. sisPCA (and the non-linear HCV [Lopez et al., 2018] ) focuses on the decomposition of a single dataset, as opposed to a pair of target and background datasets. The sense of contrast in sisPCA comes from explicit supervision, of which each subspace either aligns with or is repulsed from. If we concatenate the pair of target-background into one dataset and consider the case-control information as supervision, the associated sisPCA subspace can also reflect the difference similar to cPCA (the objective is still not the same).B Properties of the sisPCA-linear optimization problem B.1 Proof of Remark 3.1 on the equivalence of sisPCA-linear and autoencoderHere we show that thesisPCA-linear optimization 2 is equivalent to solving a regularized linear autoencoder. First note the well-known fact of PCA that maximizing the variance of the projection ||XU || 2 F is equivalent to minimizing the reconstruction error of a linear orthogonal autoencoder, since||X-XU U T || 2 F = tr(X T X -2X T XU U T + X T U U T U U T X) = tr(X T X) -tr(U T X T XU ) = ||X|| 2 F -||XU || 2 F .Now consider supervised PCA with a positive semidefinite kernel K y and its Cholesky decomposition K y := L T L. Note L := I in vanilla unsupervised PCA. That is, PCA can be viewed as aligning to a supervision kernel where every sample forms a group of itself. Without loss of generality, we assume X and K y to be centered, such that the first term in eq.2 per subspace becomestr(XU U T X T K y ) = tr(U T X T K y XU ) = ||LXU || 2 F .Maximizing the above quantity is thus equivalent to minimizing the weighted reconstruction error ||L(X -XU U T )|| 2 F where L from supervision controls the aspects to focus on.Below we describe the Algorithm 1 to solve (2).Here we show that the sisPCA-linear problem 2 is equivalent to a regularized linear regression. For simplicity, we assume both X and Y to be centered and Y ∈ R n×1 to be univariate. Since X T K Y X = X T Y Y T X is rank-one, we need only to consider the largest eigenvalue and corresponding eigenvector of X T K Y X. Therefore, maximizing the first cross-covariance termEQUATIONAlgorithm 1 Solving sisPCA-linear using alternating eigendecomposition input: Data X ∈ R n×p , kernels for m target variables{K Yj , j = 1, ..., m}, dimensions of subspaces {d j , j = 1, ..., m}, independence regularization λ.Randomly initialize and orthonormalize U j and, optionally, determine the optimal order of subspace updates through one complete update cycle. Assuming with update order {1, ..., m}. repeatK Zi ← XU i U T i X T for i ∈ {1, ..., m}; for j ∈ {1, ..., m} do Kj ← H(K Yj -λ 2 m i=1,i̸ =j K Zi )H; Update U j ← the first d j eigenvectors of X T Kj X;Update K Zj ← XU j U T j X T ; end for until converge return: Linear projections {U j } and subspace representations of data {Z j := XU j } for i ∈ {1, ..., m}.That is, maximizing the linear HSIC(Xu, Y ) solves the regression problem of approximating the target Y with u serving as the regression coefficients, while also maximizing variations of the fitted value Ŷ := Xu (or, equivalently, minimizing an additional self-supervised reconstruction loss). The particular regularization term in (3) corresponds to a zero-mean multivariate Gaussian (MVN) prior on u. To see this, we first diagonalize X T X into P T SSP where P is orthogonal and S = diag{σ 1 , ..., σ p } contains the singular values of X in decreasing order. Denote k := P u ∈ R p×1 . Now we have-||Xu|| 2 = -(P u) T S(P u) = - p i=1 σ 2 i k 2 i = p i=2 (σ 2 1 -σ 2 i )k 2 i -σ 2 1 ||k|| 2 ||k||=1 == || Xu|| 2 -σ 2 1where X := (σ 1 I -S)P . Subsequently, (3) becomesEQUATIONand solving (5) is equivalent to finding the maximum a posteriori estimator of the regression coefficients u under a MVN prior with zero mean and inverse covariance XT X. Directions corresponding to smaller singular values will receive stronger shrinkage effects. The formulation of ( 5) is also closely related to Zellner's g-prior in Bayesian regression where the coefficient β follows β ∼ M V N (0, gσ 2 (X T X) -1 ) for some positive g and noise variance σ 2 . Finally, when fixing other subspaces {Z v }, HSIC(Xu, Z v ) := (Xu) T K Zv (Xu), now a function of u, can also be integrated into (5). The new quadratic regularization term is Q(u) := u T (X T K Zv X + XT X)u, which again corresponds to a zero-centered MVN prior.C Optimization landscape of sisPCA-linear and Conjecture 3.1We first provide a motivating example to see how the first supervision term in (2) affects the optimization landscape. Assume no supervision Y (i.e. K Ym = I) and all the latent subspaces {Z m } share the same dimension d. In the absence of disentanglement regularization, every subspace will simply be the same vanilla PCA space, leading to complete overlap. The introduction of the disentanglement penalty in (2) maintains the objective's symmetry with respect to U m . In other words, subspaces are interchangeable in the pure unsupervised setting, and consequently the symmetry gives rise to multiple local optima. Nevertheless, the presented challenge is trivial, in that all local optima are also global, and the number of local optima, considered in terms of invariant groups up to column sign-flipping, depends only on the number of subspaces m.Under supervision, some subspaces may be selectively favored. For example, the scaling of subspace supervision kernels {K Ym } can break the symmetry of the overall objective, allowing information to be preferentially preserved in one subspace while being removed from others. In such scenarios (referred to as unbalanced supervision), local optima that are closer to the supervised PCA results of subspaces with larger K Ym also yield higher objective values. The balanced supervision condition can be evaluated by calculating the relative ratio of kernel scales, or equivalently, the ratio of objective gradient with respect to each subspace kernel.Another important observation is that, because the disentanglement penalty is symmetric, tuning the hyperparameter λ will not alter the ranking of local optima. This suggests a straightforward initialization strategy to navigate towards the global solution by evaluating the relative supervision strength, HSIC(Z 0 m , Y m ), for each subspace. Here, Z 0 m represents the initial, pre-disentanglement solution for subspace m, which can be efficiently computed using supervised PCA.We can directly visualize the optimization landscape in a simplified scenario (Fig. 7 ). Consider X = 1 0 -1 0 , which consists of two 2D vectors separated along the first axis. With supervisionsY 1 = (1, -1)T and Y 2 = cY 1 differing only by a scalar c, our goal is to determine the projections U = (u, u 2 ) T and V = (v, v 2 ) T that transform X into subspace representations Z 1 = XU = (u, -u) T and Z 2 = XV = (v, -v) T . Adopting a linear kernel for K Y , the sisPCA-linear objective ( 2) is now equivalent to the following simplified problemargmax u,v∈[0,1] f (u, v) = u + cv -λuv.Since X only has one effective dimension, both U and V align with that axis (i.e., u = 1 and v = 1, the unregularized solution) when the disentanglement penalty λ is low (λ < 1), leading to a single optimal solution (u, v) = (1, 1). Under strong regularization, two local maxima (1, 0) and (0, 1) emerge, and their relative order depends solely on the scalar c and is unaffected by λ. Specifically, each maximum represents a preference for capturing the significant axis in X within one subspace, allowing for an easy decision on initialization based on a comparison between u and cv. In the context of alternating optimization methods like Algorithm 1, it can be simply done by fixing the subspace V (v = 1) and updating U first to minimize the disentanglement penalty (from u = 1, the unregularized solution, to u = 0). The optimization path stays near the global optima regardless of λ. This phenomenon also resembles the situation described in Ge et al. [2017] on f (U V T ), where an additional regularizer ||U T U -V T V || 2 F is necessary to maintain balance when dealing with asymmetric matrices. The first half of Conjecture 3.1 is indeed motivated by the above result. For sisPCA, a similar regularizer in the form of||Z i Z T i -Z i Z T i || 2F may be enough to ensure balanced supervision. Moreover, given that the multi-solution challenge arises from subspace interchangeability, we hypothesize that more discriminative Y could alleviate the issue. Since there are only pairwise interactions between subspaces in (2), Ge et al. [2017] 's framework can be adapted to study the behavior of more than two subspaces. Intuitively, the proof of Conjecture 3.1 would be to show that the local gradient direction of (2) aligns with the global descent direction in a region near global optima, and proper regularization or initialization can make sure the optimizer stays in that region. Following the procedures described in section 5 of Ge et al. [2017] for asymmetric matrices, we may concatenate all subspaces into W := [Z 1 , ..., Z m ] T and reformulate the objective as a quadratic function of N := W W T plus some regularization.In this section, we discuss the sisPCA-general optimization problem in detail. Recall that it is motivated by the limitation of sisPCA-linear that zero HSIC disentanglement regularization does not ensure statistical independence. In sisPCA-general, we want to solve (1) by employing HSIC with with universal kernels, e.g. the Gaussian kernel. Definition D.1. Given a compact metric space X and a continuous kernel K on X , K is called universal if for every continuous function g ∈ C(X ) and all ϵ > 0, there exists an f ∈ H such that ||f -g|| X < ϵ.A nice result on universal kernel from the Constrained Covariance (COCO) framework [Gretton et al., 2005b] is that it can approximate any continuous function on a compact metric space. In other words, the HSIC with a universal kernel will be zero only if, among all possible continuous transformation, the covariance of two random variables are always zero. In sisPCA-general, the latent representation is still encoded by a linear projection Z j = XU j . This indicates that U j can be interpreted in the same way as PCA and sisPCA-linear to extract feature importance.We solve the sisPCA-general problem with Gaussian kernel,k(x, x ′ ) = exp(-w 2 ||x -x ′ || 2 ),using gradient descent (Algorithm 2). Since the problem is constrained, the algorithm unfortunately is not guaranteed to converge. The HSIC with Gaussian kernel also affects the optimization landscape. Empirically, we observed that sisPCA-general tends to have many local optima, requiring exhaustive tuning of training parameters such as the learning rate and stopping criteria. The Gaussian kernel scale w is also important. Ma et al. [2020] reported in their experiments that the performance of HSIC-based optimization is moderately sensitive to w and thereby proposed a multiple-scale solution that essentially aggregates results from different w. In our experiments, we followed Gretton et al. [2005a] by setting the kernel scale to the median distance between data points in the input space.Another notable disadvantage of sisPCA-general is its cubic-scaling computational complexity.Here we focus on scenarios where the number of features p is fixed and the sample size n is large (n > p), which is common when data preprocessing and filtering is in place. Specifically, the computational complexity for calculating HSIC(X, Y ) with a linear kernel for X ∈ R n×p and a given target kernel K Y for Y is O(n 2 p). If both K X and K Y are linear, the complexity can be further reduced to O(np). Additionally, the eigendecomposition of the (cross) covariance incurs a cost of O(p 3 ). This means that each update in Algorithm 1 has an overall complexity of O(n 2 ) wrt the sample size n, aligning with that of supervised PCA. In comparison, most dimensionality reduction methods requiring pair-wise distance computation have a complexity of at least O(n 2 ). On the other hand, for sisPCA-general, computing a Gaussian-kernel HSIC(X, Y ) naively is of O(n 3 ) complexity. From an engineering perspective, further optimizations for both the sisPCA-linear (Algorithm 1) and sisPCA-general (Algorithm 2) algorithms are conceivable, for instance, by exploiting data sparsity and by mini-batch or incremental update.Algorithm 2 Solving sisPCA-general using gradient descent input: Data X ∈ R n×p , kernels for m target variables {K Yj , j = 1, ..., m}, dimensions of subspaces {d j , j = 1, ..., m}, a universal kernel Φ Z (•, •) for latent representations, independence regularization λ. Randomly initialize and orthonormalize U j and, optionally, determine the optimal order of subspace updates through one complete update cycle. Assuming with update order {1, ..., m}. repeatK Zi (k, r) ← Φ Z (X k U i , X r U i ) for i ∈ {1, ..., m}; L ← - m j=1 tr(K Zj HK Yj H); L ← L + λ m j=1 m i>j tr(K Zi HK Zj H);Jointly update {U j } by minimizing L using gradient descent; Project U j to the orthonormal space; until converge return: Linear projections {U j } and subspace representations of data {Z j := XU j } for i ∈ {1, ..., m}.Recovering supervised and unsupervised subspaces in simulated data We demonstrate the performance and limitation of sisPCA-general using simulated data. The task and simulation details are described in Section 4.1. As illustrated in Fig. 8 , sisPCA-general partially recovers the two supervised subspaces S1 and S2. Challenges arise, however, in fully eliminating continuous S2 signatures from S1. This difficulty may stem from the smaller scale of HSIC values with Gaussian kernels as compared to those with linear kernels, and a more complicated optimization landscape due to the non-linearity introduced in the regularization. Additionally, the unsupervised subspace S3 appears to collapse to one dimension, which partially illustrate the identifiability issue raised in Section 3.3. The issue may be more pronounced in sisPCA-general because of the diminished HSIC scale (i.e. weaker connection between subspaces and their target variables). Following kernel PCA's methodology [Schölkopf et al., 1997] , we can extend sisPCA for nonlinear scenarios where features in X interact through a feature map Φ : R p → F . The space F may have infinite dimension, and we can only explicitly write out the projection U in the finite case F ⊂ R k , for example, Φ(x) = x 2 . In that context, the latent representation Z := Φ(X)U can be naively computed, and the sisPCA objective (1) will stay the same. By simply replacing X with Φ(X), methods such as gradient descent and Algorithm 1 still hold. However, the beauty of the kernel trick is that we can bypass the computation of Φ(X) to obtain Z. Instead of learning the projection matrix U , which is no longer possible when F is infinite-dimensional, we can learn a set of coefficients ⃗ α fromZ = Φ(X)U := n i ⃗ α i k(X, X i ),with {X i } representing the n data points and k(x, y) := ⟨Φ(x), Φ(y)⟩ the corresponding kernel of Φ. Note that k(•, •) is specified wrt to the data X, which is different from the kernel K z for latent space and the kernel K y for supervision. We leave the practical implementation of this extension to future exploration. We evaluate model performance based on two key aspects:Representation Quality: This criterion assesses the ability to learn low-dimensional representations that reflect specific and unmixed data properties. A high-quality subspace representation should contain minimal confounding information from other subspaces. We employ both quantitative metrics (e.g., information density and subspace separateness) and qualitative 2D visualizations of learned representations.Interpretability: This criterion evaluates the ease of understanding feature contributions to different data properties. For linear methods (PCA, sPCA, and sisPCA), we consider the learned projection U as feature importance and extract top features from the loading (see examples in Sections 4.2 and 4.4). Non-linear models generally lack straightforward feature-to-subspace mapping, making them inherently less interpretable. While approaches like gradient-based saliency maps are becoming standard for interpreting black-box models such as VAEs, they are less straightforward and face challenges in aggregating sample-level gradients into global feature importance scores. Therefore, we only compare sisPCA to other linear models for interpretability analysis (e.g. ranking and selecting top genes according to the loading).Table 4 provides a comprehensive comparison of linear and non-linear models used in this study, including the linear PCA, sPCA and sisPCA (proposed work) described in the Methods section, as well as sevral non-linear VAE counterparts inspired by the work of HCV [Lopez et al., 2018 ]:• VAE: Vanilla VAE with Gaussian likelihood.• supVAE: Gaussian VAE with additional predictors for target variables.• hsVAE: supVAE with additional HSIC disentanglement penalty (Gaussian kernel). It is essentially a general-purpose HCV [Lopez et al., 2018] with Gaussian likelihood.In the sisPCA package 5 , PCA and sPCA are implemented as special cases of sisPCA and solved analytically using eigendecomposition. VAE models and variational inference are implemented using the SCVI framework 6 . The VAE model is a special case of SCVI's VAE module (latent_distribution = "normal", dispersion = "gene") but with a Gaussian generative model. The supVAE is a multi-subspace extension of VAE with additional predictors to predict target supervisions from the corresponding latent representations. The last layer of the predictor is either a linear transformation for continuous supervision or linear followed by softmax for categorical supervision. The hsVAE model is an extension of supVAE where we add additional subspace disentanglement penalty (i.e., HSIC with a Gaussian kernel) to the objective, as described in HCV [Lopez et al., 2018] .For single-cell RNA-seq data, we further implement a domain-specific hsVAE model, namely hsVAEsc. Specifically, hsVAE-sc uses negative binomial as the data generative distribution with an extra autoencoder to learn the library size (the sum of total counts per cell) and is correspondingly trained on counts-level data. It is very close to a re-implementation of the scVIGenQCModel developed in the HCV paper under the latest SCVI framework.As outlined in Table 4 , sisPCA has three main hyperparameters:1. Subspace latent dimension d.2. Subspace supervision kernel K Y .3. Disentanglement penalty strength λ.The tuning of (3) will be discussed in Appendix G. Here we illustrate the selection of ( 1) and ( 2), which are largely pre-determined based on the target variable supplied for supervision.In our analysis, we consistently use the delta kernel for categorical targets and the linear kernel for continuous targets (Section 3.2). Given a supervision kernel K Y of rank k, the effective dimension of the corresponding sisPCA subspace, defined as the number of dimensions with non-zero variance, is determined by the eigendecomposition step in Algorithm 1. Specifically, we have the following upper bounds:• For sPCA (λ = 0): The effective dimension equals rank(X T K Y X) ≤ min(k, p) ≈ k,where p is the number of features in X.• For sisPCA: The effective dimension equals rank(X T KY X) ≤ min(k + d, p) ≈ k + d,where KY ← K Y -λZ T Z is the supervision kernel after the rank-d disentanglement update and d is the number of latent dimensions of Z.In practice, we observe that the effective dimensions of both sPCA and sisPCA closely approximate the target kernel rank k. For instance:• For a 1-D continuous variable with linear kernel (e.g., age in Section 4.3), the effective dimension is always 1. • For a categorical variable of k groups with delta kernel (e.g., time point in Section 4.4), the effective dimension is approximately k.Consequently, we find that the learned sisPCA subspaces remain consistent regardless of the specified dimension. Any additional axes beyond the effective dimension tend to collapse to zero and are primarily affected by numeric errors in the SVD solver. Since sisPCA orders subspace dimensions by the variance explained (i.e., eigenvalues), it is straightforward to remove these extra dimensions by examining the variance explained curve.In contrast, VAE models are more sensitive to dimension changes and generally have more hyperparameters to tune (Table 4 ). Designing optimal architectures for autoencoders and target predictors is a topic in itself and beyond the scope of our current work. For benchmarking purposes, we largely follow the default VAE architecture from SCVI:• Autoencoders for latent mean and variance: One hidden layer with 128 hidden units, ReLU activation, batch normalization, and dropout. • Predictor design (adapted from the scVIGenQCModel in the HCV paper): One hidden layer with 25 neurons, ReLU activation, and dropout. Extra softmax output head for classifier. • Training objective: Equal weighting of prediction error and reconstruction loss, minimized using variational inference.In Section 4.4, we set the subspace dimension to 10 for all VAE models. We use UMAP to project the learned subspace onto 2D and subsequently calculate the Silhouette score in Table 3 on this 2D subspace, ensuring dimensional consistency with the PCA-based spaces (which are also projected onto 2D using UMAP).Recovering Supervised and Unsupervised Subspaces in Simulated Data The code used to simulate the donut data in Section 4.1 is available in the sisPCA GitHub repository. We have verified that altering the noise distribution (e.g., from Gaussian to uniform noise, or from uniform projection matrix to Gaussian matrix) yields similar results.Learning Diagnostic Subspaces from Breast Cancer Image Features We obtained the dataset from Kaggle and scaled each real-valued feature to have zero mean and unit variance. The input for PCA includes all 30 quantitative features. The 'radius' subspace is supervised with 'radius_mean' and 'radius_sd' (effective dim = 2), while the 'symmetry' space is supervised with 'symmetry_mean' and 'symmetry_sd' (effective dim = 2). For sPCA and sisPCA, we use the remaining 26 features as input and project the data onto each subspace. Using the elbow method, we set the PCA subspace dimension to 6, and correspondingly, the dimension for sPCA and sisPCA subspaces to 3 (although their effective dimension is approximately 2). The Silhouette scores are computed on the first three principal components for all subspaces.Separating aging-dependent DNA methylation changes from tumorigenic signatures We downloaded the TCGA methylation data from GDC (controlled access) and retained only the first 5,000 non-constant and non-NA CpG sites for analysis. The data was zero-centered. All subspaces are specified with 10 latent dimensions. However, as noted in the previous section, the sPCA and sisPCA aging subspaces have an effective dimension of only one.Disentangling infection-induced changes in the mouse single-cell atlas of the Plasmodium liver stage We used the single-cell data preprocessed by Piran et al. [2024] 7 . The processed data comprises 19,053 cells and 8,203 mouse genes (all malaria genes have been filtered out), and we keep the top 2,000 variable genes using scanpy's 'sc.pp.highly_variable_genes' function. For all models except hsVAE-sc, we use the normalized and log1p-transformed expression data as inputs. For hsVAE-sc, we supply the raw data before normalization. It's important to note that the data underwent a background correction step before normalization, where the mean expression in empty wells was subtracted from the observed data. This results in the raw counts being technically non-integers. Consequently, in hsVAE-sc, a continuous extension of the negative binomial likelihood is used to compute the reconstruction loss. All subspaces are given 10 latent dimensions. Although based on the kernel rank, we know that the sPCA and sisPCA infection subspaces have an effective dimension of 2, while the time subspaces have an effective dimension of 6.We adapt the approach from Abid et al. [2018] to develop a computational pipeline for automated selection of the sisPCA hyperparameter λ. The process involves:1. Fitting sisPCA models with multiple λ values. 2. Computing pairwise similarity (affinity) between learned subspaces. 3. Identifying representative subspace clusters using spectral clustering.We compute subspace affinity based on the principal angle [Miao and Ben-Israel, 1992] :d subspace (Z 1 , Z 2 ) := d({θ 1 , ..., θ k }) := d(Θ)where Z ∈ R n×k is the orthonormal basis of a dim-k subspace with n observations, and the principal angles Θ are computed from the SVD of Z T 1 Z 2 = U cos(Θ)V T . Through empirical testing, we found the Fubini-Study Grassmann distance provides optimal performance. Thus, we define subspace affinity as:Affinity subspace (Z 1 , Z 2 ) = k i=1 cos(θ i ) = k i=1 σ i (Z T 1 Z 2 ).7 See the notebook https://github.com/nitzanlab/biolord_reproducibility/blob/main/notebooks/spatiotemporal-infection/1_spatio-temporal-infection_preprocessing.ipynb For stability, we truncate subspace latent dimension to k ′ := min(rank(K Y ), k), where K Y is the corresponding supervision kernel (see discussion on subspace effective dimension in Appendix F.3). For sisPCA models with multiple subspaces, we compute model-wise affinity as the average of subspace affinitiesAffinity model ({A 1 , ..., A m }, {B 1 , ..., B m }) = 1 m m i=1Affinity subspace (A i , B i ).Learning Diagnostic Subspaces from Breast Cancer Image Features We evaluated 20 λ values ranging from 0 to 100 and analyzed the pairwise similarity between learned sisPCA subspaces (Fig. 9 ). The sPCA solution (λ = 0) shows clear separation from other sisPCA models. As λ increases, the symmetry subspace becomes progressively less predictive of diagnostic status (Fig. 9b ), and subspaces stabilize and converge to a robust solution after λ = 1. This convergence pattern is also reflected in the elbow of the reconstruction loss curve (Fig. 9c ). Disentangling infection-induced changes in the mouse single-cell atlas of the Plasmodium liver stage We examined 10 λ values ranging from 0 to 100 to assess their effect on subspace characteristics (Fig. 10 ). Our analysis again reveals the clear separation of sPCA solution (λ = 0) from other sisPCA models. As λ increases, infected cells form a tighter cluster in the infection subspace (Fig. 10a ), while temporal dynamics become less pronounced, with increased mixing of cells from different time points (Fig. 10b ). Moreover, sisPCA results are quite robust across λ in preserving subspace structure (Fig. 10c ). Since the disentanglement effect primarily manifests in PC2 of the infection subspace, all models, even sPCA, extract similar sets of top contributor genes in PC1 of the infection subspace, thus yielding GO enrichment results comparable to those shown in Fig. 6 . H Supplementary visualizations of baseline performance in applications.For completeness, here we provide subspace visualization for baseline models not included in the main figures. Specifically, Fig. 11 shows the full sPCA performance on the simulated data, related to Fig. 3 . Fig. 12 visualizes the single-cell subspaces learned by sPCA, VAE, supVAE, and hsVAE-sc, related to Fig. 5 .  
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?Answer: [Yes] Justification: We discuss the computational complexity of sisPCA in Appendix D. The proposed method is highly efficient because of the linearity. We ran all provided notebooks (https://github.com/JiayuSuPKU/sispca/tree/main/docs/source/tutorials) using a personal M1 Macbook Air with 16GB RAM and completed most analysis steps in minutes including model training.Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines?Answer: [Yes] Justification: We confirm the work conforms with the code of ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.• The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none which we feel must be specifically highlighted here. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: This paper poses no such risks. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [NA] Justification: The paper does not involve crowdsourcing nor research with human subjects. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.A Python implementation of sisPCA is available on GitHub at https://github.com/JiayuSuPKU/sispca (DOI 10.5281/zenodo.13932660). The repository also includes notebooks to reproduce results in this paper.uciml/breast-cancer-wisconsin-data, CC BY-NC-SA 4.0 license.Data access is controlled through the NCI Genomic Data Commons (https://portal.gdc.cancer.gov/).Preprocessed data available at https://doi.org/10.6084/m9.figshare.22148900.v1 (CC BY 4.0).https://github.com/JiayuSuPKU/sispcascvi-tools v1.2.0 (https://scvi-tools.org/) 
 
 Appendix:  


Decision: All four reviewers argued to accept the paper, albeit three of them only weakly. However, the paper has improved quite a bit relative to initial reviewer concerns during the discussion period, particularly noting the additional experiment which broadens the scope from the focus on single-cell RNA sequencing data; and the study into automatic settings of \lambda and the sensitivity to the parameter choice. Multiple reviewers raised their score after the author response.

Please be certain to incorporate everything discussed in the rebuttal phase when preparing the final version of the manuscript.

Reviews: review: paper_summary: paper_summary: This paper presents a method for distentangling multiple independent linear latent subspaces that align with a set of response variables. This uses the Hilbert-Schmidt Independence Criterion (HSIC) to measure the dependence each subspace and the targe variable providing supervision for that subspace, a concept which was previously applied to infer a single subspace that was connected to a set of response variables. HSIC is additionally used to enforce independence between the subspaces, to encourage the method to find independent subspaces that correlate with each of the provided dependent variables. This expands on Supervised PCA by enabling inference of an indendent subspace for each response variable, rather than a single subspace which aligns with the vector of response variables, reflecting the fact that different responses may be interacting with different subspaces in the data, e.g., methylation changes reflecting natural human aging and changes related to a disease like cancer.

The method is evaluated on synthetic data and shown to do a good job at extracting the different subspaces. In particular, Supervised PCA is dominated by one of the two supervised subspaces. 

The method is then evaluated on real methylation data and single-cell gene expression data. In the cancer data, the new method is shown to provide subspaces with greater separatio
strength_weakness: strength_weakness: One could argue that the work is somewhat incremental. The use of the HSIC is not new and nor is the idea of providing supervision to a latent variable model to ensure the latent subspaces reflected some additional variables. However, I feel the contribution is novel enough to warrant presentation at NeurIPS
questions: questions: No questions
limitations: limitations: The authors describe the well-known limitations of unsupervised PCA, i.e., identifiability, and other limitations implied by their choice of kernel.
paper_summary: paper_summary: The work proposes a new method that determines the factors of variation bonded to different labels in a supervised way, obtaining a method akin to supervised PCA but in several independent subspaces thorough the Hilbert-Schmidt independent criterion (that it’s employed to maximize both the independence between different subspaces and the dependence between the data subspaces and the labels). The paper focus on a linear kernel application of HSIC, but other kernels are also considered. Then it’s applied to some simulated data and to DNA/RNA sequencing data.
strength_weakness: strength_weakness: The novelty is limited, it appears to be an extension of previously described methods. 

The comparison is provided only with PCA or sPCA, which are not by any means SOTA methods.  

The explanation of the simulated data is not clear, what should it be the expected result? Are you learning the dimensions of the subspaces?  

A better description of the experiments is needed. How do you choose the hyperparameters (for instance, the dimension of the subspaces is 10 in both sections 4.2 and 4.3, why?). 

In section 4.2, the use of HSIC as a quality metric is not fair, since it’s already included in your loss.  

The significance of the paper would be increased if the authors prove that the method is useful for the analysis of other types of real-world datasets. 

The discussion is too short and not too clear. I acknowledge the limitations of space, but the discussion is one of the most important sections and should be properly done (even by moving some mathematical details to SI).
questions: questions: Could you please provide a comparison of your method with other SOTA methods? 

Would you prove that the method is useful in other kinds of data types? 

Which are the hyperparameters of the method? How do you set them?
limitations: limitations: I’m not satisfied with the way the limitations are addressed. One should explain them clearly and provide examples in which the method is not working.
paper_summary: paper_summary: This paper proposed a linear dimensionality reduction and subspace extraction method based on Hilbert-Schmidt Independence Criterion (HSIC) and Supervised PCA. In specific, several interpretable subspaces, which are independent from each other, are disentangled from the data observations and the leaned subspaces are enforced to be maximally dependent on the supervision variables. The proposed approach is shown to be effective on both synthetic data and two real datasets.
strength_weakness: strength_weakness: 1.	It is not clear if the proposed method is guaranteed to extract the latent subspaces and/or under what conditions/assumptions the method would work or fail. For example, in the synthetic experiment, the mixing matrix is uniformly drawn from [0, 1]. Does it still work if the matrix is changed to Gaussian? Overall, the theoretical analysis is somewhat lacking, e.g., identifiability analysis for linearly mixed subspaces.
2.	Baselines seem missing in the real data applications. Only PCA and supervised PCA are considered here. The authors are encouraged to include more state-of-the-art methods to showcase the superiority of the proposed method, e.g., how does PCA work under such a setting? Are other HSIC-based nonlinear methods, e.g., [1], capable of extracting the independent subspace with interpretability?

[1] Lopez, Romain, Jeffrey Regier, Michael I. Jordan, and Nir Yosef. "Information constraints on auto-encoding variational bayes." Advances in neural information processing systems 31 (2018).
questions: questions: 1.	In Sec. 3.2, the linear kernel is used to target variables. Does it make a difference if a nonlinear kernel function is used?
2.	What are A and B in Fig. 2? Are they the target variables? In the same figure, why does it maximize dependence (HSIC) with the components of a subspace? I thought the only two objectives are 1) maximizing dependence between the subspace and its corresponding target variable; 2) minimizing the dependence between the subspaces. How is this maximization reflected in Eq. (2)?
3.	What does balanced/unbalanced supervision mean mathematically in Conjecture 3.1?
limitations: limitations: The authors are encouraged to discuss limitation of the work in the main paper.
paper_summary: paper_summary: Here the authors propose an extension of PCA that decomposes datasets into multiple independent subspaces that are encouraged to reflect provided covariates of interest. To enforce these dependencies, the authors leverage the Hilbert-Schmidt Independence Criterion (HSIC).
strength_weakness: strength_weakness: While my views of the manuscript are generally positive, I do have a couple of minor concerns. If the authors are able to address my concerns I would be happy to raise my score.

* **Impact of $\lambda$**: The authors' method attempts to find a balance between encouraging dependence between subspaces and their corresponding covariates while also minimizing the dependence between individual subspaces, and this balance is tuned via $\lambda$ in Equation (1). I imagine this parameter could have a significant impact on the embeddings/loadings produced by sisPCA, though this isn't explored in the manuscript. It would be great if the authors could perform additional experiments exploring these changes. For example, do the GO results presented in Figure 5 vary for different values of $\lambda$? How different do the UMAPs look for the scRNA-seq data for different values of $\lambda$?

* **Automatic selection of $\lambda**: Related to my previous point, the authors briefly mention that $\lambda$ could be selected using a procedure similar to [1]. It would be great to see an example of this in practice: forcing the user to select $\lambda$ manually seems like it could present a big obstacle in terms of usability, so I think it would be useful to see how the automated method performs.

* **Re: connections with self-supervised learning**: In the abstract the authors note "We elucidate the mathematical connections of
sisPCA with self-supervised learning and regularized linear regression". While I was able to locate the regularized linear regression results in Appendix B.1, I couldn't find any additional mathematical details on the self-supervision connection (the only other mention of self-supervision in the main text was in line 83 in the related works section). Could the authors expand upon this connection? (Please let me know if I missed anything here).

Minor points (these did not affect my score):

* It appears that the authors may have used \citet accidentally instead of \citep throughout the manuscript.
* Figure quality was relatively low-resolution when I printed out the manuscript- perhaps the authors could try increasing the DPI?
* In line 91 the authors mention "contrastive learning"- it may be worth putting explicit citations to e.g. Abid et al. [1] to clarify for the reader that "contrastive learning" here isn't referring to contrastive in the sense of e.g. SimCLR.
* I couldn't find any details on how the real-world biological datasets were preprocessed (e.g. I assume library size normalization was applied to the scRNA-seq data?)

[1] Abid et al., "Contrastive principal component analysis"
questions: questions: Did the authors try applying their method with different values of $\lambda$ to the real-world datasets?
limitations: limitations: The greatest limitations I see are (1) the method relying on linear transformations (as opposed to e.g. deep neural networks) and (2) potential difficulties in setting $\lambda$. For (1) the authors specifically note that they are intentionally trading off expressivity for increased interpretability

Comments: title: Thank you
Thank you
Thank you
New Application on Breast Cancer Diagnostic Data
Follow up
comment: Thank you for the feedback and suggestion. We will add this application as a new subsection and adjust the main text accordingly.
I thank the authors for their reply. This new application makes both the application and the method much clearer, and I would put it in the main text (just before 4.2). I will raise my score accordingly.
Thank you for your positive feedback. We very much appreciate your recognitions and all your helpful suggestions for improving our paper. We will include the automatic selection procedure as part of the package release. To further enhance usability, we will also release the re-implemented HCV (hsVAE-sc) models for single-cell data and provide corresponding tutorials on example datsets and general documentations.
Thank you for your feedback. We appreciate your recognition and the suggestion on adding the non-linear baselines which has helped improve the quality of our work. We will include relevant benchmark results in the final version.
We appreciate your feedback and apologize for our previous insufficient response. Due to space constraints, we couldn't include additional figures and results initially. To further illustrate sisPCA's versatility as a plug-in replacement of PCA for disentangling quantities, we provide below a new application, which will be included in the appendix of the final paper, and as a tutorial upon package release.

### **Problem and Dataset**
We used the Kaggle Breast Cancer Wisconsin Data Set (569 samples, 30 summary features, uciml/breast-cancer-wisconsin-data, CC BY-NC-SA 4.0 license). The 30 real-valued features are computed from imaging data of breast mass, which include the mean, standard error and the extremum of quantities like cell nuclear radius, texture, perimeter etc.
 
Here, our goals are to:

1.	Learn a representation for predicting sample disease status (Malignant or Benign, not used during training).
2.	Understand how original features contribute to data variability, the learned representation, and diagnosis potential.

### **Experiments and Results**
We focus our comparison on three linear models: PCA, sPCA, and sisPCA, using zero-centered and variance-standardized data as inputs. The diagnosis label used to measure subspace quality remains invisible to all models. Below we summarize the quantitative results, and will include subspace visualization in the final paper.

#### **PCA**
We projected all features (dim = 30) onto one PCA subspace (dim = 6, determined by the elbow rule), explaining 61.6% of total variance. Malignant and benign samples appear well separated in PC1 and PC2. 'symmetry_mean' (loading = -0.223) and 'radius_mean' (loading = -0.219) are the top 2 features that negatively contribute to PC1. That is, the higher the two feature scores, the lower PC1 score and the greater the possibility of the sample being malignant. 

#### **sPCA**
From the PC1 loadings, we sought to construct two subspaces to separately reflect nuclei size (‘radius_mean’ ) and shape ('symmetry_mean'). We set $Y_{radius}$ ('radius_mean', 'radius_sd') of $569 \times 2$ as the target variable for the radius subspace, and $Y_{symmetry}$ ('symmetry_mean', 'symmetry_sd') of $569 \times 2$ as the target for the symmetry subspace. The remaining 26 features were projected onto the two subspaces using sPCA (dim = 3, effective dim = 2). *Both subspaces better explained diagnosis status than PCA (**Table S3**) but remained highly entangled.* Specifically, the PC2 loadings of the two spaces have a Pearson correlation of 0.716, and 'perimeter_worst', which is highly correlated to 'radius_mean' (corr = 0.965), also strongly contributes to the PC2 of the symmetry subspace (loading = 0.238). 

#### **sisPCA**
We next applied sisPCA ($\lambda$ = 10) to further disentangle the two subspaces, increasing their separation (Grassmann distance increased from 1.593 to 2.058, PC2 loadings correlation decreased to 0.257). Here ‘perimeter’ features no longer contribute to the symmetry subspace. *As a result, the radius subspace remained predictive of diagnosis, while the symmetry subspace became less relevant (**Table S3**).* We confirmed the finding by directly measuring the predictive potential of the target variables $Y_{radius}$ and $Y_{symmetry}$ (Silhouette score = 0.457 and 0.092, respectively). 

**Table S3: Predictability of diagnosis status, measured by Silhouette score**

| | Radius subspace (dim = 3) | Symmetry subspace (dim = 3) | Overall subspace (dim = 6) |
|:-|:-|:-|:-|
| PCA (one subspace) | 0.294 (PC 1-3) | 0.013 (PC 4-6) | 0.160 (PC 1-6) |
| sPCA | 0.534 | 0.374 | 0.464 |
| sisPCA | 0.553 | 0.027 | 0.511 |

### **Interpretation**
1.	Our results suggest that nuclear size (radius subspace) is more informative for breast cancer diagnosis than nuclear shape (symmetry subspace), aligning with clinical observations [1].
2.	Unsupervised PCA and sPCA captured both radius and symmetry aspects as they contribute most to data variability. *Without disentanglement, the two aspects remain intertwined, potentially leading to incorrect conclusions about symmetry's predictive power for diagnosis.*
3.	*sisPCA successfully separated the two representations, revealing their distinct relationships to diagnosis.* The results are interpretable: the radius subspace is constructed using features like 'area' and 'perimeter', while the symmetry subspace uses features like 'compactness' and 'smoothness'.

This example demonstrates sisPCA's ability to disentangle different aspects of data variation and uncover underlying relationships. **Importantly, sisPCA improves upon the diagnosis potential of PCA's representation, even through the diagnosis labels were never used during training.**

[1] Kashyap, Anamika, et al. "Role of nuclear morphometry in breast cancer and its correlation with cytomorphological grading of breast cancer: A study of 64 cases." Journal of cytology 35.1 (2018): 41-45.
Dear authors,

Thank you for your thoughtful responses! I've accordingly raised my score. 

Based on the rebuttal I have a few remaining minor comments:

* I would encourage the authors to avoid the name "sVAE" when describing the supervised VAE baseline, as this name has previously been used for a different class of models [1,2]
* Two recent works [3,4] have investigated using HSIC-based penalties with VAEs for similar supervised disentanglement tasks; for completeness it would be great if the authors could include citations to these works
* Please do include your final results re: the automatic selection procedure for $\lambda$, as I believe this will make a big impact in terms of usability of the method.

[1]: Lachapelle et al. "Disentanglement via mechanism sparsity regularization: A new principle for nonlinear ICA"
[2]: Lopez et al. "Learning Causal Representations of Single Cells via Sparse Mechanism Shift Modeling"
[3]: Tu et al. "A Supervised Contrastive Framework for Learning Disentangled Representations of Cell Perturbation Data"
[4]: Qiu et al. "Isolating structured salient variations in single-cell transcriptomic data with StrastiveVI"
Thank you for the clarification and for including the non-linear VAE-based baselines, which should make this work more solid. I agree that linear methods have more benefits in terms of interpretability and simplicity.
I thank the authors for their reply and acknowledge their usefulness. I will, accordingly, raise my score. 
However, I have to say that I'm not satisfied with the reply to "other kinds of data sets" reply and it refrains me from further raising the score.