$\sigma$-PCA: a building block for neural learning of identifiable linear transformations
Authors: Fahdi Kanavati
Lucy Katsnith
Masayuki Tsuneki
Venue: Accepted by TMLR ()

Abstract: Linear principal component analysis (PCA) learns (semi-)orthogonal transformations by orienting the axes to maximize variance. Consequently, it can only identify orthogonal axes whose variances are clearly distinct, but it cannot identify the subsets of axes whose variances are roughly equal. It cannot eliminate the subspace rotational indeterminacy: it fails to disentangle components with equal variances (eigenvalues), resulting, in each eigen subspace, in randomly rotated axes. In this paper, we propose  $\sigma$-PCA, a method that (1) formulates a unified model for linear and nonlinear PCA, the latter being a special case of linear independent component analysis (ICA), and (2) introduces a missing piece into nonlinear PCA that allows it to eliminate, from the canonical linear PCA solution, the subspace rotational indeterminacy — without whitening the inputs. Whitening, a preprocessing step which converts the inputs into unit-variance inputs, has generally been a prerequisite step for linear ICA methods, which  meant that conventional nonlinear PCA could not necessarily preserve the orthogonality of the overall transformation, could not directly reduce dimensionality, and could not intrinsically order by variances. We offer insights on the relationship between linear PCA, nonlinear PCA, and linear ICA — three methods with autoencoder formulations for learning special linear transformations from data, transformations that are (semi-)orthogonal for PCA, and arbitrary unit-variance for ICA. As part of our formulation, nonlinear PCA can be seen as a method that maximizes both variance and statistical independence, lying in the middle between linear PCA and linear ICA, serving as a building block for learning linear transformations that are identifiable.

Published in Transactions on Machine Learning Research (06/2024)
o-PCA: a building block for neural learning of identifiable linear transformations
Fahdi Kanavati Invertible AI, London, U.K.
fannovati@invertibleeccccm
Lucy Katsnith
Invertible AI, London, U.K.
Masayuki Tsuneki Invertible AI, London, U.K.
mttumekki@invertiblee..ccmm
Reviewed on OpenReview: https: //tps://oet opt.reveeee....e...mmm.m net/ forum? id= wwt.kppW6660 JJ60nw/
Abstract
Linear principal component analysis (PCA) learns (semi-)orthogonal transformations by orienting the axes maximize variance. Consequently, can only identify orthogonal axes whose variances are clearly distinct, but it cannot identify the subsets of axes whose variances are roughly equal. It cannot eliminate the subspace rotational indeterminacy: it fails to disentangle components with equal variances (eigenvalues), resulting, in each eigen subspace, in randomly rotated axes. In this paper, we propose o-PCA, method that (1) formulates a unified model for linear and nonlinear PCA, the latter being a special case of linear independent component analysis (ICA), and (2) introduces a missing piece into nonlinear PCA that allows climinate, from the canonical lincar PCA solution, the subspace rotational indeterminacy without whitening the inputs. Whitening, preprocessing step which converts the inputs into unit- inputs, has generally been prerequisite step for linear ICA methods, which meant that conventional nonlinear PCA could not necessarily preserve the orthogonality of the overall transformation, could not directly reduce dimensionality, and could not intrinsically order by variances. We offer insights on the relationship between linear PCA, nonlinear PCA, and linear ICA three methods with autoencoder formulations for learning special lincar transformations from data, transformation that (stmi-)orthogonal for PCA, and arbitrary unil-variance for ICA. As part of our formulation, nonlinear PCA can be seen a method that maximizes both variance and statistical independence, lying in the middle between linear PCA and linear ICA, serving as a building block for learning linear transformations that are identifiable.
1 Introduction
Principal component analysis (PCA) (Pearson 1901 Hotelling 1933) needs no introduction. is classical, ubiquitous, perennial. is an unsupervisod learning method that can be arrived at from three paradigms of representation learning [Bengio et al. 2013): neural, manifold, and probabilistic. This paper is about learning linear transformations from data using the first paradigm: the neural, in the form of single-layer autoencoder model which encodes the data into new representation and then decodes that represcntation back into reconstruction of the original data.
From the data, PCA learns a linear point x into a new reprosentation y xW, representation with reduced dimensionality and minimal loss of information (Diamantaras & Kung 1996 Jolliffe Cadima 2016). This representation has components that uncorrelated, i.e. have linear associations, or, in other words, the covariance wetween two distinct.
transformation W that transforms a given data

---

Published in Transactions on Machine Learning Research (06/2024)
components is zero. When the variances of the components are clearly distinct, PCA has unique solution, up to sign indeterminacies, that consists in a set of principal eigenvectors (or axes) ordered by their variances. This solution can usually be obtained from the singular value decomposition (SVD) of the data matrix (or the cigendecomposition of the data covariance matrix). Such a solution not only maximizes the amount of variance but also minimizes the error of reconstruction. And it is the latter in the form of the squared difference between the data and its reconstruction that is also the loss function of a single-layer lincar autoencoder. However, minimizing such a loss will result in a solution that lies in the PCA subspace (Baldi & Hornik, 1989; Bourlard & Kamp 1988 Oja 1992b), i.e. a linear combination of the principal eigenvectors, rather than the actual principal eigenvectors. And so, linear autoencoder performs PCA but not in its expected axis-aligned form.
What happens when we insert a nonlinearity in the single-layer autoencoder? By natural extension, it must 1995 Karhunen et al. 1998). However, to yield an interesting output, was observed that the input data needs to be whitened, i.c. rotated and
be a form of nonlinear PCA (Xu 1993 Karhunen & Joutsensalo, 1994
scaled into uncorrelated components of unit variance. And after the input is whitened, applying nonlinear
PCA results in similar output (Karhunen & Joutsensalo, 1994 1995
Karhunen et al 1998 Hyvärinen &
Oja 2000
to lincar independent component analysis (ICA) (Jutten & Herault
1991
Comon
1994
Bell &
Sejnowski 1997 Hyvarinen 1999..
Linear ICA, too with a long history in unsupervised learning [Hväärinnn & Oja 2000 yields, as the name implies, not just uncorrelated but independent components. Linear ICA seeks to find a linear transformation B not necessarily orthogonal such that y = xB has its components as statistically independent as possible. It is based the idea that x is an observed linear mixing of the hidden sources y, and it aims recover these sources solely based on the assumptions that they are non-Gaussian and independent. Generally, unlike PCA, we cannot determine the variances of the components nor their order (Hyvarinen 1999 Hyvärinen & Oja 2000 Hyvärinen 2015); they are impossible to recover when we have no prior knowledge about how the data was generated. For B is a transformation that recovers the sources, then so is BA, with A an arbilrary diagonal scaling matrix. And so, the best that can be done, without any priors, is to assume the sources have unit variance. Linear ICA guarantees identifiability when its assumptions are met (Comon) 1994., meaning that it can recover the true sources up to trivial sign, scale, and permutation indeterminacics.
We thus have three methods for learning linear transformations from data: linear PCA can learn rotations that orient axes into directions that maximize variance, but, by solely relying on variance, it cannot eliminate the subspace rotational indeterminacy: it cannot separate components that have the same variance, i.e. it cannot identify the true principal eigenvectors in a given eigenspace, resulting in randomly rotated eigenvectors the eigenspaces. Nonlinear PCA (Karhunen & Joutsensalo 1994) can learn rotations that orient axes into directions that maximize statistical independence, reducing the subspace indeterminacy from rotational to trivial permutational, but it cannot be applied directly to data without preprocessing whitening. Linear ICA learns linear transformation not just rotation that points the axes into directions that maximize statistical independence under the assumption of unit variance. The relationship between all three can be best understood by the SVD of the linear ICA transformation into a sequence of rotation, scale, rotation - stated formally, WE 'V, with W (semi-)orthogomme diagonal. In one formulation of linear ICA, linear PCA learns the first rotation, W, and nonlinear PCA learns the second, V. The scale, Σ-!, simply the inverse of the standard deviations.
V orthogonal, and Σ
The problem is that, in contrast to linear PCA, conventional nonlinear PCA cannot be used directly on the data to learn the first rotation, the first being special it can reduce dimensionality and order by variances. The process of first applying the linear PCA rotation and then dividing by the standard deviations is the whitening preprocessing step, the prerequisite for nonlinear PCA. And so, nonlinear PCA, rather than being on an equal footing, is dependent on linear PCA. Conventional nonlinear PCA, by the mere introduction of a nonlinear function, loses the ability for dimensionality reduction and ordering by variances a jarring disparity between both linear and nonlinear PCA.
In the literature, for ICA, the terms linear and nonlinear refer the transformation, while, for PCA, they refer to the function. Both linear ICA and nonlinear PCA use nonlinear functions, but their resulting transformations are still linear.
2

---

Published in Transactions on Machine Learning Research (06/2024)
As part of our contributions in this paper, we have identified the reason why conventional nonlinear PCA has been unable to recover the first rotation, W: it has been missing, in the reconstruction loss, Σ. This means that the nonlinear PCA model should not just be W but WE-1, The reason why Σ is needed is simply because standardizes the components to unit variance before applying the nonlincarity - while still allowing us to compule the variances. consequence of introducing Σ that nonlinear PCA can now learn not just the second, but also the first rotation, i... can also be applied directly to the data without whitening. Another key observation for nonlinear PCA to work for dimensionality reduction that it should put an emphasis not the decoder contribution, but on the encoder contribution in contrast to conventional linear and nonlinear PCA. In light of the introduction of Σ, and to distinguish our proposed formulation from the conventional PCA model, the conventional being simply a special case with Σ = I, we call our formulation o-PCA. E
Thus, our primary contribution a reformulated nonlinear PCA mothod, which naturally derives from our o-PCA formulation, that can climinate, from the canonical linear PCA solution, the subspace rotational indeterminacy without whitening the inputs. And so, like linear PCA, nonlinear PCA can learn a semi- orthogonal transformation that reduces dimensionality and orders by variances, but, unlike linear PCA, it can also climinate the subspace rotational indeterminacy - it can disentangle components with equal variances. With its ability to learn both the first and second rotations of the SVD of the linear ICA transformation, nonlinear PCA serves as a building block for learning linear transformations that are identifiable.
2
Preliminaries
2.1
Learning linear transformations from data and the notion of identifiability
Given a set of observations, the goal is to learn a linear transformation that uncovers the true independent sources that generated the observations. It is as the observations are the result of the sources becoming linearly entangled (Bengio et al. 2013) and what we want to do is disentangle the observations and identify the sources (Comon 1994 Hyvärinen Oja 2000).
Let Y. R be the data matrix of the true independent. sources, with y. Rlik row vector representing source sample, and X E gnxp be the data matrix of the observations, with x E gxxp row vector representing an observation sample. As the observations can have a higher dimension than the sources, we have p 2 k. Without loss of generality, we assume the sources are centred. The observations were generated by transforming the sources with a true transformation B B-] E Rxxp, the left inverse of B., to obtain the
observations Y.B-1 (L))
To learn from data is to ask the following: without knowing what the true transformation was, to what extent can we recover the sources Y. from the observations X? Or other words, to what extent can we recover the left invertible transformation B.? The problem then to identify B such that
=B
(2.1)
with y as close as possible to y.. In the ideal case, we want B = B., allowing us to recover the true source as = xB. = y.B (L) B. = y.. We can note that
xB.B =.B (L) (L)
B.B (L) = y.B =x
(2.2)
implying that B, salisfies an autoencoder formulation, with B. representing the encoder and B reprosenting the docoder. This means that, in theory, the minimum of the reconstruction loss
(D)
2 Arguably, e-PCA should be the general PCA model, for it is common to divide by the standard deviations when performing the PCA transformation. Therefore, instead of using the names linear o-PCA and nonlinear o-PCA, shall simply refer to the model with Σ I as the conventional PCA model.
TTo disentangle mcans recover the sources but not necessarily their distributions, while identify means to recover not just. the sources but also their distributions. Identifiability implies disentanglement, but disentanglement does imply identifiabilitt.
3

---

Published in Transactions on Machine Learning Research (06/2024)
E[|x - xBB xBB (2)118)
(2.3)
can be reached when B = B., giving us a potential way to recover B.. Alas, B. is not the only minimizer of the reconstruction loss. Equally minimizing the loss is B,D, where E Rkxk is an invertible matrix with orthogonal columns (Baldi & Hornik 1989). This is because we can write BB (L) BDD BB))
(DB(( The matrix D is indeterminacy matrix, and it can take the form of a combination of certain special indeterminacies:
Scale
A diag( (.,.... As]) E R4*
(2.4)
Sign
I
=1]) O(k),
(2.5)
Permutation
P E O(k) permutation matrix,
(2.6)
Rotation
R E (O(k) an orthogonal matrix with det(R) = 1.
(2.7)
Any orthogonal matrix Q can be expressed as an arbitrary product of other orthogonal matrices, and, more specifically, as the product of rotation (determinant +1) and reflection (determinant -1) Both the sign and the permutation matrices can have determinant +1 -1. Consequently, can express orthogonal indeterminacy matrix Q as I but it can equally be an arbitrary product of I,, P, and R, e.g.
2= = =R" R''
(2.8)
An indeterminacy can be over the entire space or only a subspace, i.e. affecting only subset of the columns. If it affects subspace, then shall denote it with the subscript e.g. subspace rotational indeterminacy, which takes the form of block orthogonal matrix, will be denoted as R.,
The worst of the four indeterminacies the rotational indeterminacy. The other three are trivial in comparison. This is because a rotational indeterminacy does a linear combination of the sources, i... it entangles the sources, making them impossible identify is the primary reason why representations are not reproducible.
A linear transformation said to be identifiable if it contains no rotational indeterminacy only sign, scale, and permutation indeterminacies. In other words, transformation is identifiable if for any two transformations A and B satisfying y xA and xB, then A = BAPI+
2.2 Singular value decomposition (SVD)
Let B E ppxk represent linear transformation, then its reduced SVD
Trefethen & Bau III 1997 is
B USVT,
(2.9)
where U Baxk semi-orthogonal matrix with each column left eigenvector, V E O(k) orthogonal matrix with each column right eigenvector, and S B diagonal matrix of ordered eigenvalues.
Sign indeterminacy If all the eigenvectors have distinct eigenvalues, then this decomposition is, up to sign indeterminacies, unique. This is because the diagonal sign matrix can commute with S, so we can write
B = USV = UI = "
(2.10)
O(1). This implies that an orthogonal matrix can be mapped rotation trivial reflection, such as diagonal matrix of all except in the first row.
4This is because O(k) t SO(K)
4

---

Published in Transactions on Machine Learning Research (06/2024)
Subspace rotational indeterminacy If multiple cigenvectors share the same cigenvalue, then, in addition to the sign indeterminacy, there is a subspace rotational indeterminacy. The eigenvectors in each eigen subspace are randomly rotated or reflected. This means that any decomposition with an arbitrary rotation or reflection of those cigenvoctors would remain valid. To see why, let Q, = R, be a block orthogonal matrix where each block allects the subset of eigenvect.rs that share the same eigenvalue, then, given that it. can commute with S, i.c. SQ, = Q₂S (sce Appendix we can write
USV = USV'=UIA ,. = USS"
(2.11)
2.3 Linear PCA via SVD
Linear PCA can be obtained by the SVD of the data matrix. Let X E Roxp be the data matrix, where each row vector x is data sample, then we can decompose X into
X USV SU,
(2.12)
Letting Y= UR,L, and W = LR,L, and E = Ja 8, we can write
X YEW?,
(2.13)
= XWE-1
(2.14)
The columns of W the principal eigenvectors, the columns of Y the principal components, and E is ordered diagonal matrix of the standard deviations.
Because of the presence of R,, linear PCA has subspace rotational indeterminacy, which means that it cannot separate, cannot disentangle, components that share the same variance.
2.4 Degrees of identifiabiliy via variance and independence
From the SVD, we know that any linear transformation can be decomposed into a sequence of rotation, scale, rotation. A rotation - an orthogonal transformation is thus a building block. The first rotation, in particular, can be used for reducing dimensionality, so can for finding rotations from data: maximizing variance and maximizing independence. The former fails for components with equal variance, and the latter, although succeeding for components with equal variance, fails for components with Gaussian distributions. The latter fails because of two things: muit-variance Gaussian has spherical symmetry and maximizing statistical independence is akin to maximizing non- non-Gaussianins
ortoooalll.o.a.. There are two ways
Con-Gaussianity is an essential assumption of linear ICA (Comon 1994 to guarantee identifiability. With a linear transformation, identifiability theory guarantees that independent sources with unit variance can be recovered least all but one of the components have non- non-Gaussian distributions (Comon 1994 Hyvärinen & Oja 2000). Figure illustrates the identifiability rotations using maximization of variance or independence. Although maximizing independence cannot. identify sources with Gaussian distributions, maximizing variance still identify the sources for one particular case: when the true transformation is orthogonal and the variances are all distinct (the ellipse in Fig. D see Appendix D.3
Let us now consider how and to what degree we can identify the linear transformation based on whether it is (semi-)orthogona, or non-orthogonal. Although we will be expressing what is idontified in terms the truc transformations, in reality, what we end up recovering, due to noise and absence of infinite data, are but approximations of the true transformations.
iI theory, this is the case and the rotation is random, but practice, the eigenvalues exactly equal SVD solver could still numerically reproduce the fixed rotational indeterminacy (Hyvärinen 2009).
5

---

Published in Transactions on Machine Learning Research (06/2024)
01 7 d2
01 02
01
Entangled x
Identified y
Variance
Independence
Figure 1: Identifiability of rotations using variance independence maximization. The aim is to find rotation that transforms x into aligning the axes. Red indicates unaligned; green, aligned. For the rectangle and the ellipse, representing distinet variances, we can figure out rotation, except that we do not know if we should flip up or down, left or right a sign indeterminacy. For the square and the circle, representing equal variance, we cannot figure out a rotation only from the variance. For the square, however, we can use independence figure out a rotation, except that we do know if we should flip and/or permute the vertical and horizontal axes sign and permutational indeterminacies. For the circle, nothing can be done. There are no favoured directions; even if we apply rotation, it would remain the same rotational indeterminacy. A Gaussian is like a circle.
2.4.1 Semi-orthogonal transformation
Suppose that B, semi-orthogonal, i.e. BiZ) BTT with BT I. Let E, denote the matrix of independent orthonormal eigenvectors ordered by their variances, so that true transformation B, E..
The primary method for learning semi-orthogonal transformations is linear PCA. The extent of how much it can identify the true transformation depends the variances of the components:
E.L+
all the variances are distinct identifiable;
B
E.I.R,,
if subset of the variances are equal partially identifiable;
(2.15)
E. R,
all the variances are equal unidentifiable.
If B. had orthogonal but non-orthonormal columns, then there is also a scale indeterminacy, arising because we can introduce diagonal matrix A such that B B" A AB B
B B'.
Although conventional nonlinear PCA can resolve the rotational indeterminacy (Karhunen & Joutsensalo 1994 1995; Karhunen et al.) 1998). currently can only work after linearly transforming the inputs into decorrelated components with unit variance after preprocessing step called whitening. This can be better understood when we consider the non-orthogonal transformation.
2.4.2 Non-orthogonal transformation
Suppose that B, non-orthogonal. With a non-orthogonal matrix, the columns do not necessarily have unit norm, and so, without any prior information, there is a scale indeterminacy, arising because we can introduce diagonal matrix such that B B = Let us suppose that the true transformation B. U A., where A, represents the standard deviations of the independent components, and U, represents the linear transformation that results in unit variance components. The scale A, can only ever be recovered if prior knowledge is available about the norm of the columns of B, B-1 BT but that is rarely the case.
= B there is no way tell what A is.
Linear ICA allows us to recover U, up to sign and permutation indeterminacies. It exploits the fact that once the data has been whitened rotated and scaled into uncorrelated components of unit variance then any

---

Published in Transactions on Machine Learning Research (06/2024)
other subsequent rotation will maintain unit variance. What this means is that if A E R is whitening matrix, then [[(xA)xxA) = I, and if we attempt to apply another rotation V E O(k), then
B(V
XX)( V
(2.16)
Therefore, can be recovered in two stages by first finding whitening transformation A and then finding rotation V that maximizes independence, resulting in U AV. We already have method for finding whitening matrix: linear PCA, resulting in A E, + R. What is left now is to find V. For that, one method is the fast. fixed-point algorithm of FastTCA (Hyvarinen 1999); alternative is conventional
nonlinear PCA (Karhunen & Joutsensalo 1994). If U, E. V., then
V= V-RLVVL
(2.17)
resulting
U AV = E.L. -V
(2.18)
2.5 Conventional neural linear PCA
A single-layer tied linear autoencoder has a solution that lies in the PCA subspace (Baldi & Hornik 1989
Bourlard Kamp 1988 Oja 1992b). Lot X e R2xx be the centred data matrix, with x row vector X,
and W E RRXX, with k p, then we write the autoencoder reconstruction loss
((W) B(((W)) = E[[[x - xWW7 11).
(2.19)
From this we see that the autoencoder transforms the data with W (the encoder) into set of components = xW, and then it reconstructs the original data with W (the decoder) to obtain the reconstruction & = yWT. & < then it performs dimensionality reduction. The loss naturally enforces WTW and when it is enforced, the contribution to the gradient originating from the encoder is zero (see Appendix A.2.1 Given that, the linear PCA loss can be equally expressed in form that results in the subspace learning algorithm Oja 1989 weight update (see Appendix A.2.2
((W) = [[[[ x W 12)
(2.20)
where lsg is the stop gradient operator, an operator that marks the value of the weight as a constant i.e. the weight not treated as a variable that contributes to the gradient.
Let. E, R be the matrix whose columns ordered orthonormal independent eigenvectors, then the minimum of Eq. 2.19 2.20 is obtained by any orthonormal basis of the PCA subspace spanned by the columns of E.. That is, the optimal W has the form
E.R.R.R.
(2.21)
where R O(k) is an indeterminacy over the entire space, resulting in a combined indeterminacy R,R, which does not affect the loss, given that E.I,R, RR RRT R.R.RRRRTR
E,E.. R I, then the lincar PCA solution E. can be recovered. The indeterminacy R appears because the loss is symmetric: no component is favoured over the other. To eliminate R, simply need to break the symmetry (see Appendix A.3
7

---

Published in Transactions on Machine Learning Research (06/2024)
2.6 Conventional neural nonlinear PCA
A straight.orward extension [Xu 1993 linear to nonlinear PCA is to simply introduce in Eq. 2.19 a nonlinearity h, such as tanh to obtain
Karhunen & Joutsensalo, 1994 1995
Karhunen et al. 1998) from
((W) E[[[x h(xW)
(2.22)
This loss tends to work only when the input is already whitened. This means that, when the input has unit variance, this loss can recover a square orthogonal matrix that maximizes statistical independence. It can be shown that it maximizes statistical independence as exact relationships can be established between the nonlinear PCA loss and several ICA methods (Karhunen et al. 1998). In particular, these relat ionships become more apparent when noting that, under the constraint W W = I, Eq. 2.22 can also written as
L(W) E(Ily h(y)
(2.23)
which is similar to the Bussgang criterion used in blind source separation (Lambert
1996 Haykin 1996
Unlike linear PCA, conventional nonlinear PCA cannot recover useful semi-orthogonal matrix that reduces dimensionality. With the current formulation, the mere introduction of nonlinear function creates a jarring disparity between both linear and nonlinear PCA. Ideally, we would like to be able to apply nonlinear PCA directly on the input to reduce dimensionality while maximizing both variance and statistical independence, allowing us to order by variances and to eliminate the subspace rotational indeterminacy.
3 の- pCPA:: unified neural model for linear and nonlinear PCA
Instead of generalizing directly from the linear reconstruction loss from Eq. 2.19 let take step back and consider the loss with B e Rexk, & a left invertible transformation:
B[Ix - h(BBBB (2)1) 1B3.
(L)
(3.1)
Conventional nonlinear PCA sets B = W, with W orthogonal. But if we consider that (1) the conventional nonlincar PCA loss requires unit variance input, and (2) the lincar PCA transformation, whether from the SVD decomposition (Eq. 2.14) or from the whitening transformation, takes the form WΣ-1 then it seems natural that B should not be just W but WE-1, with Σ a diagonal matrix of the standard deviations, The transformation WE-1 standardizes the components to unit variance before the nonlincarity h. We thus propose general form of the unified linear and nonlinear PCA model: a-PCA. Its reconstruction loss is
L(W, Σ) = E([[x MXWE-")EW" (2)
(3.2)
with WTW = I. Formulating it this way allows us to unify lincar and nonlincar PCA. Indoed, if is linear, then we recover the linear PCA loss: E[[xx- x[[xxxxW--1 ------------ Ew- = E(Ix - xWWF l1)..
From this general form of the loss, we will see that we can derive specific losses for linear PCA, nonlinear and 3.21 respectively.
PCA, and linear ICA in equations 3.19 3.20
3.1 A closer look at the gradient
To derive the new nonlinear PCA loss from the general form in Eq. 3.2 we need to make one modification: need to omit, from the weight update, the decoder contribution. To understand why, need to have
8

---

Published in Transactions on Machine Learning Research (06/2024)
closer look at the individual contributions of both the encoder and the decoder to the update of W. Given that gradient contributions are additive, we will be computing each contribution separately then adding them up. This is as if we untie the weights of the encoder (We, Σe) and the decoder (Wa, Σa), compute their respective gradients, and then tie them back up, summing the contributions. Without loss generality, we will look at the stochastic gradient descent updale (i.e. the loss & of a single sample), and multiply by 2 for convenience, so as to compute the gradients of
1 211x MXW
(We. Wa Le, En) =
(3.3)
Let xW e, z= y x= h(z) W and y = XW. By computing the gradients, we get
0l x W
(3.4)
OW. Oi
OW = (& - x))
(3.5)
Now that we have traced the origin of cach contribution, we can casily compute the gradient when the weights are lied, i.e. W- = .W = Wa and E = E. Eg, but we can still inspect each contribution separalely:
de
x x y) (z)z)
(3.6)
OW.
Ol
OW (& x)² {()(
(3.7)
Oe Ol de
(3.8)
dw oW. aw.
at From this, we can see the following: the contribution of the decoder, aW, will be zero when & ~ x, and that of the encoder, oW. when ý ~ y. Thus, unsurprisingly, the decoder puts the emphasis on the input reconstruction, while the encoder puts the emphasis on the latent reconstruction.
In the linear case, if WTW I, the contribution of the encoder is zero. This is because we have B4 x" w III 0. In fact, apart from acting as an additional orthonornalization term, its removal TW.
does not change the solution. We can remove it from the gradient by using the stop gradient operator leg:
(W)) = -|x x(Wlag WT
(3.9)
This results in
de
= =-xTy Wy y,
(3.10)
ow
which simply the subspace learning algorithm Oja 1989),
obtained from combining variance maximization
with an orthonormalliztion constraint (see Appendix A.2.2
In the conventional extension from linear to nonlinear PCA, the contribution of the encoder also removed KKarhnen & Joutsensalo 1994), as it was deemed appropriate to simply generalize the subspace learning algorithm (Ojal 1989 Karhunen & Joutsensalo 1994). But we find that it should be the other way around for the nonlinear case: what is essential for unwhitened inputs is not the contribution of the decoder, but the contribution of the encoder.
at Indeed, in the nonlinear case, aw, no longer zero. The problem is that it is overpowered by the decoder contribution, in particular when performing dimensionality reduction (see Appendix F.6). Arguably it would

---

Published in Transactions on Machine Learning Research (06/2024)
be desirable that the latent of the reconstruction is close possible to the latent of the input, i.c. put more emphasis on latent reconstruction. This means that for the nonlinear case we want to remove aw, ae which is opposite of the linear case. We thus arrive at what the loss
or
from the update and only keep aw.
should be the nonlinear case:
ANLPCA (BLPCA(W:E) = B( h(WW I [>[W"] (
subject to [WW.j]]l = 1,
(3.11)
where W[:j]] is the jth column of W. The corresponding gradient is
@CNLPCA ow
(3.12)
Although this loss naturally enforces the orthogonality the columns of W by the presence of the multiplicative term WTW, does not maintain their unit norm bocause W and can both be scaled by the same amount and WE-1 would remain the same, justifying the need for a projective unit norm constraint, W|: < [[ [W]llll
W
Though Σ could be trainable, without any constraints, it could increase without limit. For instance, if = tanh, we have lim tanh(z/o) = 2. One option is to add an L2 regularizer, but Σ need not be trainable, as we already have natural choice for the estimated standard deviation y. Indeed, we can set Σ [diag( / var(y))) (see Appendix F We discuss what needs to be in Section 3.4
3.2 Nonlinear PCA: an ICA method that maximizes both variance and independence
In the gradient update of nonlinear PCA (Eq. 3.12). ; we can see two terms, one corresponding to maximizing variance and the other to maximizing independence both under the constraint W I. These can be clearly seen in the stochastic weight update (omitting ! as is undone by unit normalization):
variance
AW
x - WWWW W
(3.13)
independence orthogonality
The variance maximizing term, xFy, can be obtained in isolation by minimizing the loss -1yy|13 which is the negative of the variance if taken in expectation. The independence maximizing term, -x"hhyE "WWW, can be seen as a generalization of the linear ICA update to non-unit variance input; indeed, assuming a square W E O(k) and unit variance, the variance term becomes constant in expectation, and the independence term reduces to xxy)(y). The linear ICA update of Hyvärinen & Ojal (1998) for pre-whitened inputs takes the form AW x +x" fyy), with the sign mostly dependent on whether the distribution sub- or super-Gaussian, determined by the sign of B(yf(y) f'()). If set f(v) y a tanhhy/a), then we have x'y-x"'(yy) x²y(()) = x" f(y,, thus recovering the exact same update, with the setting of the scalar similarly dependent on whether the distribution is sub- super-Gaussian (see Section 3.4p. Both variance and independence terms can be seen as Hebbian update rules, the former being linear, the latter nonlinear.
As linear ICA is based on maximizing non-Gaussianity, a general formulation Comon Oja 1998) of the linear ICA loss is
1994 Hyvärinen &
J(W) = (C(F(xW))EE(PP)) -
(3. 14)
under the constraint W W = I, where n ~ 10) and F is a smooth nonlinear function. This, however, assumes that the inputs are whiteneed For nonlinear PCA, we can make a generalization of this loss that takes the variances into account:
Tx (x(W) (((((WE-)))) -E E(F(n)
(3.15)
10

---

Published in Transactions on Machine Learning Research (06/2024)
And so, instead of comparing with Gaussian with unit variance, are now comparing with Gaussian with different variances. As the term involving the Gaussian is constant (Hyvärinen & Oja 1998), the linear ICA problem reduces to finding all the extrema points of the first term [[F(xW> ")E). Now, in similar way that conventional nonlinear PCA Karhunen et al., 1998) can be put in relation to lincar ICA, we can put nonlinear PCA in relation to linear ICA by noting that if W W-W- PCA loss can be re-expressed as
enforced, then the nonlinear
B(Ilhyy> [yleo 12).
(3.16)
If we chose F(z)= - (h(=) (z) then the nonlinear PCA loss becomes equivalent to the loss in Eq. 3.15 Given that, it can be shown that the loss will lead the columns of W to converge to the independent eigenvectors if the independent components 2, satisfy ( is satisfied, then the independent cigenvectors are stable stationary points of the loss. This guarantees local convergence to the extrema points the proof of this is given in Hyvärinen Oja (1998). Global convergence, however, cannot be easily shown, but numerical simulations tend to show that the linear ICA
# O, because, if the condition
Joss tends to converge globally (Hyvärinen & Oja 1998 Karhunen et al. 1998
Nonlinear PCA can be put in relation to other methods related to linear ICA, such as blind deconvolution and Ia sparsity constraints - in particular a linear autoencoder with an I regularization term. Both connections are discussed in Appendix F.3 and F.4
3.3 Relationship between linear PCA, nonlinear PCA, and linear ICA
Linear PCA, nonlinear PCA, and linear ICA seek to find matrix B E Rpxk such that xB. Lincar PCA puts an emphasis on maximizing variance, lin- ear ICA on maximizing statistical independence, and nonlincar PCA on both. All three assume that the resulting components of y have unit variance. This B is found by being decomposed into WE-'V, with W Iap* semi-orthogonal, Σ pixk diagonal, and V E O(k) orthogonal.
W
V
1, # 02
E
WL
Orthogonality a defining aspect of PCA, so in the case both linear and nonlinear PCA, we have V I, i.c. the restriction is for B have orthogo- nal columns; there no such restriction with ICA. Simply stated,
1,01
02
W
E
V
PCA
(Non)lincar PCA Linear ICA WV-IV.
WE-1
(3.17)
Figure 2: Linear PCA, nonlinear PCA, and linear ICA.
(3.18)
Green indicates alignment of axes.
As part of the linear ICA transformation, linear PCA can only find the first rotation, W, while nonlincar PCA can find both rotations, W and V. Figure O summarizes the differences between three methods based on the orthogonality of the overall transformation and the variances of the independent components.
The reconstruction losses of linear PCA (LPCA), nonlincar PCA (NLPCA), and linear ICA (LICA) follow the same pattern of [[[]| x MxxB) (Z) 113, but with few subtle differences, which we write
CLPCA(W) B(Ix x xWW
E([xx x W w1)
(3.19)
CNLPCA (W; Σ) = -B(|x <xW> 1x Jop 122,
(3.20)
LLICA (W, V; Σ) = -B([x[WE> '199 WX WE W "Jeg V) + CNLPCA (W;
(3.21)
11

---

Published in Transactions on Machine Learning Research (06/2024)
The value of Σ is set as [diag(o))ssg where o² var(xW) (see Section 3.5 particular, requires unit. norm constraint.
The nonlinear PCA loss, in
Linear PCA emphasizes input reconstruction, while nonlinear PCA emphasizes latent reconstruction. We can see this clearly from their losses and their respective weight updates:
Linear PCA
B([[x x [W] og 12),
(3.22)
Nonlinear PCA
B([[x h(xwe-
(3.23)
Linear PCA AW x
x - x)
(3.24)
Nonlinear PCA AW x (x"y y x'& o 1'y-
x (y 8) h'
(3.25)
Unlike the nonlinear PCA and linear ICA losses, the lincar PCA loss in Eq. 3.19 does not allow to obtain solution without rotational indeterminacy, in addition to the subspace rolational indeterminacy (see Section 2.5 but there are multiple ways to obtain one by breaking the symmetry. We list many existing methods in Appendix A.3
Omitting N' from the stochastic update, as it simply scaling that is undone by unit normalization, in the weight update of linear and nonlinear PCA, we can see the variance, independence, and orthogonality terms:
variance
orthogonality
Linear PCA
AW O x
Wy y
(3.26)
Nonlinear PCA
AW
x M'E
(3.27)
variance
independence orthogonality
Suppose that E, is the matrix whose columns are the true ordered orthonormal eigenvectors, and let WL be the axis-aligned linear PCA solution, WN the nonlinear PCA solution, then we have
E.I. E.I R.y E.I.R,
if all the variances are distinct identifiable;
Linear PCA
WL
if a subset of the variances are equal - partially identifiable; if all the variances are equal umidentifiable.
(3.28)
E, I
if all the variances are distinct identifiale; if subset of the variances equal identifiable;
Nonlinear PCA
WN
E.I.P... s.
(3.29)
E.I.P.
if all the variances are equal identifiable.
Now suppose that the true transformation is non-orthogonal of the form E.Σ Linear ICA can recover it two ways, each consisting two steps, akin to a two-layer model: (1) linear PCA followed by nonlinear PCA (with unit variance), and (2) nonlincar PCA followed by nonlinear PCA (with unit variance):
(1) WLE-IV E, (2) WNE-IV E.I. E.I. VII.
Linear ICA
(3.30)
The nonlinear PCA solution, Wx. exists within the linear ICA solution space. This because if we consider the transformation B WE-IV, the nonlinear PCA solution space is restricted to block orthogonal matrices
12

---

Published in Transactions on Machine Learning Research (06/2024)
V that commute with E i.c. E-1, ,
VV; whereas the linear ICA solution space includes not just the matrices V that commute, but also the ones that do not, i.e. V + VE-I, When V does not. commute the transformation B does not necessarily have orthogonal columns. Therefore, in theory but not necessarily in practice as the linear ICA model is over-parametrized with a larger solution space the nonlinear PCA solution can be recovered using linear ICA (as linear PCA followed by conventional nonlinear PCA). If V. I, then when we apply linear ICA, we obtain
(3.31)
To reduce the whole space permutation indeterminacy, we can simply renormalize E,E-I to unit norm columns and reorder by the estimated variances to obtain WN.
3.4 Choice of nonlinearity for sub- and Gaussian input distributions
0.8
2 "yr
Sub-Ganssan
0.6
0.
2
10.2
- +-
(a) Probability densily functions of three distributions
(b) a (anh(z/a))
Figure Three distributions with unit variance (a): uniform distribution (sub- Gaussian), Gaussian dis- tribution, and Taplace distribution (super-Gaussian) Shaded grey any |ar| 2 2a. When using tanh without any scale adjustment, we can see that for the super-Gausse Gaussian distribution (or for any hcavy-tailed distribution), values beyond 2a might have their reconstruction impaired because of the squashing by tanh (b). A remedy, in this case, would be to use tanh(z/a)) with a 2 1. For sub- nab-Gaussian distribution, it is more suited to use S as the values are within 20.
The choice of nonlinearity generally depends on whether the distribution is sub- or super-Gauses (Hyvärinen) & Oja 2000 Bingham et al.] 2015). A typical choice is tanh, but to make tanh compatible with the distribution type, need to use
h(2) = tanh(z/a).
(3.32)
with a > 0 instead of tanh(z). The choice of will generally depend on whether the distribution is sub- or super-Gaussian with a S for a sub-Gaussian and a 2 for a super-Gaussian. Fundamentally, the scalar a allows controlling the sign of [x=f(z)= - fz()), where f(z) = 2 tanh(z/a/ (see Section 3.2 and Theorem in Hyvärinen & Oja (1998) From reconstruction point of view, the reason why a super- Haper-Gaussian requires I is that it is more likely to be lail-heavy, so, given that is standardized, there will be more values greater 20, which would result in tanh squashing them close 1 (see Fig. 3), making it unlikely to satisfy E(a tanh(z/a)) se IE(2) for the stationary point (see Appendix F.1 Therefore, needs to be just large enough to allow good reconstruction but not too large that the function becomes close to linear. If Σ were trainable, a regularized is more likely work with Gaussian than sub-Gaussian inputs, since will tend
6ff there cigenspaces, then, by virtuc of taking the variances into account, nonlincar PCA will tend to maximize independence within each eigenspace Sj, whereas, reducing all the eigenvectors to have the eigenvalue, linear ICA would maximize independence the entire space
Sm, which might necessarily lead the result.
13

---

Published in Transactions on Machine Learning Research (06/2024)
increase. Other options include setting to be trainable, using max(-a,minzz,d)) approximation of tanh, using an asymmetric function we briefly explored these during experiments in Appendix G
3.5
Mean and variance of data
Mean centring is an important part of PCA (Diamantaras & Kung 1996 Hyvärinen & Oja 2000 Jolliffe & Cadima 2016). So far, we have assumed that the input data is centred; if it is not, then it can casily be centred as & pre-processing step or during training. If the latter, we can explicitly write the loss function in
Eq. 3.11 as
p1y) ((W) E([[x 11 (h(
(x
a [ 1.911111)))) 1g)1)),
(3.33)
where ,4, B(x) and o² (var(x .gg- The mean and variance can be estimated on batch of data, or, is done with batch normalization (loffe & Szegedy 2015), estimated using exponential moving averages (Appendix J.. See Appendix F.8 for non-centred variant.
3.6 Ordering of the components
After minimizing the nonlinear PCA loss, order the components based Σ. It is also possible to automatically induce the order based on index position. One such method is to introduce the projective
deflation-like term P(W) - -NWW [W]g)) into Eq. 3.11 as follows
((W) - (
(3.34)
where is lower triangular without the diagonal. See Appendix F.7 for details and other methods.
4
Experiments
We applied the neural PCA models on image patches and time signals (sce Appendix G for additional experiments);
optimized using gradient descent (see Appendix H for training details).
Image patches We extracted 11 x 11px overlapping patches, with zero padding, from a random subset of 500 images from CIFAR-10 (Krizhevsky et al. 2009), resulting a total of 512K patches. We estimated Σ from the data and we used a (tanhzz/a) with a = 4, though a E [3, 14] also worked well (see Appendix Fig, [2] as a lot of natural images tend have components with super- Gaussian distributions with similar variances (Hyvärinen et al., 2009). We applied linear PCA, nonlinear PCA, and FastICA (Hyvarinen, 1999), and summarized the results in Fig. We see that the PCA filters (Fig. 4aa appear like blurred combinations of the filters found by nonlinear PCA (Fig. 4dh. This stems from the subspace rotational indeterminacy of linear PCA for components with similar variances. We see that a = does not work (Fig. 4dh in comparison to = 4 (Fig. 4dj. By slowly increasing we can have smooth transition passing from lincar PCA filters up to nonlincar PCA filters (sce Appendix Fig. 12). We also see that conventional nonlincar PCA (Fig. 4g) failed lo recover any meaningful filters. This is equally the case when the cintribttion of the decoder not removed (Fig. 4f highlighting that the presence of the decoder contribution is in fact detrimental.
We used
Time signals We adapted classic linear ICA example of signal separation (Comon 1994] from scikit-learn [Pedrogosa et al. 2011). First, we considered an orthogonal transformation. We generated three noisy signals: sinusoidal, square, and sawtooth. The last two had the same variance, and we mixed all three using a random orthogonal mixing matrix. We estimated the variance from data and used a tanh(z/a) with a = 0.8 since the distributions are sub-Gaussian. Figure shows the result of applying lincar and nonlinear PCA. As expected, linear PCA separated the sinusoidal signal, but not the square and sawtooth signals because they had the same variance. Nonlinear PCA recovered the three signals and their variances, up to sign
14

---

Published in Transactions on Machine Learning Research (06/2024)
(a) Linear PCA using SVD
(b) Linear PCA using (c)
Nonlinear PCA
with (d)
Noulinear PCA with
weighted subspace algorithm A(x) tanh(x) (Eq.
3.33
h(x) ttahhx//4)
(Eq.
[Appendix A.3.1
3.33
(e) Linear PCA using sym- (f) Nonlinear PCA (Eq. 3.2
(g)
Conventional nonlinear (h) FastICA with filters or-
metric loss (Eq. 3.2
without stop gradient on the PCA (Eq. 2.22 without dered by norm (inversely pro- decoder
whitening
portional to the standard de- viations)
Figure 4: A of 32 11x11px filters obtained on patches from the CIFAR-10 datasct. We obtained similar filters with the proposed muil-norm-preserving linear PCA loss (b) as the ones obtained via SVD (a). The filters obtained by nonlinear PCA (c, d) seem to have further separated the mixed filters from linear PCA. In particular, this is obvious by looking at the vertical and horizontal line filters at different positions that have been unmixed with nonlinear PCA. We see that setting a = 1 is not enough to separate (c), requiring a large a such as a = 4 (d). We obtained meaningless filters the PCA subspace with the symmetric linear PCA loss (e). We also obtained meaningless filters when the contribution of the decoder was included (f). This is similarly the case with conventional nonlinear PCA without whitening (g). FastICA (h) relaxes the orthogonality assumption of the overall transformation, however, there is some overlap with nonlinear PCA filters.
filters that not necessarily orthogonal;
15

---

Published in Transactions on Machine Learning Research (06/2024)
(a) Observations (mixed signal)
(b) TTme
-5
-5
400
800
1000
1200
400
000
1000
1200
(0) Linear PCA recovered signals
(d) Nonlinear PCA recovered signals
5
0
200
400
800
1000
1200
400
000
800
1000
1200
Figure Three signals (sinusoidal, square, and sawtooth) that were mixed with an orthogonal mixing matrix. Lincar PCA separated the sinusoidal signal as had a distinet variance, but did not separate the square and the sawtooth signals as they had the same variance. Nonlinear PCA separated the signals and recovered their variances.
(a) Observatious (mixed signal)
(b) True sourecs
200
400
600
800
1000
1200
200
400
600
800
1000
1200
(c) Linear PCA recovered signals
(d) Nonlineur PCA recovered signals
2.5
0.0
-25
200
400
600
800
1000
1200
200
400
800
800
1000
1200
(c) 2-layer nonlinear PCA rocuvered signals
(f) FastICA revovered signels
2.5
2.5
0.0
-2.
200
400
600
800
1000
1200
200
400
600
800
1000
1200
Figure 6: Three signals (sinusoidal, square, and sawtooth) that mixed with non-orthogonal mixing matrix. As expected, both linear (c) and nonlinear PCA (d) were not able to recover the signals. Only FastICA (f) and the 2-layer nonlinear PCA model (e) recovered the signals, but without their variances.
indeterminacies. Second, we considered non-orthogonal transformation. We mixed the three signals, this time with distinct variances, using a non-orthogonal mixing matrix. We summarized the results in Fig. 6 As expected, both linear and nonlinear PCA did not recover the signals - they simple recovered the closest orthogonal transformation in terms reconstruction loss. Both FastICA (linear PCA followed by the fast. fixed-point algorithm) and the two-layer nonlinear PCA model (Eq. 3.21 recovered the signals, but without their variances due to the scale indeterminacy.
5 Related work
previously proposed unified model of lincar PCA and lincar ICA theoretical framework from the probabilistic paradigm, where the variances of the components are modelled as separate parameters but are eventually integrated out. Our key idea the same: we model the variances as separate parameters, except that we do not integrate them out, and we flesh out the key role they play, from the neural paradigm, for nonlinear PCA that allows it learn semi-orthogonal transformations that reduce dimensionality rather than just orthogonal transformations.
Hyvärinen (2055)
16

---

Published in Transactions on Machine Learning Research (06/2024)
Neural PCA was started with the seminal work of Oja 1982)
which showed that the first principal component
can be extracted with a Hebbian (Hebb 1949) learning rule. This spurred a lot of interest. in PCA neural
networks in the 80s and 90s, resulting in different variants and extensions for the extraction of multiple
components (Foldiak) 1989 Oja 1989 Rubner & Schulton 1990
Kung & Diamantaras 1990. Oja 1992b;
Bourlard & Kamp 1988: Baldi & Hornik 1989 Oja
1992a, 1995
(see Diarnantaass & Kung (1996) for an
overview), some of which extracted the axis-aligned solution while others the subspace solution. lincar
autoencoder with a mean squared error reconstruction loss was shown to extract the subspace solution
[Burllard & Kamp 1988 Baldi & Hornik 1989
Unlike in the lincar case, the term nonlinear PCA has been applicd more broadly (Hyvärinen & Oja 2000
to single-layer [Xu 1993 Karhunen & Joutsensalo 1994; Oja 1995), multi-layer (Kramer 1991; Oja 1991
Scholz & Vigánio 2002 and kernel-based (Schölkopf et al. 1998) variants. As a single-layer autoencoder, the main emphasis has been its close connection to linear ICA, especially with whitened data (Karhunen) ct al.
1997 1998 Hyvärinen & Oja 2000). Linear ICA is a more powerful extension to PCA in that it sceks to find components that are non-Gaussian and scalistically independent, with the overall transformation not necessarily orthogonal. Many algorithms have been proposed (see Hyvärinen & Oja (2000); Choi et al. Bingham et al. (2015) for an exhaustive overview), with the most popular being FastICA (Hyvarinen 1999 Reconstruction ICA (RICA) (Le et al., 2011 is an autoencoder formulation which combines a linear reconstruction loss with L sparsity; although it was proposed as a method for learning overcomplete ICA features, we show that RICA, in fact, is performing nonlinear PCA, except that it does not preserve unit norm (sce Appendix F.4..
2005
6 Discussion and conclusion
We have proposed o-PCA, unified neural model for linear and nonlinear PCA. This model allows nonlinear PCA to be on equal footing with linear PCA: just like lincar PCA, nonlinear PCA can now learn a semi- orthogonal transformation that reduces dimensionality and orders by variances. But, unlike linear PCA, nonlinear PCA does not suffer from a subspace rotational indeterminacy: it can identify, disentangle, non-Gaussian components that have the same variance.
The mechanism by which nonlinear PCA eliminates the subspace rotational indeterminacy from the canonical lincar PCA solution can be seen directly in its resulting weight update, where there are two terms that appear: one that maximizes variance, and another that maximizes independence.
Previously, it was only possible to apply conventional nonlinear PCA after whitening the input. This meant that when we consider the lincar ICA transformation as a sequence of rotation, scale, rotation, then conventional nonlinear PCA could only be applied to learn the second rolation on unit variance input, but was not able to learn the first rotation. With our model, nonlinear PCA can now be applied to learn not just the second but also the first rotation. Although the nonlinear PCA solution implicitly exists in the linear ICA solution space, it is not necessarily explicitly recovered, because linear ICA has far larger solution space, atlempting to fit a linear transformation of the form WE-IV, whereas nonlinear PCA attempis to fit, WI-1, with W semi-orthogonal and V orthogonal.
An interesting aspect that emerged is that there is an elegant mirroring between linear and nonlinear PCA: the former puts an emphasis on input reconstruction, while the latter on latent rectnn; the former relics on the decoder contribution, while the latter relies on the encoder contribution. With the a-PCA model, we can have a smooth transition between both, rather than a harsh dichotomy. The advantage of nonlinear PCA is that it can be applied wherever linear PCA tends to be applied to learn Jothogonll transformations- - with the added benefit that it can climinate the subspace rotational indeterminacy from the canonical linear PCA solution, making it fully identifiable.
When applied on non- Non-Gaussian inputs once, nonlinear PCA can learn identifiable semi-orthogonal trans- formations that reduce dimensionality, and, when applicd twice, it overall can perform lincar ICA to learn arbitrary unit-variance identifiable linear transformations. Nonlinear PCA, as part of our proposed p-PCA formulation, thus serves as building block for learning linear transformations that are identifiable.
17

---

Published in Transactions on Machine Learning Research (06/2024)
References
Pierre Baldi and Kurt Hornik. Neural networks and principal component analysis: Learning from examples without local minima. Neural networks, 2(1):53-58, 1989.
David Barber. Bayesion reasoning and machine learning. Cambridge University Press, 2012.
Anthony J Bell and Terrence Sejnowski. The "independent components" of natural scenes are edge filters. Vision research, 33(2))33327 3338, 1997.
Yoshua Bengio, Aaron Courville, and Pascal Vincent. Reprosentation learning: A review and new perspectives. IEEE transactions on pollern analysis and machine intelligence, 35(8):1798 1828, 2013.
Ella Bingham, Samuel Kaski, Jorma Laaksonen, and Jouko Lampinen. Advonces in independent component analysis and learning machines. Academic Press, 2015.
Hervé Bourlard and Yves Kamp. Auto-association by multilayer perceptrons and singular value decomposition. Biological cybernetics, (9(4):291-22) 1294,, 1988.
J-F Cardoso and Beate H Laheld. Equivariant adaptive source separation. IEEE Transaclions on signal processing, 44((2):3011-3300 3030, 1996.
Seungjin Choi, Andrzej Cichocki, Hyung-Min Park, and Soo-Young Lee. Blind source separation and independent component analysis: A review. Newral Information Processing-Lettee Letters and Reviews, 6(1):1-57, 2005.
Pierre Comon. Independent component analysis, new concept? Signal processing, 36(3):287-311 -14, 1994.
Konstantinos 1 Diamantaras and Sun Yuan Kung. Principal component neural networks: theory and applications. John Wilcy & Sons, Inc., 1996.
Yuguang Fang, Kenneth A Loparo, and Xiangbo Feng. Inequalities for the trace of matrix product. IEEE Transactions on Automatic Control, (9(12) -490, 1994.
Peter Foldiak. Adaptive network for optimal lincar feature extraction. In Proceedings of IEEE/INNS Int. Joinl. Conf. Neural Nelworks, volume 1, pp. 401-405, 1989.
Harry Horace Harman. Modern factor analysis. University of Chicago press, 1976.
Simon Haykin. Adaptive Filter Theory (3rd Ed.). Prentice-Hall, Inc., USA, 1996.
Donald Olding Hebb. The organization of behaviour. neurophyssological theory, 1949.
Harold Hotelling. Analysis of a complex statistical variables into principal components. Journal of educational psychology, 24(6):417, 1933.
Aapo Hyvärinen. Survey on independent component analysis. Neural Computing Surveys, 2:94 1244, 1999.
Aapo Hyvarinen. Fast and robust fixed- point algorithms for independent component analysis. IEEE transactions Neural Networks, 10(3))626-633) 634, 1999.
Aapo Hyvärinen. unified probabilistic model for independent and principal component analysis. In Advances Independent Component Analysis and Learning Machines, pp. 75-82. Elsevier, 2015.
Aapo Hyvärinen and Erkki Oja. Independent component analysis by general nonlinear hebbian-lice learning rules. signal processing, 64(3):301-313) 1998.
Aapo Hyvärinen and Erkki Oja. Independent component analysis: algorithms and applications. Newral nelworks, 13(4-5):4114445 430, 2000.
Aapo Hyvärinen, Jarmo Hurri, and Patrick O Hoyer. Natural image statistics: probabilistic approach to early computational vision., volume 39. Springer Science & Business Media, 2009.
18

---

Published in Transactions on Machine Learning Research (06/2024)
Sergey loffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. In International Conference on Machine. Learning, pp. 448-456, 2015.
Tan T Jolliffe and Jorge Cadima. Principal component analysis: a review and recent developments. Philosophionl Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374(2065):201 2016.
Christian Jutten and Jeanny Herault. Blind separation of sources, part i: An adaptive algorithm based on neuromimetic architecture. Signal processing, 24(1):1-10, 1991.
Juha Karhunen and Jyrki Joutsensalo. Representation and separation of signals using nonlinear pea type learning. Newral networks, 7(1):113-13 100-111-.17, 1994.
Juha Karhunen and Jyrki Joutsensalo. Generalizations of principal component analysis, optimization problems, and neural networks. Neural Networks, 8(4):549 562, 1995.
Juha Karhunen, Erkki Oja, Liuyue Wang, Ricardo Vigario, and Jyrki Joutsensalo. A class of neural networks for independent component analysis. IEEE Transactions newral networks, 8(3):486 504, 1997.
Juha Karhunen, Petteri Pajunen, and Erkki Oja. The noulinear pea criterion in blind source separation: Relations with other approaches. Newrocompuling, 22(1-3):5 20, 1998.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. ar preprial 0rXiv:141.69960 2014.
Mark A Kramer. Nonlinear principal component analysis using autoassociative neural networks. ATCLE journal, 37(2):233 243, 1991.
Alex Krizhevsky, Geoffrey Hinton, et al. Learning multiple layers of features from tiny images. 2009.
Sun-Yuan Kung and Kostas Diamantaras. neural network learning algorithm for adaptive principal component extraction (apex). In International Conference on Acoustics, Speech, and Signal Processing, pp. 861-864. IEEE, 1990.
Russell Howard Lambert. Multichannel blind deconoolution: FIR matric algebra and separalion of mullippth mintures. University of Southern California, 1996.
Quoc Le, Alexandre Karpenko, Jiquan Ngiam, and Andrew Ng. Ica with reconstruction cost for efficient. overcomplete feature learning. Advances in neural information processing systems, 24, 2011.
Erkki Oja. Simplified neuron model as a principal component analyzer. Journal of mathematical biology, 15 (3):267 273, 1982.
Erkki Oja. Subspace Methods of Pattern Recognition. Electronic & electrical engincering rescarch studies. Research Studies Press, 1983.
Erkki Oja. Neural networks, principal components, and subspaces. International jourraaaoonalauurall of neural systems, (01):61-68, 1989.
Erkki Oja. Data compression, feature extraction, and autoassociation in feedforward neural networks. Arlificial neural networks, 1991.
Erkki Oja. Principal component analysis by homogeneous neural networks, part i: The weighted subspace criterion. IEICE Trans. on Inf. and Syst., 7563):366 375, 1992a.
Erkki Oja. Principal components, minor components, and lincar neural networks. Neural networks, 5(6): 927 935, 1992b.
Erkki Oja. Pca, ica, and nonlinear hebbian learning. In Proc. Int. Conf. on Artificial Neural Networks (ICANN'95), pp. 89 94, 1995.
19

---

Published in Transactions on Machine Learning Research (06/2024)
Erkki Oja, Hidemitsu Ogawa, and Jaroonsakdi Wangviwattana. Principal component analysis by homogeneous neural networks, part ii: Analysis and extensions of the learning algorithms. IEICE Transactions on Information and Systems, 75(3):376 382, 1992.
Bruno A Olshausen and David Field. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature, 361(6383):607 609, 1996.
Karl Pearson. Liii on lines and planes of closest fit to systems of points space. The London, Edinburgh, and Dublin philosophical magazine and journal of science, 2(11):559 572, 1901.
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825 -8330, 2011.
Oren Rippel, Michael Gelbart, and Ryan Adams. Learning ordered representations with nested dropout. In International Conference on Machine Learning, pp- 1754. PMLR, 2014.
Jeanctte Rubner and Klaus Schulton. Development of feature detectors by self-organization: A network model. Biological Cybernelics, 62(3):193- 1990.
Tim Salimans and Durk P Kingma. Weight normalization: A simple reparamcterization accelcrate training of deep neural networks. Advances in neural information processing systems, 29, 2016.
Terence D Sanger. Optimal unsupervised learning in single-layer lincar feedforward neural network. Neural networks, 2(6):459 473, 1989.
Erhard Schmidt. Zur theorie der linearen und nichtlinearen integralgleichunnen. Mathematische Annalen, 63 (4):433 476, 1907.
Bernhard Schölkopf, Alexander Smola, and Klaus-Robert Müller. Nonlinear component analysis as a kernel eigenvalue problem. Neural computation, 10(5):1299-133) 1319, 1998.
Matthias Scholz and Ricardo Vigário. Noulinear poa: a new hierarchical approach. In Esann, pp. 439 444, 2002.
Michael E Tipping and Christopher M Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society Series B: Statistical Methodology, (1(3):611-62) 1999.
Lloyd N Trefethen and David Bau III. Numerical Linear Algebra volume 50. S1AM, 1997.
Liuyue Wang and Juha Karhunen. A unified neural bigradient algorithm for robust pea and mea. International journal of neural systems, 7(01):53 67, 1996.
Liuyue Wang, Juha Karhunen, and Erkki Oja. A bigradient optimization approach for robust, pea, mea, and source separation. In Proceedings of ICNN '95- International Conference on Newral Networks, volume 4, pp. 1684-1689. IEEE, 1995.
R.J. Williams and San Diego Institute for Cognitive Science University of California. Feature Discovery Through Error-correction Learning. ICS report. Institute for Cognitive Science, University of California, San Diego, 1985.
Lei Xu. Least mean square error reconstruction principle for self-organizing neural-nets. Neural networks, 6 (6):627-648, 1993.
Appendix
A Linear PCA
In this section we go through different methods related to PCA and highlight their loss functions and weight updates where relevant. There are a variety of other neural methods for PCA; see (Diamantaras & Kung 1996) for an in-depth review.
20

---

Published in Transactions on Machine Learning Research (06/2024)
A.1 PCA eigendecomposition of covariance matrix SVD of data matrix
Let X E guxp be the centred data matrix and let C = _ X be the symmetric covariance matrix, then we can write the SVD X as
X = USV
(A.1)
where U RexP is semi-orthogonal matrix, S RPx diagonal matrix of singular values, and E Rexp is orthogonal matrix; and we can write the eigendecomposition of C as
C EAE"
(A.2)
where E E RPxP is orthogonal matrix of eigenvectors corresponding to the principal axes, and A E Rexp is a diagonal matrix corresponding to the variances.
We can relate both by noting that we have
C 1 X = U'S
(A.3)
EAE = SiV
(A.4)
n
From this we can see that E = V and A = 1 = diag o²
A.2 PCA subspace solutions
A.2.1 Linear autoencoder
Let X be the centred data matrix and x row of X, then we can write the reconstruction loss of a linear autoencoder as
L(We, Wa) E(((We, Wa)) [[[[xW. WT x11²),
(A.5)
where W. Bpxk is the encoder and WT Rpxk is the decoder. The stochastic loss yields the following gradients:
08
-x W)
(A.6)
OW
at
(XW.WW XWW
(A.7)
OWa
Let us consider the tied weights case, i.c W. We W, but keep the contributions separate. Let y xW, write
21

---

Published in Transactions on Machine Learning Research (06/2024)
dE
XxWW x)W
(A.8)
OW.
-x'xW(W W-I II
(A.9)
=x'y(W W-I I)
(A.10)
0l
(xWWT x)
(A.11)
OW
XW -XWW)) XWW -XxX
(A.12)
=WW x² XW XW- -X'W
(A.13)
=(WW" -1)x'xW
(A.14)
-WWW -I)x
(A.15)
We can note that the term -x"y appears in both the contributions of the encoder and the decoder, and this term is a variance maximization term. Indeed, given that the data is centred, we can write the variance maximization loss
J(W) B w B(11
(A.16)
which has the following gradient (omitting the expectation)
0J oW
-x xW -xx"y.
(A.17)
Thus, the encoder multiplies the variance maximization term with (WTW - I) and the decoder with (WWT - I). If W is constrained to be orthogonal, i.e. WTW I, then the contribution from the encoder is zero, and the contribution from the decoder is dominate. We can also note that the variance maximization of the latent results in a Hebbian update rule AW & x Oja. 1982
Let = w and & = xW, then we can rewrite the contributions as
Oe x y)
(A.18)
OW de
OW = (& - x)T y,
(A.19)
highlighting that the encoder puts an emphasis on latent reconstruction while the docoder puts an emphasis input reconstruction. Rewritten weight updates, have
AW, xx'y xfxW xx-y-x'xW AW Waxxxyy Wy'y. y.
(A.20)
(A.21)
A.2.2 Subspace learning algorithm
The subspace learning algorithm Oja 1983 Williams & University California, 1985 Oja 1989
Hyvärinen
& Oja 2000) was proposed a generalization of Oja's single component neural learning algorithm Ojaa
1982
The weight update has the following form:
22

---

Published in Transactions on Machine Learning Research (06/2024)
AW = x - W(W"x'xx) -x'y-Wy'S y.
(A.22)
(A.23)
This update is obtained by maximizing the variance (minimizing the negative of the variance) under an orthonormality constraint:
((W) (E((W)) -
(A.24)
subject to W W= =
The variance maximizing term results in the gradient
O
xTy
(A.25)
Once combined with an orthonormalization update, this results the following algorithm:
xW W ((11) W + x (0)
(A.26)
W (111) (W+nn(WT ((+1) W ((44)))
(A.27)
where 7 is the learning rate. After expanding both equations as power series of v up to first order, and assuming WT, (() = I, we obtain
(n W(HD)) (I (WT (i+1) ((D
(A.28)
(W + n/x xW (6) w/ (1) xx)xxxWW 1 0(1))
XW
(A.29)
-W xX x 7/ WWoo
xW +
(A.30)
Finally, this results in
AW -x y.
(A.31)
The additional term -WyTy y thus serves as an orthonormalization term
We see that this update exactly the same as Eq. A.21 is simply the decoder contribution from a linear reconstruction loss. As we have seen from Eq. A.9 the contribution of the encoder zero if W is semi-orthogonal, in the subspace learning algorithm the encoder contribution is simply omitted. The subspace learning algorithm can be obtained as a reconstruction loss with the stop gradient operator placed on the encoder:
((W) | WT -x
(A.32)
A.3 Finding axis-aligned principal vectors
Here we look at the main idea for obtaining axis-aligned solutions: symmetry breaking. We saw in the previous section that the weight updates in the linear case are symmetric no particular component is favoured over the other. Therefore, any method that breaks the symmetry will tend to lead to the axis-aligned PCA solution.
23

---

Published in Transactions on Machine Learning Research (06/2024)
A.3.1 Weighted subspace algorithm
One straightforward way to break the symmetry is to simply weigh each component differently. This has been
shown to converge Oja 1992a Oja al. 1992 Xu 1993 Hyvärinen & Oja 2000 to the PCA eigenvectors.
The weighted subspace algorithm does exactly that, and it modifies the update rule of the subspace learning algorithm into
AW x'y yy'yA- AW x Wy"y.
(A.33)
(A.34)
where A = (diag(A1:... Ax) such that d1 > ... > As > O.
The second form can be written as loss function:
( (W)= [[[[x[[],,, ag W - x112 - +|[W)
(A.35)
A.3.2 Weighted subspace algorithm with unit norm
Although the above updates (Eq. A.33 and A.34 do converge to the eigenvectors, they no longer maintain unit norm columns. To see why, we can look at the stationary point, where we know that AW should 0 and that the components of y should be uncorrelated, i.c. y(y"))= Σ is diagonal. We can write:
AW
(A.36)
E(x"y Wy YA -1 yA- 0
(A.37)
E(WT x y W Wyy'' YA -1
(A.38)
E(y"y) y) W WE y)A -1 0
(A.39)
ER W WE A -1= 0
(A.40)
W WA - "22 $
(A.41)
w = A.
(A.42)
As we know W has orthogonal columns, W W must be diagonal, and so the norm of each ith column becomes equal to V we want unit norm columms, a straightforward remedy is to normalize at the end of training. A closer look allows us to see that the main contributor to the norm is the diagonal part of y y -1 this suggests to us two options that could maintain unit norm columns.
As first option, we simply counteract the effect of A-1 the diagonal by multiplying by its inverse the other side of yT to obtain
AW xTy - W(Y1)"yA-1, AW = =x W(yA)
(A.43)
(A.44)
for which we now have
AW = 0
(A.45)
$2 - W W A 0
(A.46)
WWA $2 $2
(A.47)
WTW-1
(A.48)
24

---

Published in Transactions on Machine Learning Research (06/2024)
As second option, remove the diagonal part of y This consists in adding
-Wddiag(diag Soclagg-(0 '(y"))) diag(diag '(y" YA 1))
(A.49)
Eq. A.33 or
Wddiag(diag""(>> (y" yA)) (o
(A.50)
to Eq. A.34
The latter can be derived from the gradient of
(1] sg 1
| (1.
(A.51)
We thus arrive at a total loss that maintains unit norm columns:
f( WW [ W x11
W Σ
A
11
|
1 ||xWA + ||>]
(A.52)
We see that this combines reconstruction, weighted regularization, and weighted variance maximization. We can also note from this that X, should be < 1 to avoid the variance maximization term overpowering the other terms.
A.3.3 Generalized Hebbian Algorithm (GHA)
The generalized Hebbian algorithm (GHA) Sanger 1989
learns multiple PCA components by combining Oja's update rule (Oja 1982) with Gram-Schmid-llik orthogonalisation term. It breaks the symmetry in the subspace learning a algorithm by taking the lower triangular part of y"'y, allowing the weights to converge to the true PCA eigenvectors.
y =
(A.53)
AW x xxyy Wa(y"')
(A.54)
One way to write this as a loss function to use variance maximization and the stop gradient operator to get the following:
((W) = - x xw 2 + +1 o [Wi(yy =) yW
(A.55)
A.3.4 Autoencoder with GHA update
We can also consider combining the orthogonalisation term of GHA with the autoencoder reconstruction to get
((W) = - 1W @ (W y)
(A.56)
where & the operation that takes the lower triangular without the diagonal; otherwise, the term "'diag(diag""y"" y)) would be doubled unless we also include another xTy to counteract it. The double term might still work in practice, but it becomes detrimental when combined with projective unit norm
25

---

Published in Transactions on Machine Learning Research (06/2024)
constraint. If W is orthogonal, then this is similar to combining the subspace update rule (Eq. A.23 with the orthogonalisation term of GHA to result in
AW -x-y--y'' Wy"y Wo(r"y),
(A.57)
which can be obtained from the following loss
((W) |xx x
[W] ag W 1W @ [Wa(y") [Wa(y")))op1)
(A.58)
Nonetheless, keeping the encoder contribution can still help maintain the orthogonality of W, as the term W is similar to the one derived from symmetric orthogonalisation regularizer (see Appendix B
Another variant consider is to make triangular the full gradient of the linear autoencoder update, which includes not just the contribution of the decoder but also that of the encoder, i.e. by laking
Oe aw = =x'y(W' 1) +
(A.59)
and changing into
al aw x'y( - -) x"y y + Wa(y"'). y).
(A.60)
Though the contribution from the encoder is negligible when W is close orthogonal, the term WWTW also acts as an approximation to Gram-Schmidt orthogonalisation (see Appendix B) and when W deviates from being orthogonal, it can help bring it back more quickly.
A.3.5 Nested dropout
Another way to break the symmetry and enforce an ordering is nested dropout (Rippel et al. 2014). a procedure for randomly removing nested sets of components. can be shown (Rippel et al.] 2014) that it leads to an exact equivalence with PCA in the case of single-layer linear autoencoder. An idea similar in spirit was the previously proposed hierarchical PCA method (Scholz & Vigário 2002) for ordering of the principal components; however, it was limited in number of dimensions and was not stochastic. Nested dropout works by assigning a prior distribution p(-) over the component indices 1...k, and then using it to sample an index and drop the units C + 1,..., k. A typical choice of distribution a geometric distribution with p(i) = p - p), 0 <<p<l.
[1 ifi then the reconstruction loss becomes
Let mii E f0, 13 be vector such that m,
otherwise'
((W) = E( (xW o mij W x((2).
(A.61)
where the mask mjlj is randomly sampled during each update step. One issue with nested dropout is that the higher the index the smaller the gradient update. This means that training has to run much longer for them to converge. One suggested remedy for this (Rippel et al., 2014) is to perform gradient sweeping, where once an earlier index has converged, it can frozen and the nested dropout index can be incremented.
A non-stochastic version, which is an extension of the hicrarchical PCA method (Scholz & Vigário 2002 sums all the reconstruction loss to result in
26

---

Published in Transactions on Machine Learning Research (06/2024)
((W) = Σ B[ W - x x1
(A.62)
where ily is the truncated matrix that contains the first j columns of W. However, this can gel more computationally demanding when there is large number of components compared to the stochastic version.
A.3.6 Weighted variance maximization
We can formulate a variant of the weighted subspace algorithm more explicitly as a reconstruction loss combined with weighted variance maximization term:
1 C(W) - xx]]xI -
(A.63)
A consists of linearly spaced values between 0 and 1, and a E [-1,1]. We can also pick A to be proportional to the variances, but the term should not be trainable we can do that with the stop gradient operator: A
where A = diag(A),,... ...,,,,.
B(y²)
max: (B(x)) lag-
The associated gradient of the loss
0l
-x''y Wy y I) yA - x x + Wy y + yww'w--- I). 1).
(A.64)
dW
(A.65)
Following similar analysis as in Appendix A.3.2 it can be shown that at the stationary point we have WTW I+ 2 A.
Alvernative, we can use stochastic weighting by applying nested dropout on the variance regularizer:
((W) B[[[x xww all xW mm 12)
(A.66)
with randomly sampled mask per data sample.
B Enforcing orthogonality
There are few ways enforce orthogonality either via regularizer projective constraint, and this can be symmetric asymmetric.
B.1 Symmetric
A symmetric orthogonalisation scheme can be used when there are no favoured components.
Regularizer A straightforward regularizer, which also enforces unit norm, is
J(W) a - w
(B.1)
and it has the gradient
DJ oW
Wo WID.
(B.2)
27

---

Published in Transactions on Machine Learning Research (06/2024)
Combining the above with PCA ICA updates has been previously referred to the bigradient rule (Wang
et al. 1995 Wang & Karhunen 1996 Karhunen et al. 1997 where is at most I g: for which the justification
can be derived from its corresponding projective constraint.
Projective constraint An iteralive algorithm, as used in FastICA
(Hyvarinen 1999 is the following:
Wi
W;
wi
Repeat until convergence:
W W - wwTw.
3
2
Where the update W WWW stems from the Taylor expansion of W(W' around WTW W We can note that
w 1 WWTW = w W(WTW-D
1
(B.3)
W 1 OJ
(B.4)
8a 8W
which allows to to see that we can set a 8 in Eq. B.2 to derive the projective constraint from the stochastic
gradient descent update.
The iterative algorithm can be shown converge (Hyvarinen) 1999 Hyvärinen & Oja 2000) by analysing the evolution of the eigenvalues of WT W. The initial normalization guarantees that, prior to any iterations,
the eigenvalues of W"W W are I 1. Let ET AE be the eigendecomppsition of WTW, and let consider the general case
W = W WW(W W-I
(B.5)
= =1+ Wwww.
(B.6)
After one iteration we have
- w +(1 w) (2 (WTW3)3
(B.7)
=) AE 2(1 (1++0)8( (3" AE) + (E" AE) = (1++8)) ΛE -((1 8) A'E + "E A
(B.8)
(B.9)
=ET 3)² A 2(1+3)3A* + A
(B.10)
We now have the eigenvalues after one iteration a function of the cigenvalues of the previous itcration:
2(1+ x = (a(A) = (1+ 3)31 + 33
(B.11)
(B.12)
We know that the eigenvalues are the interval J0, 1]. For = we have
fa(1) =
(B.13)
((1) = (1+ - 4(1 888 433 -23.
(B.14)
We that there is stationary point at. B= 2- Noting that fo(A) = can plot fg:
28

---

Published in Transactions on Machine Learning Research (06/2024)
1.5
(c(2)
1
0.5
-0.5
0.5
8
2
We that if ejo, 1] and 33 e)0, 21, then the cigenvalues will always be in J0, 1]; and given that f(A) X, all the eigenvalues will eventually converge to 1.
Thus, for 3 = 1/2 we get the fastest convergence where À does not exceed J0, 1]. This means that when using stochastic gradient descent with a learning rate of and with symmetric regularizer (Eq. B.1], an optimal value for a is 3. This analysis does not consider potential other interactions from the use of additional losses that could affect W the use of an adaptive optimizer like Adam Kingma & Ba (014)).
B.2 Asymmetric
Regularizer From the symmetric regularizer we can derive an asymmetric version that performs Gram- Schmidt-like orthogonalisation. We simply need the lower (or upper) triangular part of W while using stop gradient on either W WT. This means that the loss, which also maintains unit norm, is
1 ( I J(W)=
(B.15)
2
or without the diagonal (to remove the unit norm regularizer):
J(W) (w
(B.16)
where A refers to the lower triangular matrix without the diagonal. From this, we have
0J = W W I).
(B.17)
OW
JJ W
(B.18)
oW
The use of such updale has been previously referred to as the hierarchic version of the bigradient algorithm (Wang et al. 1995).
Recall that in the weight update of GHA (Appendix A.3.3 we have the term Wi(y y); the asymmetric update can benefit from being a similar scale order to be effective, especially when the variances are large. To do this, we can weight the loss by the non-trainable standard deviations to obtain
J(W) = ( (W
1
(B.19)
Gram-Schmidt projective constraint After each gradient update step, we can apply the orthogonalisation procedure (Schmidt 1907].
29

---

Published in Transactions on Machine Learning Research (06/2024)
B.3 Encoder contribution of linear reconstruction
As noted in Appendix A.2.1 in linear autoencoder the contribution of the encoder to the gradient is zero when W is orthogonal, and so used in isolation it can serve as an implicit orthonormality regularizer:
J(W) = B[[|x xW 129)))
(B.20)
having
OJ ow
x yyWWWWII
(B.21)
as gradient. We see that simply replaces in the symmetric regularizer the first W with xTy.
C Enforcing unit norm
If not already taken care of by the orthogonality constraint, there are a few ways to enforce unit norm.
C.1 Projective constraint
The columns w/ of W can be normalized to unit norm after each update step:
(C.1)
C.2 Regularization
We can enforce unit norm by adding a regularizer on the norm of the column vectors w... This results in the loss
J(w!) = (1 112)
(C.2)
which has the following gradient
0J
(1
ww
(C.3)
Ow!
|w
C.3 Differentiable weight normalization
A third way to do this is to use differentiable weight normalization (Salimans & Kingma 2016 except without the scale parameter. This consists in replacing the column vectors with
W
(C.4)
This has the following Jacobian
Ow?" Ow!
w
(C.5)
which results for given loss L in
30

---

Published in Transactions on Machine Learning Research (06/2024)
0c ow? dL
(C.6)
dw? dw? Ow!
Linear ICA
D.1 Overview
The goal of ICA (Jutten & Herault 1991 Comon 1994
Bell Sejnowski 1997: Hyvarinen 1999)
linearly transform the data into set components that are statistically independent as possible. That is, if Rlxk is vector, the goal is to find unmixing matrix B RR such that
y xB
(D.1)
has its components V1,... Us as independent possible.
In an alternative equivalent formulation, we assume that the observed 21, were generated by mixing k independent sources using mixing matrix A, i.e. we have x = SA with A =BB)) to recover the mixing matrix.
(L) and s = y, and the goal is
To be able to estimate B, al least. - of the independent components must have non-Gaussian distributions. Otherwise, if the independent components have Gaussian distributions, then the model is not identifiable.
Without any other assumptions about how the data was mixed, ICA has two ambiguities: it is not possible to determine the variances nor the order of the independent sources.
There are few different algorithms for performing ICA (Hvärinen & Oja 2000 Bingham al. 2015 the most popular being FastICA (Hyvarinen 1999). In the ICA model, the mixing matrix can be arbitrary, meaning that it is not necessarily orthogonal; however, to go about finding it, ICA algorithms generally decompose it into a sequence of transformations that include orthogonal matrices, for there is a guarantee that any matrix has an SVD that encodes the sequence: rolation (orthogonal), scaling (diagonal), rotation (orthogonal). FastICA uses PCA as preprocessing step whiten the data and reduce dimensionality, The whitening operation performs the initial rotation and scaling. What remains after that is to find the last rotation that maximizes non- Non-Gaussianity If USVT is the SVD of the covariance matrix, then the whitening transformation is VS-1, and the resulting overall ICA transformation has the following form
B = VS-
(D.2)
with and W orthogonal matrices and diagonal.
Ordering of components As ICA makes assumption about the transformation, all the ICs assumed to have unit variance. And so there is no order implicd. Without prior knowledge about how the mixing of the sources occurred, it is impossible to resolve the ICA ambiguity. However, if we have reason to assume that the unmixing matrix is close to orthogonal, or simply has unit norm columns (akin to relaxing the orthogonality constraint in PCA while maintaining unit norm directions) then we can in fact order the components by their variances. In this case, the norm of the columns of the unmixing matrix obtained by FastICA inversely proportional to the standard deviations of the components. On the other hand, also possible assume that is the mixing matrix, rather than the unmixing matrix, that has unit norm columns (Hyvärinen 1999). It is also possible to base the ordering on the measure of non- non-Gaussianity [Hväärinen 1999).
D.2 Equivariant adaptive separation via independence
Equivariant adaptive separation via independence (EASI) [Cardoso & Laheld 1996) is serial updating algorithm for source separation. Tt combines both a whitening term
31

---

Published in Transactions on Machine Learning Research (06/2024)
w(y"y-11 1)
(D.3)
and a skew-symmetric term
-(y)h())-h())3>))
(D.4)
obtain the following global relative gradient update rule
AW = - W I+ -+ h(y)"
(D.5)
with 7) the learning rate. The skew-symmetric term originates from skew-symmetrising Wh(y)" in order to roughly preserve orthogonality with each update. To see why, suppose we have WTW = I, and we modify it into W + WE, then we can expand it as
(W + WE)* (W + WE) = I+ ET + + 3.
(D.6)
If want ++ WE to remain orthogonal up first-order, we must also have (W+ WWWE) I+o(E). This implics that & must be skew-symmetric with E E.
D.3 Non-identifiability of a Gaussian distribution
When the criterion used depends on the variance and/or the independence of components, the only case where it is possible to identify a Gaussian distribution is when the transformation is orthogonal and all the variances are clearly distinct. This means that lincar PCA can identify the sources if they have distinct variances and if they were transformed by a rotation. is impossible to identify any other case.
Given y (O, >²) and yA, with x diag(o²), then, from the affine property of the multivariate Gaussian, have
x X~NO
(D.7)
If A is orthogonal, then if all the variances are distinct then is possible to identify up to sign indeterminacy. But components have equal variances, then they are not possible to identify.
If A is non-orthogonal,, then given the symmetry of A AT the spectral thcorem allows us to write A Q" AQ, where Q is orthogonal and A is diagonal. Therefore, there is way to recover the non-orthogonal A even if the sources have distinct variances because it could have been equally transformed by the orthogonal Q.
E Probabilistic PCA and Factor Analysis
Probabilistic PCA (pPCA) Tipping & Bishop 1999) and Factor Analysis (FA) (aaman 1976) are two other methods for learning linear transformations. The former is related to linear PCA whereas the latter in loose sense to linear ICA. The primary difference is that, unlike the standard linear PCA and ICA models, the pPCA and FA models include an additional term for modelling the noise.
Without loss of generality, we will assume that the inputs have zcro mean. Let RIxp be input sample and y E IR. Rxlk its corresponding latent. The probabilistic PCA model is of the form
YA Xe,
(E.1)
32

---

Published in Transactions on Machine Learning Research (06/2024)
where I)) and À is scalar. The covariance of the input can be expressed C =E(x"x) A A + A The pPCA solution takes the form
A = EUT
(E.2)
where Σ = (max ( S is a diagonal matrix of the eigenvalues of C, and U is an orthogonal matrix of the eigenvectors of C. We can express the latents as function of the input using
B B-A-I-UE-
(E.3)
The solution is thus different from that of linear PCA, except that the standard deviations truncated by that of the noise.
FA further generalizes probabilistic PCA by modelling the scale of each component of the noise independently, and so its data model
x= + cA,
(E.4)
where A is diagonal matrix. FA goes about finding both A and A in an itcrative fashion (Barber 2012 This originates from the fact that if we start from initial guess of A, estimated from the variances of the input data, then we can reduce the FA model to that of pPCA by dividing by A to obtain
XA YAA +E.
(E.5)
This can now be solved in the same way as pPCA. This process is repeated a number of times until the estimate of the log likelihood of the data, under the Gaussian assumption, no longer changes. The obtained solution ends up being of the form
A = EUTA,
(E.6)
where Σ = (max(0, S I)) and S and correspond, respectively, to the eigenvalues and eigenvectors of This means that, unlike PCA, the matrix docs not nocessarily have orthogonal rows bocause AA T ΣV" VVV.. Some FA methods proceed with doing an additional rotation by multiplying with another orthogonal matrix V E O(k) to obtain
A -'CAA CA
A - A
(E.7)
This means that the overall transformation of from the inputs to the latents
B A
(F.8)
We this see that, similarly to linear ICA, FA estimates two orthogonal transformations, and so the resulting overall transformation is not necessarily orthogonal. Where they differ is in the presence of an additional scaling of the inputs, and notably, in the process of finding the second rolation, linear ICA maximizes independence, whereas FA maximizes variance (Hyvärinen Oja 2000).
From the point of view learning linear transformation such that y xB, we can summarize the models learnt by the different methods in Tab. 1
33

---

Published in Transactions on Machine Learning Research (06/2024)
B WE
Criterion
Linear PCA
variance
Probabilistic PCA
WE
variance
Nonlinear PCA
WE
variance and independence
Linear ICA WE V
variance (W); independence (V)
Factor Analysis
A W-
variance
Table 1: Linear transformation models learnt by the different methods. Both W and V are orthogonal, whereas both A and E are diagonal. The optimal values found for each of the variables are not necessarily the same. We can note that FA has the most parameters.
F Nonlinear PCA
F.1
Stationary point
F.1.1 Modified derivative
Let us consider the case where h' h1(x)) = 1. This approximation does not seem to affect the learned filters when Σ is estimated from the data (sce Appendix G.4 We can write our gradient as
de aw = x" (byyE"))w'
-) =x=xxx =x=x -)))..
(F.1)
We can take a look at what happens at the stationary point. Taking Eq.
F.1 and assuming that WTW =I, I,
we have
AW = x 1'
(F.2)
0 E(WT, x-y (WW'"W)
(F.3)
0 Σ² ('EY"-")EEWW
(F.4)
² (oyTh(yE-11)
(F.5)
Σ Bly b(y=)
(F.6)
Given that E(y" h(y=)(yE11)) is expected to be diagonal, we end up with o, = E(y,b(y,/,))) be able to satisfy the above, we would need h to adapt such that / /var(h(yE-1)) R We can do this by using a function that adjusts the scale, such as l(x) a tanh(x/a). For instance, if x ~ (0, 1), we have / //va(ltah(z)) & 0.62 and vvar(3 tanh( 0.90. So by setting a value of at least a 3, we get closer back to unit variance.
Alternatively, we can consider an asymmetric nonlinear function such as a ttanh()), where 1/ var(Ianh((z)). If the standardized x is roughly standard normal, then 1.6.
We also have another option instead of scaling: compensate by adding additional loss term to the encoder contribution as follows:
L(W) = B[[xx M(xWE-
(F.7)
+ +OE(IX -
(F.8)
M(x) tanh(x)
(F.9)
6'(x)
(F.10)
34

---

Published in Transactions on Machine Learning Research (06/2024)
This results in
AW = a(xx -
y yW W) + +)'' " WEW"
(F.11)
w W (aE + + "= + byy
(F.12)
hyy
(F.13)
(F.14)
a->xc
And so if set large enough, compensate for the scale. In practice 2 2 is enough to induce effect.
F.1.2 Unmodified derivative
If the derivative of h is not modified, and assuming ha tanh(x/a), (2(x) -22(x)//², the stationary point yields
=== By a
(F.15)
We that the larger is, the more the equality will be satisfied. The choice of will generally depend the type of distribution, which can be seen more clearly in Fig. 3
F.2 Gradient of trainable Σ
From the loss in Eq. 3.3 have
de
diag(diag
(F.16)
0Σ.
0
diag((* xx--x))xe W)
(F.17)
aΣd
where diag creates a diagonal matrix from vector and diag takes the diagonal part of the matrix as vector. If we tie the weights but keep the contributions separate, we obtain
08 diag(diag (y y) M(y
(F.18)
σΣ.
Oi
(F.19)
ΣΣ,
Given Σ = diag(or), can write
Oc
y(y yy)''(yo
(F.20)
di
(ÿ y)h(yo-1)
(F.21)
do
df
-y(s yyo -000- - + y--yo
(F.22)
do
Let consider h tanh, which results in
35

---

Published in Transactions on Machine Learning Research (06/2024)
0R
do - y(8 - y)(1 -) -
- -y)h(yo-)))
(F.23)
(y y) y)(- yo yo 118(000 hlyo
(F.24)
(9 y)(-yo f(yo '))
(F.25)
oly"o -2 yo /(yo - yyyo-yy-
s(yo
(F.26)
alz² - 2f(z) - 2z + 2 (2)),
(F.27)
where f(z) = (x) h(2), a - yo and à - yo 1. We can plot f(2)z, x² and f(z)z
=f(2)
=f(2) 22
1.5
0.5
2
2
Given that y is centred, have E(y2) var(y), and if we take the expectation for the gradient, will obtain
E(z²o) = -ily2 var(y)o
(F.28)
For |a| < we also have =f(2) R 22. And so can see that at takes the form of terms that are roughly do x var(y)o-1 but there is no guarantee that Σ would stop growing unless the reconstruction differences are zero.
F.3 Reconstruction loss in latent space and relationship to blind deconvolution
When performing the reconstruction in latent space, we need to add something that maintains the orthogonality of W. If we multiply both sides of the difference in Eq. 3.11 by [W]sg and include symmetric orthonormality regularizer, obtain
xW
C(W. Σ) = x h(
o Jo Jo W/ (2) + B w
(F.29)
Although we can assume that [T sg = is maintained by an orthonormality regularizer, and might be tempted to simply remove it, we find that keeping it helps make it more stable, especially during the initial
period of training where W might be ioning via a non-orthogonal matrix. Given that [WT W7] is
square, we can also use its lower triangular form N N[""W]e, to obtain an automatic ordering:
XW
L(W, Σ) E([x og hi
oo
W
W W
11)
+
w
(F.30)
36

---

Published in Transactions on Machine Learning Research (06/2024)
If indeed do assume that [WTW]s sy = we find that it generally requires much stronger orthonormality regularizer B >> 1, with an asymmetric regularizer being better than the symmetric, and it works better when E is trainable rather than estimated from data.
C(W, Σ) [[]][[]ee - h(
J (2) + -
WTw
(F.31)
Instead of B||I W W can also use the asymmetric version B A
W W Wesgl (see Appendix B.2
Alternatively, we can use the linear reconstruction loss B(]|xW[W']]] -x (2) for maintaining the orthogonality:
L(W, Σ) = E[[[x h
E x W[WW]
x x112)-
(F.32)
Putting this in context of blind deconvolution
LLmberrt
1996
Haykin
1996), we can note the first term in
Eq. F.29
E(I y o (ylgg11))
(F.33)
can be seen modified form of the Bussgang criterion
B([[(()) - y113.
(F.34)
F.4 Connection with L1 sparsity
ICA is closely related to sparse coding (Olshausen & Field 1996; Bell & Sejnowski 1997 This is because maximizing sparsity can be seen as method for maximizing Gaus-aassyanityy (Hyvärinen & Oja 2000) which is a particular ICA method. Reconstruction ICA (RICA) (Le et al., 2011h is a method that combines L spansity with a linear reconstruction loss, and it was initially proposed for learning overcomplcte sparse Teatures, in contrast to conventional ICA which does not model overcompletennes. RICA is, in fact, simply performing nonlinear PCA: induces a latent reconstruction term the weight update, and, by tying the weights of the encoder and decoder, it enforces W to have orthogonal columns. RICA is thus a special case of the more general ICA model which can learn arbitrary transformation rather than just an orthogonal transformation. The RICA loss is
((W) - xww -
1 + Xx,]'
Σ
(F.35)
where w/' is a column of W. This works equally well if we use the subspace variant, i.c.:
((W) = ||x - x|W],,, W + 3 Σ [ww]] T
(F.36)
If we compute the gradient contributions from Eq. F.36 we obtain
af 3x sign(y))
(F.37)
OW. c
(X x)"
(F.38)
OW,
37

---

Published in Transactions on Machine Learning Research (06/2024)
which results in the combined gradient
de
(&t x)²y 3xx sign(y)))
(F.39)
oW
x"y + x² (Bsign(y) y).
(F.40)
The main relevant part in Eq. F F.40 is
x x (Bsign(y) - y),
(F.41)
and this is similar to the latent, reconstruction term in the nonlinear PCA gradient, except we now have a dependency on the value of 3. No suggestion for the best value of 3 is given by the authors (Le et al.) 2011); however, previous works on sparse coding (Olshausen & Field 1996) have
to be proportional the standard deviation of the input. We can perhaps gain some insight why by looking more closely at the influence Eq. F.40 of the norm of the input. For that extract ||x||2 to get
Bsign(y) - -y ||xx|12( sign(y)
xW )
(F.42)
|xx|22
|x||
We now that if ||x||2 increases, then the norm of sign(y) constant. Indeed, have
will decrease, while that of xW Tx]l2
will remain
sign(y)]l
VE
(F.43)
|
|
||x||a-000
and
X
(F.44)
|x||2
because, from the reconstruction term, we have
xxW xWW -xx xx" ||
(F.45)
And so to make the relative difference By)gn(y) y invariant to the norm of the input, we can set B to be proportional to it, or more simply proportional to the standard deviation of the input as it is already centred, i.e. B = 80 ((() As to what the value of 30 should be, we can plot the functions fro (y) = Bosign(y) y and fa. (y) = Bo tanh(y) y:
(y)
3g 0.25
tanh(y)
2
0.25
Ag
A
( 2.0
y
2
38

---

Published in Transactions on Machine Learning Research (06/2024)
For tanh, if Bo S 1, then it has single inflection point, while if Bo > 1, then it has two inflection points. A reasonable range for So appears to be 1.00. This range is in line with (Olshausen & Field 1996 where 30 was set to 0.14.
The additional term the gradient does not maintain unit norm exactly; this is because, if we follow a similar analysis as Appendix A.3.2 we have for AW
W JEly" sign(y))
(F.46)
if using the subspace RICA loss, or
_ WwW -1 2
Ely sygn
(F.47)
if using the original RICA loss. Therefore, the norm of the columns is less than one.
Another thing we can note is that lima-+o tanh(ax) = sign(x), so in Eq. F.40 we could potentially replace sign with tanh to get
Ol x +x' (Btanh(y) - y).
(F.48)
aw
We know that tanh is the derivative of log cosh, which none other than the function used by FastICA for the negentropy approximation Hyvarinen 1999 This can be expressed the following loss function:
((W) = - x [[
ag W
1 + + Σ
log cosh(xw;)) ).
(F.49)
F.5 Skew-symmetric update
combine the gradient obtained from linear reconstruction loss with the term MyTh(y)) However, to avoid similar problem as in the previous section of the gradient update not maintaining unit norm columns, we can simply remove the diagonal part of h(y), resulting in the weight update
Similar to the derivation of the FAST algorithm (see Appendix D.2
AW xx Wy'y sW(y"h(y) diag(h(y)))))))) Q))..
(F.50)
Or, similar to EASI, use it in its skew-symeetric form:
AW x x - M(y)T))),
-
(F.51)
which we can also write the loss function
((W) = _ 1|x - - W WT 31W o (WyyTh(y) - h(y) y)))l
(F.52)
As justified by the previous section, we can set 8 = [((xx)2), or we can simply use Σ to compensate:
AW x - Wy"y
- diag(y @ 1))
(F.53)
In all of these variants, we have h sign h tanh.
39

---

Published in Transactions on Machine Learning Research (06/2024)
F.6 Modified decoder contribution
The gradient updates of the encoder and decoder (as derived in 3.1 are
di
= y) l1
(F.54)
OW
de
aW (x - (xxx) h(z)E
(F (F.55)
The main issue is that, when both are used in the tied case, the decoder contribution overpowers the encoder contribution. This evident in the lincar case, where, when W semi-orthogonal, the encoder contribution
is zero. Indoed, let W e Roxk with p > k, recall that in the lincar case (Appendix A.2.1 we have
OR = =x W- I)
(F.56)
Ow
OW, (WW" - 1
(F.57)
If we consider their Frobineus norms, write
de || OW.
(F.58)
de ow
Amiu (x (x3y) y)vp
k
||e,
(F.59)
where the lower bound on the decoder contribution can be derived from the fact that for A E Rx and B E Rex,, we have Fang et al. 1994.
Amin
(A)][B| < ||A||| < Amax [|)
(F.60)
where Amin (A) and Amax(A) refer to the smallest and largest eigenvalues of A, respectively, and
WWT 1|; =tr((W" I -I) WrwwW'WW' 2WW +I) = -w' p
(F.61)
(F.62)
(F.63)
In the nonlinear case, the gradient updates are approximately close to the lincar update (especially if h==) = a tanh(z/a) with > 3), except that the encoder contribution is no longer zero, and still overpowered by the decoder contribution. And so if omit the latter, we can better see the effect of the former.
An alternative would be to modify the relative scale that the decoder contribution does not overpower the encoder contribution. The variance of the components plays role in this, given that this is not a problem in the case of conventional nonlinear PCA with whitened input, where all the components have unit variance. This suggests three options: (1) scale the encoder contribution by Σ; (2) drop Σ from the decoder contribution and, optionally, h'(z) from the encoder contribution; (3) scale down the decoder contribution by a constant, which can simply be the inverse of the spectral, Frobineus, or nuclear norm of Σ.
For the third option we can write it a loss:
C(W) E[[ - -X--MXWE-
Xx
[4(WE -IWEE"))....
E) 1
(F.64)
89
Σ
40

---

Published in Transactions on Machine Learning Research (06/2024)
F.7 Ordering the components based on index position
Here we look at few ways for ordering the components automatically based on index position.
F.7.1 Regularizer
One way to do this is via Stam-Schmidt-lite regularizer [Wagg et al. 1995 by adding
J(W) IN ( W
(F.65)
where & refers to the lower triangular matrix without the diagonal element.
F.7.2 Triangular weight update
In the GHA (see Appendix A.3.3), we can order the components by index position by taking the lower (or upper) triangular part of the term Wy"y, y, a term which originates from the linear decoder contribution. Similarly to the GHA, we can also include a triangular term taken from the nonlinear decoder contribution, which is
W(M())))"<()E,,
(F.66)
where yE-1 However, as have seen in Section F.6 need to scale this term appropriately lest it overpowers the encoder contribution. We can consider variations where we drop E from either left or right, or make the approximation of ~ h(z)E:
w ((M2)E)) h(z))
(F.67)
or
WT ((z))
(F.68)
or
w ()(2)Ty)
(F.69)
W (zz((z))
(F.70)
or W (x"y)
(F.71)
or W (y'2) 2)
(F.72)
If E is non-trainable and h = tanh, then Σ = aL where a 2 3 and is the estimated standard deviation of
y.
F.7.3 Weighted latent reconstruction
Another option is similar to the weighted subspace algorithm (see A.3.2 where we insert a linearly-spaced A into the weight update derived Eq. F.T
AW
(F.73)
where A diag (Ai.g()i......aaaa Ax) such that > ... 0.
As loss function we can have:
L(W)- [E(IIy),g
-1yE-1) ΣA W ( W W)
(F.74)
41

---

Published in Transactions on Machine Learning Research (06/2024)
F.7.4 Embedded projective deflation
Let us consider the loss
eWe. Wa. Le Ed) = -1 3xxx-(xW.1P(W.)2 ')&
(F.75)
where
P(W.) = +I Sw 22 Wes.
(F.76)
and W. = Wa = Wa We for elucidating the contribution of each part to the gradient. Taking the gradient of the loss with respect to the weights, we have
08
dWel (((X=
(F.77)
0e - (o<(8-x)WaN((E>"))
(F.78)
OW c2
08
3Way (W ($=x) 1),
(F.79)
OW 63
Setting W. W. W and E. Σg E, obtain
de
T ((9 ()) OY
(F.80)
oW
OL
SW((S - y) o
(F.81)
aw
c2
08
-wy (
(F.82)
aW e3
Let ys = O'yy )) then we can write
x yoP(W)
(F.83)
OW
BWy.y
(F.84)
oW
BWy yo.
(F.85)
aW 63
which results in the overall gradient
de
= xx' y SW M ys))
(F.86)
aW.
If we set P(W) I W W sy: obtain
al
= x ysP(WW) W (yTy).
(F.87)
JW
This introduces an asymmetric term W yyTy), similar to GHA, resulting in an ordering of components by index position.
42

---

Published in Transactions on Machine Learning Research (06/2024)
F.7.5 Nested dropout
ii<j
Similarly to the linear case (Appendix A.3.5,, let mil E 10, be vector such that m, then we can order components using the loss
otherwise'
L(W) [[[ h(xWE-)) m11 E[W"]] ssg + W w
(F.88)
One thing to note about nested dropout is that the larger the index, the less frequently the component receives gradient update, training has to for longer. It also benefit from the addition of an orthonormality regularizer to provide gradient signal to the components receiving less frequent updates from the reconstruction.
F.8 Non-centred nonlinear PCA
When using rero-centred function like tanh, its input should also be zero-centred, so, to oblain a non-cent version of nonlinear PCA, we can simply subtract the mean, and then add it back after passing through the nonlinear function:
xW (W (hl
uy jo + [WTgg] σ
(F.89)
We are simply standardizing (or normalizing to zero mean and unit variance) the components before the nonlincarity and then undoing the standardization after the nonlinearity. This could also be seen as applying a batch normalization layer çlofle & Szegedy 2015) and then undoing afler the nonlinearity.
We can find an upper bound of the loss by writing
) =[([x (N( xW By
ja + [W 1sg (1)
(F.90)
o
xW My
=[[[xx μx + 1x (hl
1o + Myy (T] ]og1
(F.91)
σ
xW (m(---- By (a)l ( + E(lix - Byy Py["" (
< [[[[x fxx - (n(
(F.92)
The upper bound consists of the centred nonlinear PCA loss and a linear mean reconstruction loss. Due to the stop gradient operator, the latter has no effect when W orthogonal (see Appendix A.2.1 instead change to
We can
E(IA Fx WW
(F.93)
resulting into the following non-centred loss:
L(W.o) = E([xx - F - (h( xW My
()) (1 + B(l
(F.94)
43

---

Published in Transactions on Machine Learning Research (06/2024)
G Additional Experiments
G.1 Nonlinear PCA
G.1.1 Highlighting the inability of linear PCA to separate equal variance components
A straightforward way to see how linear PCA is not able to separate components with equal variance is ZCA-whiten the imputs. This involves transforming the imputs with W.Σ ''E-WW which effectively transforms all the components into unit-variance components before projecting them back into the imput space. Figure 7 shows that linear PCA is not able to relearn the original transformation without rotational indeterminacy over the entire space; this is not the case with nonlinear PCA.
(a) Linear PCA filters on unwhitened data
(b) Linear PCA filters ZCA-whitened data
(c) Noulinear PCA filters ZCA-whitened data
Figure 7: PCA filters ZCA-whitened data. Linear PCA is not able to recover the principal eigenvectors since all the variances are equal.
G.1.2 Resolving the linear PCA rotational indeterminacy
Here we atlempt to recover the rotational indeterminacy R from Wr. WvR. For that we use gradient descent to recover R using the loss || WyR + ||RRR - I[]. Figure displays the first 64 linear and nonlincar PCA filters, and Fig. shows the matches for a given lincar PCA filter where the corresponding value in Reyl 0.2.
44

---

Published in Transactions on Machine Learning Research (06/2024)
11:29
7.38
6.7
4.3
3.68
3.3
2.96
2.7
11.21
6.91
6,4
4.54
3.52
3.34
3.3
3.04
2.63
2.21
1.00
1,82
1.8
1.75
1,62
1.62
2.68
1 1.86
1.79
1.79
1.78
1.76
1.07
1.00
1.49
1.36
1:32
1.31
1.3
1.29
1,25
1.19
1.64
1.58
1.52
1.51
1.49
1.44
1.42
1.38
1.14
1.11
1,07
1.06
1.012
11,99
0:99
0,99
1,38
1.34
1.33
1,322
1.17
1.01
1.01
0.59
0.97
0.96
0..96
0.8
0.79
0,77
0.76
0.76
0.97
0.89
0.84
D.8
0.8
0.8
0.8
0.78
0.75
0.56
0.05
0.64
0.64
0.63
0:62
0.61
0.75
0.75
0.74
0.73
0.72
0.61
0.61
0.35
0.32
0.52
0.51
0:51
0.71
0.7
0.69
0.68
0.67
0.07
0.06
0.65
0.5
0.5
0.47
0.42
0.42
0.41
0.41
0.64
0.63
0.02
0.02
0.39
0.58
Linear PCA
(b) Nonlinear PCA
Figure 8: PCA filters ordered by their standard deviations
45

---

Published in Transactions on Machine Learning Research (06/2024)
21 GAMY
1.3
12:1
150
180 5430
1405400
175 CASS
1 24
2016
AN60000
IN6000
18 (4t
(
120/000)
12658.0)
Figure 9: The filters with the red bounding box the top 16 linear PCA filters. The filters adjacent to each one of the filters in red are the nonlincar PCA filters that could be used as linear combination to generate the linear PCA filters. We can note that the filters generally have similar variances, but there are a few outliers which are undoubtedly due to the unconstrained fitting of the rotational indeterminacy matrix R.
46

---

Published in Transactions on Machine Learning Research (06/2024)
G.1.3 Trainable Σ
Epoch 3 without regularization
(b) Epoch 10 - without regularization (c) Epoch 10 with Le Le regular-
ization
Figure 10: Without any regularization, a will tend to keep increasing and result in the degeneration of the filters. This does no occur when using the estimated standard deviations for a.
G.1.4 Scaling of the input
(a) 0.01x
(b) 100x
(c) 1000x
Figure 11: When using nonlincar PCA with estimated standard deviations, it automatically adapts to the scale of the input.
G.1.5 Scale of tanh
47

---

Receandil Losuunng Maccing on Luaabecionna in Fubbbatit
(000)000
C-1 (a)
4-2 (b)
00-3 (0)
-
-
-
10
-hteF 1- (
(
22 -
(k)
26
10
30 =
12: Fiqure start gertiurs berter defimed filters from = 2, but quite a few
Shov-ing the effect of warying the scalar a in
iuaiiciii on
the Obkaiina dilens Audeetreeiiitt
We see dhau We ivvaraaniooooi eatuaageli
we ged
iunenseelaaaatt hotteer soparttion from at least a = 3. The fittes stant leggonatuugg a
14. still are >
48

---

Published in Transactions on Machine Learning Research (06/2024)
0.4
0.6
0.3
0.4
0.4
0.2
0.2
0.1
0.2
0.2
0.0
0.0
0.0
2.5
0.0
2.5
0.0
-
0.0
0.4
0.6
0.4
0.4
0.1
0.2
0.2
0.2
0.2
0.6
0.4
01.4
0.2
0.2
0.2
0.2
Figure 13: The distributions the first 15 filters ordered by their variances. We see that they mostly tend be super-Gass Gaussian distributions, which is why larger value of in atanh(x/a)) was needed.
49

---

Published in Transactions on Machine Learning Research (06/2024)
G.1.6 Latent space reconstruction
Here we look at the reconstruction loss in latent space and look the filters obtained using the losses described in Section
F.3
which we repeal here for convenience:
C(W, Σ) = B([xx[W]sg h XW
)o
W
W
B I w w/ 1gg (B)
F 7.29
h xW
C(W, Σ) = []]
0g (3) + - wTw W
F.30
L(W, Σ) [[[[ 1/ XW
F.31
C(W, Σ) = -<[xx -(XWE)=|1) + [[[[ x(
F F.32
(a) 8=1
(b) B =0
(c) 5 = & T W
Figure 14: Using Eq. F.29 we see that we do not require a strong gain on the orthonormality regularizer.
(a) 8 ==1
(b) 8 = 10
(c) 8 100
Figure 15: Using Eq. F.31
see that it benefits from having a larger gain on the orthonormality regularizer.
50

---

Published in Transactions on Machine Learning Research (06/2024)
Figure Using Eq. F.32 we see that adding the encoder contribution from the linear reconstruction
maintains the orthogonality of the components.
G.1.7 Asymmetric nonlinear function with modified derivative
(a) CIFAR-10
(b) CIFAR-10
((2) = 1.6 tanh(z)
h(=) tanh(2)/ var(tanlh(z)) h'(2) =
k n((x) =1
(c) MNIST
(d) MNIST
h(z) 1.6 tanh(2)
h(2) tanh(z)/ / var(tanh(2))
N'(x)
N' N'(x) =
Figure 17: Filters obtained on CIFAR-10 and MNIST with an asymmetric activation function where a is either constant or adaptive.
G.1.8 First 64 filters
51

---

Published in Transactions on Machine Learning Research (06/2024)
CIFAR-10 symmetric
(b) CIFAR-10 asymmetric
(c) MNIST -symmetric
(d) MNIST asymmetric
Figure 18: The first 64 filters, sorted by variance, obtained on CIFAR-10 and MNIST with cither symmetric activation function (h(z) = 4 lanh(z/4)) or an adaptive asymmetric activation function (N(z) = tanh(2)/ (at(za///varri /var(tanh())), h'()) = 1). We can note that the obtained filters seem to more localized in the asymmetric case compared to the symmetric case.
G.2 Li sparsity RICA
52

---

Published in Transactions on Machine Learning Research (06/2024)
(a) 0.01x, g =
(b) 100x, | 1
1000x, | 1
(d) 0.01x, B 22
(e) 100x, [(
(f) 1000x, 3 E([x
Figure 19: When using RICA F the strength of the sparsity regularizer needs to be adjusted to adapt to the scale of the imput. Here we show the effect of the input scale on the obtained filters between adaptive and nonadaptive 3. (a-c) have / = 1, while (d-f) have 8 = E(I [([xx|2)) We see that making proportional to (2])x makes it invariant to the scale the input.
(a) 8 = 0.01 (0 (12)
(b) 8 = 0.1 ((
(c) / = 1.0F(
(d) 8 x((22)
Figure 20: Filters obtained CIFAR-10 using RICA F.4
with varying L1 sparsity regularization intensity.
Unit weight normalization used (see Appendix C
G.3
Linear PCA
Here we summarize variations of linear PCA methods.
Method
Filters
1: PCA via SVD
A.1
53

---

Published in Transactions on Machine Learning Research (06/2024)
2: Linear autoencoder A.2.1
L(W) = [[[x xww
3: Subspace learning algorithm A.2.2
AW Xxxxy- Wy'y
4: Weighted subspace learning algorithm (WSLA) variant (A.3.1
AW x xxxxy-- - Wy yA -1
12 > 0
5: WSLA variant 2 A.3.1
AW x yA Wy'y y
1 12 \1 > >x>0
6: WSLA variant 3 A.3.1
AW x T - W(yA ) Y
12 A1 >
A
7: Weighted subspace algorithm with unit norm A.3.2
((W) B(I| W
T 2
+ +IW
113
W
- -
12
| :
54

---

Published in Transactions on Machine Learning Research (06/2024)
8: Generalized Hebbian algorithm (GHA) A.3.3
xW A W x x xxyy-Waa. Wa(yTy)
9: GHA with encoder contribution A.3.4
AW x Wi(yr")) y) x xyy((W"W)-
10: GHA + subspace learning algorithm A.3.4
AWa x y - Wa(y"y) y) - -Wly" y - )iag(y²))
11: Reconstruction + GHA A.3.4
((W) E[||x xWWT 1 + 1w @ (Wa(yTy))e,,1?
12: Nested dropout A.3.5
C(W) (x m11 WT -
x(13)
13: Variance maximizer with weighted regularizer A.3.6)
J(W) = - xww
X A
12 0
55

---

Published in Transactions on Machine Learning Research (06/2024)
14: Variance maximizer regularizer with nested dropout A.3.6
C(W)= E[[|x - xWW 12
|| m1151 ()
G.4
Nonlinear PCA
Here we summarize variations nonlinear PCA methods.
Method
Filters
1: Conventional without whitening 2.6
L(W) = [[[[x h(x
2: Trainable without stop gradient 3
C(W, Σ) E([[x N(XWE-/EW'
3: Trainable a (Eq. 3.11
C(W, = - ,,, 13)
4: Trainable a latent reconstruction with reconstruction orthogonal regu-
larizer
F.3
C(W, m) [([xx|[]]s)) ay h(>WE)>|1)
+ W(( xl13)
56

---

Published in Transactions on Machine Learning Research (06/2024)
5: Trainable o latent reconstruction with symmetric orthogonal regularizer
C(W, Σ) = X W h(xWE)=|
B|| w'
B = 10
6: Trainable latent reconstruction with asymmetric orthogonal regularizer
F.3
C(W, Σ) [[[x|[W]s ag h(xWE)=I3)
+ [W][[(WT[)))))))))
10
7: Trainable latent reconstruction with [WTW] sg and symmetric orthog- onal regularizer F.3
C(W, Σ) [([[x W| kg h xW
jo W/
+ || - w
8: EMA tanh
C(W) - (
h(x) tanh(x)
9: EMA or scaled tanh
3.4
") C(W) = 18([x - hi(WW -/X"W""))[W")
h(x)) = Lanh( "
10: EMA o scaled tanh with modified derivative (Eq. F.1
((W) - 1s0
Mxx) tanh( 2
h((x) =1
57

---

Published in Transactions on Machine Learning Research (06/2024)
11: EMA o with hard tanh approximation 3.4
C(W) |8[[x - h(xWX ///W"" 1,01 h(r) max( 22, (2,
(((x) =
12: Asymmetric activation adaptive F.1.1
C(W) E[||x - D/ lag
M(x) tanh(x)
M((x)
1//V//rl//// var(h(z))].
13: Asymmetric activation constant F.1.1
((W) = E[||x - -XW lgg
M(x) 1.6 tanh(x)
h'(x) 1
14: EMA o latent reconstruction with symmetric orthogonal regularizer
scaled tanh
F 3
C(W, Σ) E(Ixx W ag h(xWE)EW'W
W|],,
|| W W
Mxx tanh
15: EMA σ trainable scaled tanh 3.4
C(W, a) [(1 (xxxE- 11) ha (z(z) tanh( )
16: EMA σ baked-in GS. No need to order the components as it is done automatically. (3.6
C(W)- = E([x - MXWP 12g1
P(W) (W)
58

---

Published in Transactions on Machine Learning Research (06/2024)
17: EMA o scaled tanh with nested dropout (F.7.5
C(W) |[[[x - ho-whh((W))))
mir w
E log1
/|W/ -
M(x) tanh
18: EMA a modified tanh derivative with linear reconstruction F.1.1
((W) B([|x 1o "xW"= 1mgl
- >W[WW)=,ll M(x) tanh(x))
h((x) 1
2
19: RICA F.4
C(W) E[[|x www 12 +.3 Σ [xwIll
8 = 12(1x112)
20: Reconstruction with log cosh regularizer F.4
C(W) XxmIxxxxwWTTPP Σ log cosh(xw)])
21: Reconstruction skew symmetric (h(yy)
F.5
<(W) E[[|x - xWw 1
+ 31W O [W h(y) - ((ly)Ty)))a
tanh
0 E([[x|12)
59

---

Published in Transactions on Machine Learning Research (06/2024)
22: Reconstruction with symmetric hy) without diagonal F.5
C(W) B[[|x xWw 1
+ 31W @ (y (VyyThhi h(y) diag(h(y) ,y)llbag
h tanh
8 12
23: Reconstruction with symmetric y² without diagonal F.5
((W) B([x - xww
1W O hy
-iaggh(yy y1))l,
sign
24: EMA σ scaled tanh with triangular update F.7.2
C(W) B[[|x - - ") lagl
1W WW yE )) )
h(x) tanh
G.5
Time Series
G.5.1
Orthogonal same variance
(a) Obscrvations (wixed signal)
(b) True sourees
0
-2
-2
200
800
800
1000
1200
200
400
600
800
1000
1200
[c] Lincar PCA recovered signals
Nonlinear PCA recovered signals
0-
0
600
600
(e) FastICA recovered signals
2
200
100
800
800
1000
1200
Figure 21: Three signals (sinusoidal, square, and sawtooth) that mixed with orthogonal mixing matrix. Linear PCA was unable to separate the signals as they all had the same variance.
60

---

Published in Transactions on Machine Learning Research (06/2024)
G.5.2 Orthogonal all distinct variances
(a) Observations (mixxd signal)
(b) sourees
0-
0-
-5-
800
[c] Linear PCA recovered signals
(d) Nonlinear PCA recovrrd signals
0-
-5
200
100
800
800
1000
1200
200
100
600
1000
1200
[e] FastICA recovered signals
-5-
200
400
800
800
1000
1200
Figure 22: Three signals (sinusoidal, square, and sawtooth) with distinct variances that were mixed with an orthogonal mixing matrix. All three managed to separate the signals.
G.6 Sub- and Super- Gaussian 2D points
Here we look at the effect applying linear PCA, nonlinear PCA, and linear ICA on 2D points (n 1000) sampled from either uniform distribution (sub-Gaussian) or a Laplace distribution (super-Gaussian) Figures
X3 and 24 summarize the results.
61

---

Published in Transactions on Machine Learning Research (06/2024)
1.5
1.-
1.5
0.5
0.5
-0.5
-1.0
-1.5
-1.5
-1.3
-05
0.0
0.5
1.0
1.5
(a) Original
(b) Rotated by x/4
2.0
1.5
1.4
1.0
1.0
14
0.3
0.5
0.0
0.
0.3
-0.5
-10
-1.0-
-2.0
1.5
2.0
(c) Linear PCA
(d) Nonlinear PCA
(e) PastICA
((2) tanh(z)
1.5
1.5
1.5
1.0
1.0
05
055
00
00
0.5
1.0
-100
-10
-1.
-20
-20
0.5
(f) Nonlinear PCA
(g) Nonlinear PCA
(h) Nonlinear PCA
h(2) = 0.5 tznh(2/0.5)
h(2) = 2 tanh(2/2)
((2) = 3 (tanh22/3)
Figure 23: Uniform distribution (sub-aausiian) with equal variance. The original data (a) was rotated by
x/4 (b), then we attempted to recover the original data from the rotated data. We see that linear PCA recover the original data (up to sign and permutational indeterminacies. given the equal variance). We also that
(c) failed to recover the original data. Both nonlinear PCA (d) and FastICA (e) managed
h(z) = a tanh(2/a) with < worked best for the uniform distribution with nonlinear PCA.
62

---

Published in Transactions on Machine Learning Research (06/2024)
0-
Original
(b) Rotated by x/4
4-
0-
(c) Linear PCA
(d) Nonlinear PCA
(e) FastICA
h(z) ttanh2z/3)
4-
2-
-2-
-4
(f) Nonlinear PCA
(g) Nonlinear PCA h(2) = 2 tanh(z/2)
(h) Nonlinear PCA
h(=) = tan)(=)
M22) = 4 tanh(2/44)
Figure 24: Laplace distribution Ssuper-Gaussian with equal variance. The original data (a) was rotated by
</4 (b), then we attempted recover the original data from the rotated data. We see that linear PCA
(c) failed to recover the original data. Both nonlinear PCA (d) and FastICA (e) managed to recover the original data (up to sign and permutational indeterminacies. given the equal variance). We also see that h(2) = a tanh(2/a) with 2 3 work best for the Laplace distribution with nonlinear PCA.
63

---

Published in Transactions on Machine Learning Research (06/2024)
H Training details
Patches We used the Adam optimizer (Kingma & Ba 2014) with learning rate either of 0.01 or 0.001, 31 0.9 and 32 - 0.999. We used & batch size of 128. For the patches, each epoch consisted of 4K iterations, and, unless stated otherwise, we trained for a total of three epochs. With the majority of the nonlinear PCA methods, we used a projective unit norm constraint (see Appendix C The unit norm constraint with the Adam optimizer generally behaved well on unwhitened inputs, but not as well on whitened inputs. With whitened inputs, we found that it was better to either use the vanilla SGD optimizer with unit-norm constraint, the Adam optimizer with unit-norm orthogonality regularizer. This potentially might be due to the adaptive point-wise scaling that the Adam optimizer performs, but we did not investigate this further.
Time signals We used the vanilla SGD optimizer with a learning rate of 0.01 and momentum 0.9. We used a batch size of 100 and trained for total of 200 epochs and picked the weights with the lowest reconstruction loss. When using the Adam optimizer, we noted that it scemed better also use differentiable weight normalization in addition to the projective unit norm constraint.
2D points We used the vanilla SGD optimizer with learning rate of 0.01 and momentum 0.9. We used batch size of 100 and trained for total of 100 epochs.
Block rotation matrix RS = SR
81
1222 0
R R2x2
1)
81
R2x2
R R2x2
1)
($1I2x2
(1)
SR.
-RS
(I.1)
$2
$2
0
$2
J Exponential moving average (EMA)
# Σ yi
(J.1)
0² Σ (y: (yy:-it =2)²
-
(J.2)
ù of + (1 a) j
(J.3)
8² a0² + (1 a) 0
(J.4)
Batch mB(y) y-h 4
(J.5)
Vaa +e
EMA ng(y) y- f
(J.6)
+e
where lag is the stop gradient operator, a is the momentum, and b the batch size.
64


Reviews:
[{"review_id": "4n1J9wxjG1", "reviewer": "Reviewer_5HSt", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths\n- The paper is well written. The introduction and background sections provide a clear discussion of the related area. The connection between linear and nonlinear PCA, as well as linear ICA, is well discussed and insightful.\n- The proposed method is a simple and effective improvement over existing nonlinear PCA method.\n\nWeaknesses\n- While I understand that the existing nonlinear PCA method relies on a whitening step, and that the proposed method is an \"end-to-end\" (or single-step) procedure, the benefit of the latter is not well discussed and does not seem clear. Specifically, does such an end-to-end procedure have a better computational complexity (or run time), or does it lead to better empirical performance (or sample complexity)?\n- It would be great to give a proposition/theorem that explicitly formulates the identifiability conditions and results for the method (see below for more details).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper identifies the reason why conventional nonlinear PCA cannot recover the first rotation in the transformation (which is crucial for dimensionality reduction), and accordingly develops a solution by incorporating a (learnable) diagonal scaling term in the objective function of conventional nonlinear PCA. This results in an \"end-to-end\" nonlinear PCA procedure, without the need of first whitening (i.e., applying linear PCA rotation and standardizing) the data. The paper further discusses about the connection between linear PCA, nonlinear PCA, and linear ICA, and provide empirical studies to validate the proposed method.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- As discussed above, what is the benefit of the proposed \"end-to-end\" nonlinear PCA method, as compared to the existing two-step procedure (whitening + existing nonlinear ICA method)? Does the former have better computational complexity, run time, empirical performance (or sample complexity), or even interpretability?\n- The paper (e.g. conclusion) states that the method has \"added benefit of identifiability\". Although it might be clear to the authors, it would be great to give an explicit statement (e.g. proposition/theorem) to formulate it. E.g. what are the conditions needed for Eq. 23 to achieve identifiability and what is the indeterminacy involved.\n- For Eq. 23, it may be worth discussing whether the global minimum can be reached with the proposed procedure.\n\nSome minor points:\n- It seems that most of the indeterminacies in Sec. 2.1 also hold in factor analysis, which is known to share similarities with PCA and ICA. It is worth making this connection to enrich the paper.\n- Most of the PCA methods discussed in the introduction and related works are quite old. The paper could benefit from also briefly discussing more recent nonlinear PCA methods (such as https://arxiv.org/abs/1402.0119 and https://www.sciencedirect.com/science/article/pii/S089360802300669X).\n- In Sec. 2.4, introduction, or related works, it is stated that \"The only case where. . all distinct\" or \"solely based on ... non-Gaussian\". This could be rephrased because Gaussian components may still be identified in some cases (e.g. https://arxiv.org/abs/1804.00408, https://openreview.net/forum?id=kJIibP5bq2, https://ieeexplore.ieee.org/document/374166), and may be relevant to the identifiability of the proposed method.\n- The motivation of Sec. 3.1 does not seem very clear. It would be great to give a high-level overview in the beginning of Sec. 3.1 to explain the motivation.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "ZSN9MBO2aY", "reviewer": "Reviewer_SXR6", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n\n- *A new model*: This paper introduces σ-PCA, a novel neural model that integrates linear and nonlinear PCA through single-layer autoencoders.  A distinctive aspect of σ-PCA is its ability to handle both rotation and scale, focusing on variances. This approach effectively narrows the gap between linear and nonlinear PCA. \n- *Good properties*: σ-PCA can execute a semi-orthogonal transformation akin to linear PCA, facilitating dimensionality reduction and variance-based ordering. However, it surpasses linear PCA by avoiding the issue of rotational indeterminacy, marking a significant advancement in PCA methodology.\n\n### Weaknesses\n- Connections with modern deep learning progress: I appreciate the new model. However, it is unclear whether the proposed model can be utilized in designing or analyzing modern deep learning models like Transformers and CNNs.\n- Comprehensive analysis of the algorithm: The detailed computational cost of the proposed algorithm is not provided, and thus the scaling ability of the proposed algorithm is unclear.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces σ-PCA, a neural model that unifies linear and nonlinear principal component analysis (PCA) within the framework of single-layer autoencoders. The model is notable for addressing both rotation and scale, emphasizing variances. It provides a comprehensive view, bridging the gap between linear and nonlinear PCA methodologies.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- It is suggested to add discussions about the connections of the proposed models with modern deep learning models.\n- It is suggested to analyze the scalability of the proposed algorithm. \n- It is interesting to explore the robustness especially adversarial robustness of the proposed $\\sigma$-PCA.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "zX5Afsv2wI", "reviewer": "Reviewer_3yA9", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The paper is extremely well written and limpid.\n- It is also one of the clearest expositions of PCA, non-linear PCA and ICA I've come across.\n- Reproducibility is high thanks to the notebooks, which I hope will be posted online.\n- The contribution is simple, consequential and well-motivated.\n\nWeaknesses:\n- The paper takes for granted that the goal of learning a unique rotation in equal variance subspaces of PCA is actually worthwhile. It's not intuitive to the naive reader why this is desirable and the connections to identifyability and statistical independence only come across later as if they were self-evident.\n- The paper proposes an algorithm that is essentially mid-way between PCA -- a robust, well-established technique, useful and widely used in real-world contexts, and ICA -- the poster child for methods that look great theoretically, give great results on synthetic data, but largely fail at any real world applications due to their sensitivity to even mild violations of their underlying assumptions. This begs the question of which category \\sigma-PCA falls into. It would be really strengthening the paper if it provided evidence of the method improving performance on a real-world downstream task, not artificially constructed and not merely producing prettier filters.\n\nUnknowns:\n- I have not looked at the PCA literature in 10 years+ and would not be in a position to evaluate novelty.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents a simple modification of the PCA loss function and optimization algorithm which solves for both the objectives of PCA and non-linear PCA, solving for both variance reduction and statistical independent maximization. It provides extensive analysis of the resulting algorithm on a number of synthetic tasks, and code to reproduce it.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "None mandated. Consider whether a real-world validation on a downstream task can be fit within the scope of this (already lengthy) paper.", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper well studies relationship between, linear PCA, nonlinear PCA, and linear ICA, presenting an interesting nonlinear PCA formulation. Analysis on this relation, in the perspective of SVD of linear transformation is nice. However, the emphasis on trying to relate PCA, nonlinear PCA, and ICA, and developing a method that can simultaneously find the maximal variance subspace while also learning independent components, requires a lot to take in, since each of the mentioned problems has different objectives and underlying assumptions. 
BSS is  the common problem formulation that ties them together? Make it clear what is a common problem. If BSS is a common problem, then  clearly explains what is the benefit of the proposed formulation in solve BSS. 
Since the paper has interesting studies, so I do not want to kill this work.
Decision: Accept (poster)