NEURAL NETWORKS FOR PRINCIPAL COMPONENT ANALYSIS: A NEW LOSS FUNCTION PROVABLY YIELDS ORDERED EXACT EIGENVECTORS
Year: 2019

Abstract:
In this paper, we propose a new loss function for performing principal component analysis (PCA) using linear autoencoders (LAEs). Optimizing the standard L2 loss results in a decoder matrix that spans the principal subspace of the sample covariance of the data, but fails to identify the exact eigenvectors. This downside originates from an invariance that cancels out in the global map. Here, we prove that our loss function eliminates this issue, i.e. the decoder converges to the exact ordered unnormalized eigenvectors of the sample covariance matrix. For this new loss, we establish that all local minima are global optima and also show that computing the new loss (and also its gradients) has the same order of complexity as the classical loss. We report numerical results on both synthetic simulations, and a real-data PCA experiment on MNIST (i.e., a 60, 000×784 matrix), demonstrating our approach to be practically applicable and rectify previous LAEs’ downsides.

heading: 1 INTRODUCTION
text: Ranking among the most widely-used and valuable statistical tools, Principal Component Analysis (PCA) represents a given set of data within a new orthogonal coordinate system in which the data are uncorrelated and the variance of the data along each orthogonal axis is successively ordered from the highest to lowest. The projection of data along each axis gives what are called principal components. Theoretically, eigendecomposition of the covariance matrix provides exactly such a transformation. For large data sets, however, classical decomposition techniques are infeasible and other numerical methods, such as least squares approximation schemes, are practically employed. An especially notable instance is the problem of dimensionality reduction, where only the largest principal components—as the best representative of the data—are desired. Linear autoencoders (LAEs) are one such scheme for dimensionality reduction that is applicable to large data sets.
An LAE with a single fully-connected and linear hidden layer, and Mean Squared Error (MSE) loss function can discover the linear subspace spanned by the principal components. This subspace is the same as the one spanned by the weights of the decoder. However, it failure to identify the exact principal directions. This is due to the fact that, when the encoder is transformed by some matrix, transforming the decoder by the inverse of that matrix will yield no change in the loss. In other words, the loss possesses a symmetry under the action of a group of invertible matrices, so that directions (and orderings/permutations thereto) will not be discriminated.
The early work of Bourlard & Kamp (1988) and Baldi & Hornik (1989) connected LAEs and PCA and demonstrated the lack of identifiability of principal components. Several methods for neural networks compute the exact eigenvectors (Rubner & Tavan, 1989; Xu, 1993; Kung & Diamantaras, 1990; Oja et al., 1992), but they depend on either particular network structures or special optimization methods. It was recently observed (Plaut, 2018; Kunin et al., 2019) that regularization causes the left singular vectors of the decoder to become the exact eigenvectors, but recovering them still requires an extra decomposition step. As Plaut (2018) point out, no existent method recovers the eigenvectors from an LAE in an optimization-independent way on a standard network — this work fills that void.
Moreover, analyzing the loss surface for various architectures of linear/non-linear neural networks is a highly active and prominent area of research (e.g. Baldi & Hornik (1989); Kunin et al. (2019); Pretorius et al. (2018); Frye et al. (2019)). Most of these works extend the results of Baldi & Hornik (1989) for shallow LAEs to more complex networks. However, most retain the original MSE loss, and they prove the same critical point characterization for their specific architecture of interest. Most notably Zhou & Liang (2018) extends the results of Baldi & Hornik (1989) to deep linear networks and shallow RELU networks. In contrast in this work we are going after a loss with better loss surface properties.
We propose a new loss function for performing PCA using LAEs. We show that with the proposed loss function, the decoder converges to the exact ordered unnormalized eigenvectors of the sample covariance matrix. The idea is simple: for identifying p principal directions we build up a total loss function as a sum of p squared error losses, where the ith loss function identifies only the first i principal directions. This approach breaks the symmetry since minimizing the first loss results in the first principal direction, which forces the second loss to find the first and the second. This constraint is propagated through the rest of the losses, resulting in all p principal components being identified. For the new loss we prove that all local minima are global minima.
Consequently, the proposed loss function has both theoretical and practical implications. Theoretically, it provides better understanding of the loss surface. Specifically, we show that any critical point of our loss L is a critical point of the original MSE loss but not vice versa, and conclude that L eliminates those undesirable global minima of the original loss (i.e., exactly those which suffer from the invariance). Given that the set of critical points of L is a subset of critical points of MSE loss, many of the previous work on loss surfaces of more complex networks likely extend. In light of the removal of undesirable global minima through L, examining more complex networks is certainly a very promising direction.
As for practical consequences, we show that the loss and its gradients can be compactly vectorized so that their computational complexity is no different from the MSE loss. Therefore, the loss L can be used to perform PCA/SVD on large datasets using any method of optimization such as Stochastic Gradient Descent (SGD). Chief among the compellingly reasons to perform PCA/SVD using this method is that, in recent years, there has been unprecedented gains in the performance of very large SGD optimizations, with autoencoders in particular successfully handling larger numbers of high-dimensional training data (e.g., images). The loss function we offer is attractive in terms of parallelizability and distributability, and does not prescribe any single specific algorithm or implementation, so stands to continue to benefit from the arms race between SGD and its competitors.
More importantly, this single loss function (without an additional post hoc processing step) fits seamlessly into optimization pipelines (where SGD is but one instance). The result is that the loss allows for PCA/SVD computation as single optimization layer, akin to an instance of a fully differentiable building block in a NN pipeline Amos & Kolter (2017), potentially as part of a much larger network.

heading: 2 THE PROPOSED LOSS FUNCTION AND REVIEW OF FINAL RESULTS
text: Let X ∈ Rn×m and Y ∈ Rn×m be the input and output matrices, where m centered sample points, each n-dimensional, are stacked column-wise. Let xj ∈ Rn and yj ∈ Rn be the jth sample input and output (i.e. the jth column of X and Y , respectively). Define the loss function L(A,B) as
L(A,B) := p∑ i=1 m∑ j=1 ‖yj −AIi;pBxj‖22 = p∑ i=1 ‖Y −AIi;pBX‖2F (1)
where 〈·, ·〉F and ‖·‖F are the Frobenius inner product and norm, Ii;p is a p× p matrix with all elements zero except the first i diagonal elements being one. (Or, equivalently, the matrix obtained by
setting the last p− i diagonal elements of a p×p identity matrix to zero, e.g. I2;3 = [
1 0 0 0 1 0 0 0 0
] .) In
what follows, we shall denote the transpose of matrix M by M ′. Moreover, the matrices A ∈ Rn×p, and B ∈ Rp×n can be viewed as the weights of the decoder and encoder parts of an LAE. The results are based on the following standard assumptions that hold generically:
Assumption 1. For an input X and an output Y , let Σxx := XX ′, Σxy := XY ′, Σyx := Σ′xy and Σyy = Y Y ′ be their sample covariance matrices. We assume
• The input and output data are centered (zero mean). • Σxx, Σxy , Σyx and Σyy are positive definite (of full rank and invertible).
• The covariance matrix Σ := ΣyxΣ−1xxΣxy is of full rank with n distinct eigenvalues λ1 > λ2 > · · · > λn. • The decoder matrix A has no zero columns.
Claim. The main result of this work proved in Theorem 2 is as follows:
If the above assumptions hold then all the local minima of L(A,B) are achieved iff A and B are of the form
A = U1:pDp B = D−1p U ′ 1:pΣyxΣ −1 xx ,
where the ith column of U1:p is the unit eigenvector of Σ := ΣyxΣ−1xxΣxy corresponding to the i th largest eigenvalue and Dp is a diagonal matrix with nonzero diagonal elements. In other words, A contains ordered unnormalized eigenvectors of Σ corresponding to the p largest eigenvalues. Moreover, all the local minima are global minima with the value of the loss function at those global minima being
L(A,B) = p Tr(Σyy)− p∑ i=1 (p− i+ 1)λi,
where λi is the ith largest eigenvalue of Σ := ΣyxΣ−1xxΣxy . In the case of autoencoder (Y = X): Σ = Σxx. Finally, while L(A,B) in the given form containsO(p) matrix products, we will show that it can be evaluated with constant (less than 7) matrix products independent of the value p.

heading: 3 NOTATION
text: In this paper, the underlying field is always R, and positive semidefinite matrices are symmetric by definition. The following constant matrices are used extensively throughout. The matrices Tp ∈ Rp×p and Sp ∈ Rp×p are defined as
(Tp)ij = (p− i+ 1) δij , i.e. Tp = diag (p, p− 1, · · · , 1) , (2)
(Sp)ij=p−max(i, j)+1, i.e. Sp=  p p− 1 · · · 2 1 p− 1 p− 1 · · · 2 1 ... ... . . . 2 1
2 2 2 2 1 1 1 1 1 1
 , e.g. S4= 4 3 2 13 3 2 12 2 2 1
1 1 1 1 . (3) Another matrix that will appear in the formulation is Ŝp := T−1p SpT −1 p . Clearly, the diagonal
matrix Tp is positive definite. As shown in Lemma 2, Sp and Ŝp are positive definite as well.

heading: 4 MAIN THEOREMS
text: The general strategy to prove the above claim is as follows. First the analytical gradients of the loss is derived in a matrix form in Propositions 1 and 2. We compare the gradients with that of the original Minimum Square Error (MSE) loss. Next, we analyze the loss surface by solving the gradient equations which yields the general structure of critical points based on the rank of the decoder matrix A. Next, we delineate several interesting properties of the critical points, notably, any critical point of the loss is also a critical point for the MSE loss but not the other way around. Finally, by performing second order analysis on the loss in Theorem 2 the exact equations for local minima are derived which is shown to be as claimed.
Let L̃(A,B) and L(A,B) be the original loss, and the proposed loss function, respectively, i.e.,
L̃(A,B) := m∑ j=1 ‖yj −ABxj‖22
= ‖Y −ABX‖2F
L(A,B) := p∑ i=1 m∑ j=1 ‖yj −AIi;pBxj‖22
= p∑ i=1 ‖Y −AIi;pBX‖2F
The first step is to calculate the gradients with respect to A and B and set them to zero to derive the implicit expressions for the critical points. In order to do so, first, in Lemma 5, for a fixed A, we derive the directional (Gateaux) derivative of the loss with respect to B along an arbitrary direction W ∈ Rp×n, denoted as dBL(A,B)W , i.e.
dBL(A,B)W = lim‖W ‖F→0 L(A,B + W )− L(A,B) ‖W ‖F .
As shown in the proof of the lemma, dBL(A,B)W is derived by writing the norm in the loss as an inner product, opening it up using linearity of inner product, dismiss second order terms in W (i.e. O(‖W ‖2)) and rearrange the result as the inner product between the gradient with respect to B, and the direction W , which yields
dBL(A,B)W = −2 Tr (W ′ (TpA′Σyx − (Sp ◦ (A′A))BΣxx)) = −2〈TpA′Σyx − (Sp ◦ (A′A))BΣxx,W 〉F , (4)
where, ◦ is the Hadamard product and the constant matrices Tp and Sp, were defined in the beginning. Second, the same process is done in Lemma 6, to derive dAL(A,B)V ; the derivative of L with respect to A in an arbitrary direction V ∈ Rn×p, for a fixed B, which is then set to zero to derive the implicit expressions for the critical points. The results are formally stated in the two following propositions. Proposition 1. For any fixed matrix A ∈ Rn×p the function L(A,B) is convex in the coefficients of B and attains its minimum for any B satisfying the equation
(Sp ◦ (A′A))BΣxx = TpA′Σyx, (5) where ◦ is the Hadamard (element-wise) product operator, and Sp and Tp are constant matrices defined in the previous section. Further, if A has no zero column, then L(A,B) is strictly convex in B and has a unique minimum when the critical B is
B = B̂(A) = (Sp ◦ (A′A))−1TpA′ΣyxΣ−1xx , (6) and in the autoencoder case it becomes
B = B̂(A) = (Sp ◦ (A′A))−1TpA′. (6′)
The proof is given in appendix A.2. Remark 1. Note that as long as A has no zero column, Sp ◦ (A′A) is nonsingular (we will explain the reason soon). In practice, A with zero columns can always be avoided by nudging the zero columns of A during the gradient decent process. Proposition 2. For any fixed matrix B ∈ Rp×n the function L(A,B) is a convex function in A. Moreover, for a fixed B, the matrix A that satisfies
A (Sp ◦ (BΣxxB′)) =ΣyxB′Tp (7) is a critical point of L(A,B).
The proof is given in appendix A.3.
The pair (A,B) is a critical point of L if they make dBL(A,B)W and dAL(A,B)V zero for any pair of directions (V ,W ). Therefore, the implicit equations for critical points are given below, next to the ones derived by Baldi & Hornik (1989) for L̃(A,B).
For L̃(A,B):
A′ABΣxx = A ′Σyx, ABΣxxB ′ = ΣyxB ′.
For L(A,B):
(Sp ◦ (A′A))BΣxx = TpA′Σyx, A (Sp ◦ (BΣxxB′)) = ΣyxB′Tp.
Remark 2. Notice the similarity, and the difference only being the presence of Hadamard product by Sp in the left and by diagonal Tp in the right. Therefore, practically, the added computational cost of evaluating the gradients is negligible compare to that of MSE loss.
The next step is to determine the structure of (A,B) that satisfies the above equations, and find the subset of those solutions that account for local minima. For the original loss, the first expression (A′ABΣxx = A′Σyx) is used to solve for B and put it in the second to derive an expression solely based on A. Obviously, in order to solve the first expression for B, two cases are considered separately: the case where A is of full rank p, so A′A is invertible, and the case of A being of rank r < p. Here we do the same but there is a twist; for us there is only one case. The reason is as long as (not necessarily full rank) A has no zero column, Sp ◦ (A′A) is positive definite and hence, invertible. This is discussed in detail in Lemma 2 and we briefly explain it here. As shown in the lemma, Sp is positive definite and by Shur product theorem for any A (of any rank), Sp ◦ (A′A) is positive semidefinite. However, as a result of Oppenheim inequality (Horn & Johnson (2012), Thm 7.8.16), that in our case translates to det(Sp) ∏ i(A ′A)ii ≤ det(Sp ◦ (A′A)), as long as A has no
zero column, ∏ i(A ′A)ii > 0 and therefore, det(Sp ◦ (A′A)) > 0. Here, we assume A of any rank r ≤ p has no zero column (since this can be easily avoided in practice) and consider Sp ◦ (A′A) to be always invertible. Therefore, (A,B) define a critical point of losses L̃ and L if
For L̃(A,B) and full rank A:
B = B̂(A) = (A′A)−1A′ΣyxΣ −1 xx , ABΣxxB ′ = ΣyxB ′.
For L(A,B) and no zero column A:
B = B̂(A) = (Sp ◦ (A′A))−1TpA′ΣyxΣ−1xx , A (Sp ◦ (BΣxxB′)) = ΣyxB′Tp.
Before, we state the main theorem we need the following definitions. First, a rectangular permutation matrix Πr ∈ Rr×p is a matrix that each column consists of at most one nonzero element with the value 1. If the rank of Πr is r with r < p then clearly, Πr has p − r zero columns. Also, by taking away those zero columns the resultant r × r submatrix of Πr is a standard square permutation matrix.
Second, under the conditions provided in Assumption 1, the matrix Σ := ΣyxΣ−1xxΣxy has an eigenvalue decomposition Σ = UΛU ′, where the ith column of U , denoted as ui, is an eigenvector of Σ corresponding to the ith largest eigenvalue of Σ, denoted as λi. Also, Λ = diag(λ1, · · · , λn) is the diagonal vector of ordered eigenvalues of Σ, with λ1 > λ2 > · · · > λn > 0. We use the following notation to organize a subset of eigenvectors of Σ into a rectangular matrix. Let for any r ≤ p, Ir = {i1, · · · , ir}(1 ≤ i1 < · · · < ir < n) be any ordered r−index set. Define UIr ∈ Rn×p as UIr = [ui1 , · · · ,uir ]. That is the columns of UIr are the ordered orthonormal eigenvectors of Σ associated with eigenvalues λi1 < · · · < λir . Clearly, when r = p, we have UIr = [ui1 , · · · ,uip ] corresponding to an p−index set Ip = {i1, · · · , ip}(1 ≤ i1 < · · · < ip < n). Similarly, we define ΛIr ∈ Rp×p as ΛIr = diag(λi1 , · · · , λir ). Theorem 1. Let A ∈ Rn×p and B ∈ Rp×n such that A is of rank r ≤ p. Under the conditions provided in Assumption 1 and the above notation, The matrices A and B define a critical point of L(A,B) if and only if for any given r−index set Ir, and a nonsingular diagonal matrix D ∈ Rr×r, A and B are of the form
A = UIrCD, (8) B = B̂(A) = D−1ΠCU ′ IrΣyxΣ −1 xx , (9)
where, C ∈ Rr×p is of of full rank r with nonzero and normalized columns such that ΠC := (Sp ◦ (C ′C))−1 TpC ′ is a rectangular permutation matrix of rank r and CΠC = Ir. For all 1 ≤ r ≤ p, such C always exists. In particular, if matrix A is of full rank p, i.e. r = p, the two given conditions on ΠC are satisfied iff the invertible matrix C is any squared p × p permutation matrix Π. In this case (A,B) define a critical point of L(A,B) iff they are of the form
A = UIpΠD, (10)
B = B̂(A) = D−1Π′U ′IpΣyxΣ −1 xx . (11)
The proof is given in appendix A.4.
Remark 3. The above theorem provides explicit equations for the critical points of the loss surface in terms of the rank of the decoder matrixA and the eigenvectors of Σ. This explicit structure allows us to further analyze the loss surface and its local/global minima.
Here, we provide a proof sketch for the above theorem to make the claims more clear. Again as a reminder, the EVD of Σ := ΣyxΣ−1xxΣxy is Σ = UΛU
′. For both L̃ and L, the corresponding B̂(A) is replaced by B on the RHS of critical point equations. For the loss L(A,B), as shown in the proof of the theorem, results in the following identity
U ′A ( Sp ◦ ( B̂ΣxxB̂ ′ )) A′U = Λ∆, (12)
where ∆ := U ′ATp(Sp ◦ (A′A))−1TpA′U is symmetric and positive semidefinite. The LHS of eq. (12) is symmetric so the RHS is symmetric too, so Λ∆ = (Λ∆)′ = ∆′Λ′ = ∆Λ. Therefore, ∆ commutes with the diagonal matrix of eigenvalues Λ. Since eigenvalues are assumed to be distinct, ∆ has to be diagonal as well. By Lemma 2 Tp(Sp ◦ (A′A))−1Tp is positive definite and U is an orthogonal matrix. Therefore, r = rank(A) = rank(∆) = rank(U ′∆U), which implies that the diagonal matrix ∆, has r nonzero and positive diagonal entries. There exists an r−index set Ir corresponding to the nonzero diagonal elements of ∆. Forming a diagonal matrix ∆Ir ∈ Rr×r by filling its diagonal entries (in order) by the nonzero diagonal elements of ∆, we have
U∆U ′ = UIr∆IrU ′ Ir Def of ∆ ====⇒
ATp(Sp ◦ (A′A))−1TpA′ = UIr∆IrU ′Ir , (13) which indicates that the matrix A has the same column space as UIr . Therefore, there exists a full rank matrix C̄ ∈ Rr×p such that A = UIrC̄. Since A has no zero column, C̄ has no zero column. Further, by normalizing the columns of C̄ we can write A = UIrCD, where D ∈ Rp×p is diagonal that contains the norms of columns of C̄.
Baldi & Hornik (1989) did something similar for full rank A for the loss L̃ to derive (AL̃ = UIpC̃). But their C̃ can be any invertible p × p matrix. However, in our case, the matrix C ∈ Rr×p corresponding to rank r ≤ p matrix A, has to satisfy eq. (13) by replacing A by UIrCD and eq. (12) by replacing B̂(A) by B̂(UIrCD). In the case of Baldi & Hornik (1989), for the original loss L̃, equations similar to eq. (13) and eq. (12) appear but they are are satisfied trivially by any invertible matrix C̃. Simplifying those equations by using A = UIrCD after some algebraic manipulation results in the following two conditions for C:
CTp (Sp ◦ (C ′C))−1 TpC ′ =∆Ir and (14) C ( Sp ◦ ( (Sp ◦ (C ′C))−1TpC ′ΛIrCTp(Sp ◦ (C ′C))−1 )) C ′ =ΛIr∆Ir . (15)
As detailed in proof of Theorem 1, solving for C leads to its specific structure as laid out in the theorem. Remark 4. Note that when A is of rank r < p with no zero columns then the invariant matrix C is not necessarily a rectangular permutation matrix but ΠC := (Sp ◦ (C ′C))−1 TpC ′ is a rectangular permutation matrix with CΠC = Ir. It is only when r = p that the invariant matrix C becomes a permutation matrix. Nevertheless, as we show in the following corollary, the global map is always ∀r ≤ p : G = AB = UIrU ′IrΣyxΣ−1xx . It is possible to find further structure (in terms of block matrices) for the invariant matrix C when r < p. However, this is not necessary as we soon show that all rank deficient matrix As are saddle points for the loss and ideally should be passed by during the gradient decent process. Based on some numerical results our conjecture is that when r < p the matrix C can only start with a r × k rectangular permutation matrix of rank r with r ≤ k ≤ p and the rest of p− k columns of C is arbitrary as long as none of the columns are identically zero. Corollary 1. Let (A,B) be a critical point of L(A,B) under the conditions provided in Assumption 1 and rankA = r ≤ p. Then the following are true
1. The matrix BΣxxB′ is a p× p diagonal matrix of rank r.
2. For all 1 ≤ r ≤ p, for any critical pair (A,B), the global map G := AB becomes G = UIrU ′ IrΣyxΣ −1 xx . (16)
For the autoencoder case (Y = X) the global map is simply G = UIrU ′ Ir .
3. (A,B) is also the critical point of the classical loss L̃(A,B) = ∑p i=1 ‖Y −ABX‖ 2 F .
The proof is given in appendix A.5. Remark 5. The above corollary implies that L(A,B) not only does not add any extra critical point compare to the original loss L̃(A,B), it provides the same global map G := AB. It only limits the structure of the invariance matrix C as described in Theorem 1 so that the decoder matrix A can recover the exact eigenvectors of Σ. Lemma 1. The loss function L(A,B) can be written as L(A,B) = pTr(Σyy)− 2 Tr (ATpBΣxy) + Tr (B′ (Sp ◦ (A′A))BΣxx) . (17) The above identity shows that the number of matrix operations required for computing the loss L(A,B) is constant and thereby independent of the value of p.
The proof is given in appendix A.6. Theorem 2. Let A∗ ∈ Rn×p and B∗ ∈ Rp×n such that A∗ is of rank r ≤ p. Under the conditions provided in Assumption 1, (A∗,B∗) define a local minima of the proposed loss function iff they are of the form
A∗ = U1:pDp (18) B∗ = D−1p U ′ 1:pΣyxΣ −1 xx , (19)
where the ith column of U1:p is a unit eigenvector of Σ := ΣyxΣ−1xxΣxy corresponding the i th largest eigenvalue and Dp is a diagonal matrix with nonzero diagonal elements. In other words, A∗ contains ordered unnormalized eigenvectors of Σ corresponding to the p largest eigenvalues. Moreover, all the local minima are global minima with the value of the loss function at those global minima being
L(A∗,B∗) = p Tr(Σyy)− p∑ i=1 (p− i+ 1)λi, (20)
where λi is the ith largest eigenvalue of Σ.
The proof is given in appendix A.7. Remark 6. Finally, the second and third assumptions we made in the beginning in Assumption 1 can be relaxed by requiring only Σxx to be full rank. The output data can have a different dimension than the input. That is Y ∈ Rn×m and X ∈ Rn′×m, where n 6= n′. The reason is that the given loss function structurally is very similar to MSE loss and can be represented as a Frobenius norm on the space of n × m matrices. In this case the covariance matrix Σ := ΣyxΣ−1xxΣxy is still n × n. Clearly, for under-constrained systems with n < n′ the full rank assumption of Σ holds. For the overdetermined case, where n′ > n the second and third assumptions in Assumption 1 can be relaxed: we only require Σxx to be full rank since this is the only matrix that is inverted in the theorems. Note that if p > min(n′, n) then ΛIp : the p× p diagonal matrix of eigenvalues of Σ for a p-index-set Ip bounds to have some zeros and will be say rank r < p, which in turn, results in the encoder A with rank r. However, the Theorem 1 is proved for encoder of any rank r ≤ p. Finally, following theorem 2 then the first r columns of the encoder A converges to ordered eigenvectors of Σ while the p − r remaining columns span the kernel space of Σ. Moreover, Σ need not to have distinct eigenvectors. In this case ∆Ir becomes a block diagonal matrix, where the blocks correspond to identical eigenvalues ΣIr . In this case, the corresponding eigenvectors in A
∗ are not unique but they span the respective eigenspace.

heading: 5 EXPERIMENTS

heading: 5.1 EXPERIMENTAL SETUP
text: LAEs with Two Loss functions We will verify the loss function L(A,B) defined in eq. (1) by setting the input matrix X ∈ Rn×m equal to the output matrix Y ∈ Rn×m (Y = X), where the
linear autodecoder (LAE) becomes a solution to PCA. In order for comparison, we train another LAE using the MSE loss L̃(Ã, B̃) defined as L̃(Ã, B̃) = ∥∥∥Y − ÃB̃X∥∥∥2
F , where Y = X is also
applied in our experiments.
The weights of networks are initialized to random numbers with a small enough standard deviation (10−7 in our case). We choose to use the Adam optimizer with a scheduled learning rate (starting from 10−3 and ending with 10−6 in our case), which empirically benefits the optimization process. The two training processes are stopped at the same iteration at which one of the models firstly finds all of the principal directions. As a side note, we feed all data samples to the network at one time with batch size equal to m, although mini-batch implementations are apparently amendable.
Evaluation Metrics We use the classical PCA approach to get the ground truth principal direction matrix A∗ ∈ Rn×p, by conducting Eigen Value Decomposition (EVD) to XX ′ ∈ Rn×n or Singular Value Decomposition (SVD) to X ∈ Rn×m. As a reminder, A ∈ Rn×p stands for the decoder weight matrix of an trained LAE given a loss function L. To measure the distance between A∗ and A, we propose an absolute cosine similarity (ACS) matrix inspired by mutual coherence (Donoho et al., 2005), which is defined as:
ACSij = |〈A∗i ,Aj〉| ‖A∗i ‖ · ‖Aj‖ , (21)
where A∗i ∈ Rn×1 denotes the ith ground truth principal direction, and Aj ∈ Rn×1 denotes the jth column of the decoder A, i, j = 1, 2, . . . , p. The elements of ACS ∈ Rp×p in eq. (21) take values between [0,1], measuring pair-wise similarity across two sets of vectors. The absolute value absorbs the sign ambiguity of principal directions.
The performances of LAEs are evaluated by defining the following metrics:
RatioTP = p∑ i=1 I[ACSii > 1− ]/p (22)
RatioFP = p∑ i,j=1 i6=j I[ACSij > 1− ]/p, and (23)
RatioTotal = RatioTP + RatioFP , (24)
where I is the indicator function and is a manual tolerance threshold ( = 0.01 in our case). If two vectors have absolute cosine similarity over 1 − , they are deemed equal. Considering some columns of decoder may be correct principal directions but not in the right order, we introduce RatioTP and RatioFP in eqs. (22) and (23) to check the ratio of correct in-place and out-of-place principal directions respectively. ThenRatioTotal in eq. (24) measures the total ratio of the correctly obtained principal directions by the LAE regardless of the order.
Datasets As a proof-of-concept, both synthetic data and real data are used. For the synthetic data, 2000 zero-centered data samples are generated from a 1000-dimension zero mean multivariate normal distribution with the covariance matrix being diag(Np). For the real data, we choose to use MNIST dataset (LeCun et al., 1998), which includes 60,000 grayscale handwritten digits images, each of dimension 28× 28 = 784.

heading: 5.2 EVALUATION AND ANALYSIS
text: Synthetic Data Experiments In our experiment, p, the number of desired principal components (PCs), is set to 100, i.e. the dimension is to be reduced from 1000 to 100. Figures 1 and 2 demonstrate a few conclusions. First, during the training process, the loss ratio of both losses continuously decreases to 1, i.e. they both converge to the optimal loss value. However, when both get close enough, L require more iterations since the optimizer is forced to find the right directions: it fully converges only after it has found all the principal directions in the right order.
Second, using the loss L results in finding more and more correct principal directions, with RatioTP continuously rising; and ultimately affords all correct and ordered principal directions,
with RatioTP ending with 100%. Notice that occasionally and temporarily, some of the principal directions is found but not at their correct position, which is indicated by the rise of RatioFP in the figure. However, as optimization continues they are shifted to the right column, which results in RatioFP going back to zero, and RatioTP reaching one. As for L̃, it fails to identify any principal directions; both RatioTP and RatioFP for L̃ stay at 0, which indicates that none of the columns of the decoder Ã, aligns with any principal direction.
Third, as shown in the figure, while the optimizer finds almost all the principal directions rather quickly, it requires much more iterations to find some final ones. This is because some eigenvalues in the empirical covariance matrix of the finite 2000 samples become very close (the difference becomes less than 1). Therefore, the loss has to get very close to the optimal loss, making the gradient of the loss hard to distinguish between the two.
Real Data: MNIST Experiments We set the number of principal components (PCs) as 100, i.e., the dimension is to be reduced from 784 to 100. We also try to reconstruct with the top-10 columns found in this case. As in Fig. 3, the reconstruction performance of L is consistently better than L̃. That also reflects that L̃ does not identify PCs, while L is directly applicable to performing PCA without bells and whistles.

heading: 6 CONCLUSION
text: In this paper, we have introduced a loss function for performing principal component analysis and linear regression using linear autoencoders. We have proved that the optimizing with the given loss results in the decoder matrix converges to the exact ordered unnormalized eigenvectors of the sample covariance matrix. We have also demonstrated the claims on a synthetic data set of random samples drawn from a multivariate normal distribution and on
MNIST data set. There are several possible generalizations of this approach we are currently working on. One is improving performance when the corresponding eigenvalues of two principal directions are very close and another is generalization of the loss for tensor decomposition.

heading: APPENDIX

heading: A PROOFS

heading: A.1 PRELIMINARIES
text: Before we present the proof for the main theorems, the following two lemmas introduce some notations and basic relations that are required for the proofs. Lemma 2. The constant matrices Tp ∈ Rp×p and Sp ∈ Rp×p are defined as (Tp)ij = (p− i+ 1) δij , i.e. Tp = diag (p, p− 1, · · · , 1) ,
(Sp)ij = p−max(i, j) + 1, i.e. Sp =  p p− 1 · · · 2 1 p− 1 p− 1 · · · 2 1 ... ... . . . 2 1
2 2 2 2 1 1 1 1 1 1
 , e.g. S4 = 4 3 2 13 3 2 12 2 2 1
1 1 1 1  . Clearly, the diagonal matrix Tp is positive definite. Another matrix that will appear in the formulation is Ŝp := T−1p SpT −1 p
( Ŝp ) ij = ( T−1p SpT −1 p ) ij = 1 p−min(i, j) + 1 i.e. T −1 p SpT −1 p =  1 p
1 p · · · 1p 1p
1 p 1 p−1 · · · 1p−1 1p−1 ... ... . . . ...
... 1 p 1 p−1 · · · 12 12 1 p 1 p−1 · · · 12 1
 ,
e.g. Ŝ4 =  1 4 1 4 1 4 1 4 1 4 1 3 1 3 1 3
1 4 1 3 1 2 1 2 1 4 1 3 1 2 1  . The following properties of Hadamard product and matrices Tp and Sp are used throughout:
1. For any arbitrary matrix A ∈ Rn×p, p∑ i=1 Ii;p = Tp, and (25)
p∑ i=1 Ii;pA ′AIi;p = Sp ◦ (A′A) , (26)
where, ◦ is the Hadamard (element-wise) product. 2. For any matrices M1,M2 ∈ Rp×p and diagonal matrices D ,E ∈ Rp×p,
D (M1 ◦M2)E = (DM1E ) ◦M2 = M1 ◦ (DM2E ) . Moreover, if Π1,Π2 ∈ Rp×p are permutation matrices then
Π1 (M1 ◦M2) Π2 = (Π1M1Π2) ◦ (Π1M2Π2) .
3. Sp is invertible and its inverse is a symmetric tridiagonal matrix
(S−1p )ij =  1 i = j = 1
2 i = j 6= 1 −1 |i− j| = 1 0 otherwise , i.e. S−1p =
 1 −1 · · · 0 0 −1 2 −1 0 0 ... ... . . . ...
... 0 0 −1 2 −1 0 0 0 −1 2
 .
4. Sp is positive definite.
5. For any matrix A ∈ Rn×p, Sp ◦ (A′A) is positive semidefinite. If (not necessarily full rank) A has no zero column then Sp ◦ (A′A) is positive definite.
6. For any diagonal matrix D ∈ Rp×p
Sp ◦D = TpD , and (27) Ŝp ◦D = T−1p D . (28)
7. Let D ,E ∈ Rp×p be positive semidefinite matrices, where E has no zero diagonal element, and D is of rank r ≤ p. Also, let for any r ≤ p, Jr = {i1, · · · , ir}(1 ≤ i1 < · · · < ir < n) be any ordered r−index set. Then D and E satisfy
E ( Ŝp ◦D ) = ( Ŝp ◦ E ) D ,
if and only if, the following two conditions are satisfied:
(a) The matrix D is diagonal with p − r zero diagonal elements and r positive diagonal elements indexed by the set Jr. That is for any i ∈ Jr : (D )ii > 0 and the rest of elements of D are zero.
(b) For any i, j ∈ Jr and i 6= j we have (E )i,j = 0.
Clearly, if D is positive definite then Jr = Np and hence, both D and E are diagonal.
Proof. . The proof of the properties are as follows.
1. eq. (25) is trivial. For eq. (26) note that AIi;p selects the first i columns of A (zeros out the rest), and similarly, Ii;pA′ selects the first i rows of A (zeros out the rest). Therefore, Ii;pA
′AIi;p is a p×pmatrix that its Leading Principal Submatrix of order i (LPSi) 1 is the same as the LPSi of A′A (and the rest of the elements are zero). Hence, ∑p i=1 Ii;pA
′AIi;p (counting backwards) adds LPSp of A′A (i.e. A′A itself) with LPSp−1 that doubles LPSp−1 part of the result and then adds LPSp−2 that triples the LPSp−2 part of result, the process continues until by the last addition LPS1is added to the result for the pthtimes. This is exactly the same as evaluating Sp ◦ (A′A).
2. This is a standard result (Horn & Johnson, 2012), and no proof is needed.
3. Directly compute SpS−1p :
( SpS −1 p ) ij = p∑ k=1 (Sp)ik(S −1 p )kj ∀|k−j|>1:(S−1p )kj=0 =============⇒
=  (Sp)i,j−1(S−1p )j−1,j+(Sp)i,j(S −1 p )j,j+(Sp)i,j+1(S −1 p )j+1,j 2 ≤ j& j ≤p− 1 (Sp)i,p−1(S−1p )p−1,p + (Sp)i,p(S −1 p )p,p j = p
(Sp)i,1(S −1 p )1,1 + (Sp)i,2(S −1 p )2,1 j = 1
=  −(Sp)i,j−1 + 2(Sp)i,j − (Sp)i,j+1 2 ≤ j ≤ p− 1 −(Sp)i,p−1 + 2(Sp)i,p j = p (Sp)i,1 − (Sp)i,2 j = 1
=  max(i, j − 1)− 2 max(i, j) + max(i, j + 1) 2 ≤ j ≤ p− 1 −(p−max(i, p− 1) + 1) + 2(p−max(i, p) + 1) j = p −max(i, 1) + max(i, 2) j = 1
1For a p × p matrix, the leading principal submatrix of order i is an i × i matrix derived by removing the last p− i rows and columns of the original matrix (Horn & Johnson (2012), P17)
=  max(i, j − 1)− 2 max(i, j) + max(i, j + 1) 2 ≤ j ≤ p− 1 1− p+ max(i, p− 1) j = p max(i, 2)−max(i, 1) j = 1
=  { 1 i = j 0 i 6= j 1 < j < p{ 1 i = p 0 i 6= p j = p{ 1 i = 1
0 i ≥ 2 j = 1
= (Ip)ij .
4. Firstly, note that S−1p is symmetric and nonsingular so all the eigenvalues are real and nonzero. It is also a diagonally dominant matrix (Horn & Johnson (2012), Def 6.1.9) since
∀i ∈ {1, · · · , p} : Ci := |(S−1p )ii| ≥ ∑
j=1,j 6=i |(S−1p )ij | =: Ri,
where the inequality is strict for the first and the last row and it is equal for the rows in the middle. Moreover, by Gersgorin circle theorem (Horn & Johnson (2012), Thm 6.1.1) for every eigenvalue li of S−1p there exists i such that li ∈ [Ci − Ri, Ci + Ri]. Since ∀i : Ci ≥ Ri we have all the eigenvalues are non-negative. They are also nonzero, hence, S−1p is positive definite, which implies Sp is also positive definite.
5. For any matrix A ∈ Rn×p, A′A is positive semidefinite. Also, Sp is positive definite so by Schur product theorem (Horn & Johnson (2012), Thm 7.5.3(a)), Sp ◦ (A′A) is positive semidefinite. Moreover, if all diagonal elements of A′A are positive (i.e. A has no zero column) by the extension of Schur product theorem (Horn & Johnson (2012), Thm 7.5.3(b)) it is positive definite. This can also be easily deduced using the Oppenheim inequality (Horn & Johnson (2012), Thm 7.8.16); that is for positive semidefinite matrices Sp and A′A: det(Sp) ∏ i(A ′A)ii ≤ det(Sp◦(A′A)). Since, Sp is positive definite, det(Sp) > 0
(in fact it is 1 for any p) and if A′A has no zero diagonal then det(Sp ◦ (A′A)) > 0 and therefore, Sp ◦ (A′A) is positive definite.
6. Clearly, the matrix Tp is achieved by setting the off-diagonal elements of Sp to zero. Hence, for any diagonal matrix D ∈ Rp×p: Sp ◦D = Tp ◦D . For the diagonal matrices Hadamard product and matrix product are interchangeable so the latter may also be written as TpD . The same argument applies for the second identity.
7. This property can easily be proved by induction on p and careful bookkeeping of indices.
Lemma 3 (Simultaneous diagonalization by congruence). Let M1,M2 ∈ Rp×p, where M1 is positive definite and M2 is positive semidefinite. Also, let D ,E ∈ Rr×r be positive definite diagonal matrices with r ≤ p. Further, assume there is a C ∈ Rr×p of rank r ≤ p such that
CM1C ′ =D and CM2C ′ =DE .
Then there exists a nonsingular C̄ ∈ Rp×p that its first r rows are the matrix C and C̄M1C̄ ′ =D̄ and
C̄M2C̄ ′ =D̄Ē ,
where, D̄ = D̄ ⊕ Ir−p is a p× p diagonal matrix and Ē = E ⊕E is another p× p diagonal matrix, in which E ∈ Rp−r×p−r is a nonnegative diagonal matrix. Clearly, the rank of M2 is r plus the number of nonzero diagonal elements of E .
Proof. The proof is rather straightforward since this lemma is the direct consequence of Theorem 7.6.4 in Horn & Johnson (2012). The theorem basically states that if M1,M2 ∈ Rp×p is symmetric
and M1 is positive definite then there exists an invertible S ∈ Rp×p such that SM1S′ = Ip and SM2S
′ is a diagonal matrix with the same inertia as M2. Here, we have M2 that is positive semidefinite and C ∈ Rr×p of rank r ≤ p such that(
D −1 2 C ) M1 ( D −1 2 C )′ =Ir and(
D −1 2 C ) M2 ( D −1 2 C )′ =E .
Therefore, since S is of full rank p and D −1 2 C is of rank r ≤ p, there exists p − r rows in S that are linearly independent of rows of D −1 2 C. Establish C̄ ∈ Rp×p by adding those p− r rows to C. Then C̄ has p linearly independent rows so it is nonsingular, and fulfills the lemma’s proposition that is
C̄M1C̄ ′ =D̄ and C̄M2C̄ ′ =D̄Ē ,
where, D̄ = D̄ ⊕ Ir−p is a p× p diagonal matrix and Ē = E ⊕E is another p× p diagonal matrix, in which E ∈ Rp−r×p−r is a nonnegative diagonal matrix.
Lemma 4. Let A and B define a critical point of L. Further, let V ∈ Rn×p and W ∈ Rp×n are such that ‖V ‖F , ‖W ‖F = O(ε) for some ε > 0. Then
L(A + V ,B + W )− L(A,B) =〈V TpBΣxxB′,V 〉F −2〈ΣyxW ′Tp −A (Sp ◦ (BΣxxW ′ + WΣxxB′)) ,V 〉F +〈(Sp ◦ (A′A))WΣxx,W 〉F +O(ε3). (29)
Further, for W = W̄ := (Sp ◦ (A′A))−1 TpV ′ΣyxΣ−1xx , the above equation becomes
L(A + V ,B + W̄ )− L(A,B) = Tr (V ′V TpBΣxxB′)− Tr ( V ′ΣV Tp (Sp ◦ (A′A))−1 Tp ) +2 Tr ( V ′A ( Sp ◦ ( BΣxyV Tp (Sp ◦ (A′A))−1
+ (Sp ◦ (A′A))−1 TpV ′ΣyxB′ ))) +O(ε3). (30)
Finally, in case the critical A is of full rank p and so, (A,B) = (UIpΠD, B̂(UIpΠD)), for the encoder direction V with ‖V ‖F = O(ε) and W = W̄ we have,
L(A + V ,B + W )− L(A,B) = Tr ( V ′V Π′ΛIpΠTpD −2)− Tr (V ′ΣV TpD−2) +2 Tr ( V ′UIpΠD ( Sp ◦ ( D−1Π′U ′IpΣV D −2 )))
+2 Tr ( V ′UIpΠD ( Sp ◦ ( D−2V ′ΣUIpΠD −1))) +O(ε3). (31)
Proof. As described in appendix B.1, the second order Taylor expansion for the loss L(A,B) is then given by eq. (63), i.e.
L(A + V ,B + W )− L(A,B) =dAL(A,B)V + dBL(A,B)W + 1
2 d2AL(A,B)V 2
+dABL(A,B)V W + 1
2 d2BL(A,B)W 2 +RV ,W (A,B).
If ‖V ‖F , ‖W ‖F = O(ε) then ‖R(V ,W )‖ = O(ε3). Moreover, when A and B define a critical point of L we have dAL(A,B)V = dBL(A,B)W = 0. By setting the derivatives d2AL(A,B)V 2, dABL(A,B)V W , d2BL(A,B)W 2 that are given by eq. (69), eq. (68), and eq. (66) respectively, the above equation simplifies to
L(A + V ,B + W )− L(A,B) =〈V (Sp ◦ (BΣxxB′)) ,V 〉F −2〈ΣyxW ′Tp −A (Sp ◦ (BΣxxW ′ + WΣxxB′)) ,V 〉F +〈(Sp ◦ (A′A))WΣxx,W 〉F +O(ε3).
Now, based on the first item in Corollary 1, BΣxxB′ is a p×p diagonal matrix, so based on eq. (27): Sp◦(BΣxxB′) = TpBΣxxB′. The substitution then yields eq. (29). Finally, in the above equation replace W with W̄ = (Sp ◦ (A′A))−1 TpV ′ΣyxΣ−1xx . We have
L(A + V ,B + W̄ )− L(A,B) = = 〈V TpBΣxxB′,V 〉F − 2〈ΣyxΣ−1xxΣxyV Tp (Sp ◦ (A′A)) −1 Tp,V 〉F
+2〈A ( Sp ◦ ( BΣxxΣ −1 xxΣxyV Tp(Sp◦(A′A)) −1 +(Sp◦(A′A))−1TpV ′ΣyxΣ−1xxΣxxB′ )) ,V〉F
+〈(Sp ◦ (A′A)) (Sp ◦ (A′A))−1TpV ′ΣyxΣ−1xxΣxx, (Sp ◦ (A′A)) −1 TpV ′ΣyxΣ −1 xx 〉F +O(ε3)
= Tr (V ′V TpBΣxxB ′)− Tr ( V ′ΣV Tp (Sp ◦ (A′A))−1 Tp ) + 2 Tr ( V ′A ( Sp ◦ ( BΣxyV Tp (Sp ◦ (A′A))−1 + (Sp ◦ (A′A))−1 TpV ′ΣyxB′ ))) +O(ε3),
which is eq. (30). For the final equation, we have
TpBΣxxB ′ =TpD −1Π′U ′Ip ΣyxΣ −1 xxΣxxΣ −1 xxΣxy︸ ︷︷ ︸UIpΠD−1
=TpD −1Π′U ′IpΣUIp︸ ︷︷ ︸ΠD−1 = TpD−1 Π′ΛIpΠ︸ ︷︷ ︸D−1 =Π′ΛIpΠTpD −2, and (32)
Tp (Sp ◦ (A′A))−1 Tp =Tp ( Sp ◦ ( DΠ′U ′IpUIpΠ︸ ︷︷ ︸D ))−1 Tp
=Tp ( Sp ◦D2 )−1 Tp = TpT −1 p D −2Tp = TpD −2. (33)
Replace the above in eq. (30) and simplify: L(A + V ,B + W )− L(A,B) = Tr (V ′V TpBΣxxB′)− Tr ( V ′ΣV Tp (Sp ◦ (A′A))−1 Tp ) +2 Tr ( V ′A ( Sp ◦ ( BΣxyV Tp (Sp ◦ (A′A))−1
+ (Sp ◦ (A′A))−1 TpV ′ΣyxB′ )))
+O(ε3) eq. (32) ====⇒ eq. (33)
L(A + V ,B + W )− L(A,B) = Tr ( V ′V Π′ΛIpΠTpD −2)− Tr (V ′ΣV TpD−2) +2 Tr ( V ′A ( Sp ◦ ( BΣxyV D −2 + D−2V ′ΣyxB ′)))
+O(ε3) A=UIpΠD
=========⇒ B=B̂(UIpΠD)
L(A + V ,B + W )− L(A,B) = Tr ( V ′V Π′ΛIpΠTpD −2)− Tr (V ′ΣV TpD−2) +2 Tr ( V ′UIpΠD ( Sp ◦ ( D−1Π′U ′IpΣV D −2 )))
+2 Tr ( V ′UIpΠD ( Sp ◦ ( D−2V ′ΣUIpΠD −1))) +O(ε3),
which finalizes the proof.

heading: A.2 PROOF OF PROPOSITION 1
text: For this proof we use the first and second order derivatives for L(A,B) wrt B derived in Lemma 5. From eq. (66), we have that for a given A the second derivative wrt to B of the cost L(A,B) at B,
and in the direction W is the quadratic form
d2B2L(A,B)W 2 = 2 Tr (W ′ (Sp ◦A′A)WΣxx) .
The matrix Σxx is positive-definite and by Lemma 2, Sp ◦ A′A is positive-semidefinite. Hence, d2B2L(A,B)W
2 is clearly non-negative for all W ∈ Rp×n. Therefore, L(A,B) is convex in coefficients of B for a fixed matrix A. Also the critical points of L(A,B) for a fixed A is a matrix B that satisfies ∀W ∈ Rp×n : dBL(A,B)W = 0 and hence, from eq. (64) we have
−2〈TpA′Σyx − (Sp ◦ (A′A))BΣxx,W 〉F = 0. Setting W = TpA′Σyx − (Sp ◦ (A′A))BΣxx we have
TpA ′Σyx − (Sp ◦ (A′A))BΣxx = 0.
For a fixed A, the cost L(A,B) is convex in B, so any matrix B that satisfies the above equation corresponds to a minimum of L(A,B). Further, if A has no zero column then by Lemma 2, Sp ◦ A′A is positive definite. Hence, ∀W ∈ Rp×n : d2B2L(A,B)W 2 = 2 Tr (W ′ (Sp ◦A′A)WΣxx) is positive. Therefore, the cost L(A,B) becomes strictly convex and the unique global minimum is achieved at B = B̂(A) as defined in eq. (6).

heading: A.3 PROOF OF PROPOSITION 2
text: For this proof we use the first and second order derivatives for L(A,B) wrt A derived in Lemma 6. For a fixed B, based on eq. (69) the second derivative wrt to A of L(A,B) at A, and in the direction V is the quadratic form
d2A2L(A,B)V 2 = 2〈V (Sp ◦ (BΣxxB′)) ,V 〉F = 2 Tr (V (Sp ◦ (BΣxxB′))V ′) .
The matrix Σxx is positive-definite and by Lemma 2, Sp ◦ (BΣxxB′) is positive-semidefinite. Hence, d2A2L1(A,B)V
2is non-negative for all V ∈ Rn×p. Therefore, L(A,B) is convex in coefficients of A for a fixed matrix B. Based on eq. (67) the critical point of L(A,B) for a fixed B is a matrix A that satisfies for all V ∈ Rn×p
dAL(A,B)V = 〈−2 (ΣyxB′Tp −A (Sp ◦ (BΣxxB′))) ,V 〉F = 0 =⇒ ΣyxB
′Tp = A (Sp ◦ (BΣxxB′)) , which is eq. (7).

heading: A.4 PROOF OF THEOREM 1
text: Before we start, a reminder on notation and some useful identities that are used throughout the proof. The matrix Σ := ΣyxΣ−1xxΣxy has an eigenvalue decomposition Σ = UΛU ′, where the ith column of U , denoted as ui, is an eigenvector of Σ corresponding to the ith largest eigenvalue of Σ, denoted as λi. Also, Λ = diag(λ1, · · · , λn) is the diagonal vector of ordered eigenvalues of Σ, with λ1 > λ2 > · · · > λn > 0. We use the following notation to organize a subset of eigenvectors of Σ into a rectangular matrix. Let for any r ≤ p, Ir = {i1, · · · , ir}(1 ≤ i1 < · · · < ir < n) be any ordered r−index set. Define UIr ∈ Rn×p as UIr = [ui1 , · · · ,uir ]. That is the columns of UIr are the ordered orthonormal eigenvectors of Σ associated with eigenvalues λi1 < · · · < λir . The following identities are then easy to verify:
U ′IrUIr =Ir,
ΣUIr =UIrΛIr , (34) U ′IrΣUIr =ΛIr . (35)

heading: The sufficient condition:
text: Let A ∈ Rn×pof rank r ≤ p and no zero column be given by eq. (8), B ∈ Rp×n given by eq. (9), and the accompanying conditions are met. Notice that U ′IrUIr = Ir implies that DC
′CD = DC ′U ′IrUIrCD = A ′A, so
B = D−1ΠCU ′ IrΣyxΣ −1 xx
ΠC :=(Sp◦(C′C))−1TpC′ =================⇒
D−1D=Ip
B = D−1 (Sp ◦ (C ′C))−1 D−1DTpC ′U ′IrΣyxΣ−1xx Lemma 2-2
=======⇒ DTp=TpD
B = ( Sp ◦ (DC ′CD)︸ ︷︷ ︸ )−1 TpDC
′U ′Ir︸ ︷︷ ︸ΣyxΣ−1xx A′=D′C′U ′Ir=========⇒DC′CD=A′A B = (Sp ◦ (A′A))−1 TpA′ΣyxΣ−1xx = B̂(A),
which is eq. (6). Therefore, based on Proposition 1, for the given A, the matrix B defines a critical point of L(A,B). For the gradient wrt to A, first note that with B given by eq. (9) we have
BΣxxB ′ =D−1ΠCU ′ IrΣyxΣ −1 xxΣxxΣ −1 xxΣxyUIrΠ ′ CD −1
=D−1ΠC U ′ IrΣyxΣ −1 xxΣxyUIr︸ ︷︷ ︸Π′CD−1 eq. (35)===⇒
BΣxxB ′ =D−1ΠCΛIrΠ ′ CD −1. (36)
The matrix ΠC is a rectangular permutation matrix so ΠCΛIrΠ ′ C is diagonal so as D−1ΠCΛIrΠ ′ CD −1. Therefore, BΣxxB′ is diagonal and by eq. (27) in Lemma 2-6 we have
Sp ◦ (BΣxxB′) =TpBΣxxB′ = BΣxxB′Tp =D−1ΠCΛIrΠ ′ CD −1Tp A× ==⇒
A (Sp ◦ (BΣxxB′)) =AD−1ΠCΛIrΠ′CD−1Tp A=UIrCD=======⇒
A (Sp ◦ (BΣxxB′)) =UIrCDD−1ΠCΛIrΠ′CD−1Tp A=UIrCD=======⇒ =UIr CΠC︸ ︷︷ ︸ΛIrΠ′CD−1Tp CΠC=Ir======⇒ A (Sp ◦ (BΣxxB′)) =UIrΛIr︸ ︷︷ ︸Π′CD−1Tp eq. (34)===⇒
=ΣUIrΠ ′ CD −1Tp =ΣyxΣ −1 xxΣxyUIrΠ ′ CD −1Tp
=Σyx ( D−1ΠCU ′ IrΣyxΣ −1 xx )′︸ ︷︷ ︸Tp =ΣyxB ′Tp,
which is eq. (7). Therefore, based on Proposition Proposition 2, for the given B, the matrix A define a critical point of L(A,B). Hence, A and B together define a critical point of L(A,B).

heading: The necessary condition:
text: Based on Proposition 1 and Proposition 2, for A (with no zero column) and B, to define a critical point of L(A,B), B has to be B̂(A) given by eq. (6), and A has to satisfy eq. (7). That is
A ( Sp ◦ ( B̂ΣxxB̂ ′ )) =ΣyxB̂ ′Tp B̂(A) on RHS =======⇒
A ( Sp ◦ ( B̂ΣxxB̂ ′ )) =ΣxyΣ −1 xxΣyxATp(Sp ◦ (A′A))−1Tp
×A′ ==========⇒ Σ=ΣxyΣ −1 xx Σyx
A ( Sp ◦ ( B̂ΣxxB̂ ′ )) A′ =ΣATp(Sp ◦ (A′A))−1TpA′ Σ=UΛU ′′
======⇒ ×U ,U ′×
U ′A ( Sp ◦ ( B̂ΣxxB̂ ′ )) A′U =U ′UΛU ′ATp(Sp ◦ (A′A))−1TpA′U U ′U=In=====⇒
U ′A ( Sp ◦ ( B̂ΣxxB̂ ′ )) A′U =Λ∆, (37)
where, ∆ := U ′ATp(Sp ◦ (A′A))−1TpA′U is symmetric and positive semidefinite. The LHS of the above equation is symmetric so the RHS is symmetric too, so Λ∆ = (Λ∆)′ = ∆′Λ′ = ∆Λ. Therefore, ∆ commutes with the diagonal matrix of eigenvalues Λ. Since, eigenvalues are assumed to be distinct, ∆ has to be diagonal as well. By Lemma 2 Tp(Sp ◦ (A′A))−1Tp is positive definite and U is an orthogonal matrix. Therefore, r = rank(A) = rank(∆) = rank(U ′∆U), which implies that the diagonal matrix ∆, has r nonzero and positive diagonal entries. There exists an
r−index set Ir corresponding to the nonzero diagonal elements of ∆. Forming a diagonal matrix ∆Ir ∈ Rr×r by filling its diagonal entries (in order) by the nonzero diagonal elements of ∆ we have
U∆U ′ = UIr∆IrU ′ Ir Def of ∆ ====⇒
UU ′ATp(Sp ◦ (A′A))−1TpA′UU ′ = UIr∆IrU ′Ir UU ′=In=====⇒
ATp(Sp ◦ (A′A))−1TpA′ = UIr∆IrU ′Ir , (38) which indicates that the matrix A has the same column space as UIr . Therefore, there exists a full rank matrix C̃ ∈ Rr×p such that A = UIrC̃. Since A has no zero column, C̃ has no zero column. Further, by normalizing the columns of C̃ we can write A = UIrCD, where D ∈ Rp×p is diagonal that contains the norms of columns of C̃. Therefore, A is exactly in the form given by eq. (8). The matrix C has to satisfy eq. (38) that is
ATp(Sp ◦ (A′A))−1TpA′ = UIr∆IrU ′Ir A=UIrC=====⇒
UIrCDTp(Sp ◦ (A′A))−1TpDC ′U ′Ir = UIr∆IrU ′Ir ×UIr ,UIr×=========⇒
A′A=DC′CD
CDTp(Sp ◦ (DC ′CD))−1TpC ′D = ∆Ir Lemma 2-2 =====⇒
CTpDD −1(Sp ◦ (C ′C))−1D−1DTpC ′ = ∆Ir =⇒
CTp(Sp ◦ (C ′C))−1TpC ′ = ∆Ir . (39) Now that the structure of A has been identified, evaluate B̂(A) of eq. (6) by setting A = UIrCD, that is
B =B̂(A) = (Sp ◦ (A′A))−1TpA′ΣyxΣ−1xx =(Sp ◦ (DC ′CD))−1TpDC ′U ′IrΣyxΣ−1xx Lemma 2-2 =====⇒
B =D−1(Sp ◦ (C ′C))−1TpC ′U ′IrΣyxΣ−1xx , which by defining ΠC := (Sp ◦ (C ′C))−1 TpC ′ gives eq. (34) for B as claimed. While C has to satisfy eq. (39), A and B in the given form have to satisfy eq. (37) that provides another condition for C as follows. First, note that
Sp ◦ ( B̂ΣxxB̂ ′ ) = Sp ◦ ( D−1(Sp ◦ (C ′C))−1TpC ′U ′IrΣUIrCTp(Sp ◦ (C ′C))−1D−1 ) = Sp ◦ ( D−1(Sp ◦ (C ′C))−1TpC ′ΛIrCTp(Sp ◦ (C ′C))−1D−1 ) Lemma 2-2 =====⇒
= D−1 ( Sp ◦ ( (Sp ◦ (C ′C))−1TpC ′ΛIrCTp(Sp ◦ (C ′C))−1 )) D−1
Now, replace A and B in eq. (37) by their respective identities that we just derived. Performing the same process for eq. (37) we have
U ′A ( Sp ◦ ( B̂ΣxxB̂ ′ ))
A′U = Λ∆ A=UIrCD=======⇒ ×U ′,U×
UIrC ( Sp ◦ ( (Sp ◦ (C ′C))−1TpC ′ΛIrCTp(Sp ◦ (C ′C))−1 )) C ′U ′Ir = UΛ∆U
′ ×UIr===⇒ U ′Ir×
C ( Sp ◦ ( (Sp ◦ (C ′C))−1TpC ′ΛIrCTp(Sp ◦ (C ′C))−1 )) C ′ = U ′IrUΛ∆U
′UIr =⇒ C ( Sp ◦ ( (Sp ◦ (C ′C))−1TpC ′ΛIrCTp(Sp ◦ (C ′C))−1 )) C ′ = ΛIr∆Ir . (40)
Now we have to find C such that it satisfies eq. (39) and eq. (40). To make the process easier to follow, lets have them in one place. The matrix C ∈ Rr×p have to satisfy
CTp (Sp ◦ (C ′C))−1 TpC ′ =∆Ir and (41) C ( Sp ◦ ( (Sp ◦ (C ′C))−1TpC ′ΛIrCTp(Sp ◦ (C ′C))−1 )) C ′ =ΛIr∆Ir . (42)
Since C is a rectangular matrix, solving above equations for C in this form seems intractable. We use a trick to temporarily extend C into an invertible square matrix as follows.
• Temporarily, let M1 = Tp (Sp ◦ (C ′C))−1 Tp, and M2 = Sp ◦( (Sp ◦ (C ′C))−1TpC ′ΛIrCTp(Sp ◦ (C ′C))−1 ) . Then M1 is positive definite and
M2 is positive semidefinite, so they are simultaneously diagonalizable by congruence that is based on Lemma 3 and eq. (41) and eq. (42), there exists a nonsingular C̄ ∈ Rp×p such that C consists of the first r rows of C̄ and
C̄Tp (Sp ◦ (C ′C))−1 TpC̄ ′ =∆̄Ir , (43) C̄ ( Sp ◦ ( (Sp ◦ (C ′C))−1 TpC ′ΛIrCTp (Sp ◦ (C ′C)) −1 )) C̄ ′ =Λ̄Ir∆̄Ir , (44)
where, ∆̄Ir = ∆Ir ⊕ Ir−p is a p× p diagonal matrix and Λ̄Ir = ΛIr ⊕Λ is another p× p diagonal matrix, in which Λ ∈ Rr−p×r−p is a nonnegative diagonal matrix. • Substitute ∆̄Ir from eq. (43) in eq. (44), then left multiply by C̄ ′−1, and right multiply by C̄ ′Ir;p:
C̄ ( Sp ◦ ( (Sp ◦ (C ′C))−1 TpC ′ΛIrCTp (Sp ◦ (C ′C)) −1 )) C̄ ′ =
Λ̄IrC̄Tp (Sp ◦ (C ′C)) −1 TpC̄ ′ C̄ ′Ir;p× =====⇒ ×C̄′−1
C̄ ′Ir;pC̄ ( Sp ◦ ( (Sp ◦ (C ′C))−1 TpC ′ΛIrCTp (Sp ◦ (C ′C)) −1 )) =
C̄ ′Ir;pΛ̄IrC̄Tp (Sp ◦ (C ′C)) −1 Tp.
• Now we can revert back everything to C again. Since C consists of the first r rows of C̄ we have C̄ ′Ir;pC̄ = C ′C, and C̄ ′Ir;pΛ̄IrC̄ = C
′ΛIrC, which turns the above equation into
C ′C ( Sp ◦ ( Ip (Sp ◦ (C ′C))−1 TpC ′ΛIrCTp (Sp ◦ (C ′C)) −1 Ip )) =
IpC ′ΛIrCTp (Sp ◦ (C ′C)) −1 Tp.
• In the above equation, replace Ip by T−1p Tp in LHS and by T−1p (Sp ◦ (C ′C))T−1p Tp (Sp ◦ (C ′C))−1 Tp in the RHS. Use ΠC := (Sp ◦ (C ′C))−1 TpC ′ to shrink it into : C ′C ( Sp ◦ ( T−1p TpΠCΛIrΠ ′ CTpT −1 p )) =T−1p (Sp ◦ (C ′C))T−1p TpΠCΛIrΠ′CTp.
• By the second property of Lemma 2 we can collect diagonal matrices T−1p ’s around Sp to arrive at
(C ′C) ( Ŝp ◦ (TpΠCΛIrΠ′CTp) ) = ( Ŝp ◦ (C ′C) ) (TpΠCΛIrΠ ′ CTp) ,
where, Ŝp := T−1p SpT −1 p .
• Define p× p matrices E r := C ′C and D r := TpΠCΛIrΠ′CTp. Substitute in the above to arrive at: E r ( Ŝp ◦D r ) = ( Ŝp ◦ E r ) D r.
Both D r and E r in the above identity are positive semidefinite. Moreover, since by assumption C has no zero columns, E r has no zero diagonal element. Then the 7th property of Lemma 2 implies the following two conclusions:
1. The matrix D r is diagonal. The rank of D r is r so it has exactly r positive diagonal elements and the rest is zero. This argument is true for T−1p D rT −1 p =
ΠCΛIrΠ ′ C . Since ΛIr is a diagonal positive definite matrix, the p × r matrix ΠC := (Sp ◦ (C ′C))−1 TpC ′ of rank r should have p − r zero rows. Let Jr be an r−index set corresponding to nonzero diagonal elements of ΠCΛIrΠ′C . Then the matrix ΠC [Jr,Nr] (r × r submatrix of ΠC consist of its Jr rows) is nonsingular. 2. For every i, j ∈ Jr and i 6= j, (E r)i,j = 0. Since E r := C ′C and so (E r)i,j is the inner product of ith and jth columns of C, we conclude that the columns of C[Nr, Jr] (r × r submatrix of C consist of its Jr columns) are orthogonal or in other words C[Nr, Jr]′C[Nr, Jr] is diagonal. The columns of C are normalized. Therefore, C[Nr, Jr]′C[Nr, Jr] = Ir and hence, C[Nr, Jr] is an orthogonal matrix.
• We use the two conclusions to solve the original eq. (41) and eq. (42). First use ΠC := (Sp ◦ (C ′C))−1 TpC ′ to shrink them into :
CTpΠC =∆Ir , (45)
C (Sp ◦ (ΠCΛIrΠ′C))C ′ =ΛIr∆Ir . (46) Next, by the first conclusion, the matrix T−1p D rT −1 p = ΠCΛIrΠ ′ C is diagonal and so eq. (46) becomes
CTpΠC︸ ︷︷ ︸ΛIrΠ′CC ′ =ΛIr∆Ir eq. (45)===⇒ ∆IrΛIrΠ ′ CC
′ =ΛIr∆Ir =⇒ Π′CC
′ = CΠC =Ir, (47) which is one of the two claimed conditions. What is left is to show that ΠC is a rectangular permutation matrix. From the first conclusion we also have ΠC has exactly r nonzero columns indexed by Jr so
C[Nr, Jr]ΠC [Jr,Nr] =Ir. By the second conclusion C[Nr, Jr] is an orthogonal matrix therefore, ΠC [Jr,Nr] is the orthogonal matrix C[Nr, Jr]′. Moreover, we had T−1p D rT−1p = ΠCΛIrΠ′C is a p × p diagonal matrix with exactly r nonzero diagonal elements. Hence, ΠC [Nr, Jr]ΛIrΠ′C [Nr, Jr] is an r × r positive definite diagonal matrix with ΛIr having distinct diagonal elements, and ΠC [Nr, Jr] being orthogonal. Therefore, ΠC [Jr,Nr] (as well as C[Nr, Jr]) should be a square permutation matrix. Putting back the zero columns, we conclude that C should be such that ΠC := (Sp ◦ (C ′C))−1 TpC ′ is a rectangular permutation matrix and CΠC = Ir. Note that it is possible to further analyze these conditions and determine the exact structure of C. However, this is not needed in general for the critical point analysis of the next theorem except for the case where r = p and C is a square invertible matrix. In this case, square matrix ΠC is of full rank p, Jr = Np and therefore, C[Nr, Jr] = C[Np,Np] = C. Hence, C is any square permutation matrix Π, C ′C = Π′Π = Ip and ΠC := (Sp ◦ (C ′C))−1 TpC ′ = T−1p TpΠ′ = Π′, which verifies eq. (10) and eq. (11) for A and B when A is of full rank p.

heading: A.5 PROOF OF COROLLARY 1
text: 1. We already show in the proof Theorem 1 that for critical (A,B) the matrix BΣxxB′ is given by eq. (36) that is
BΣxxB ′ =D−1ΠCΛIrΠ ′ CD −1.
The matrix ΠC is a p × r rectangular permutation matrix so ΠCΛIrΠ′C is diagonal as well as D−1ΠCΛIrΠ ′ CD
−1. Therefore, BΣxxB′ is diagonal. The diagonal matrix ΛIr is of rank r therefore, BΣxxB′ is of rank r.
2. Again by Theorem 1 critical (A,B) is of the form given by eq. (8) and eq. (9) with the proceeding conditions on the invariance C. Therefore, the global map is
G = AB = UIrCDD −1ΠCU ′ IrΣyxΣ −1 xx
= UIrCΠCU ′ IrΣyxΣ −1 xx CΠC=Ir======⇒ G = UIrU ′ IrΣyxΣ −1 xx .
3. Based on Baldi & Hornik (1989) (A,B) define a critical point of L̃(A,B) =∑p i=1 ‖Y −ABX‖ 2 F iff they satisfy
A′ABΣxx =A ′Σyx and (48) ABΣxxB ′ =ΣyxB
′. (49) Again by assumption (A,B) define a critical point of L(A,B) so by Theorem 1 they are of the form given by eq. (8) and eq. (9) with the proceeding conditions on the invariance C. Hence,
A′ABΣxx =DC ′U ′IrUIr︸ ︷︷ ︸CDD−1︸ ︷︷ ︸ΠCU ′IrΣyx Σ−1xxΣxx︸ ︷︷ ︸
=DC ′CΠC︸ ︷︷ ︸U ′IrΣyx CΠC=Ir======⇒ A′ABΣxx =DC ′U ′IrΣyx = A ′Σyx.
Hence, eq. (48) is satisfied. For the second equation we use the first property of this corollary that is BΣxxB′ is diagonal and satisfy eq. (7) of Proposition 2 that is
A (Sp ◦ (BΣxxB′)) =ΣyxB′Tp BΣxxB ′ is diagonal ===========⇒
ATpBΣxxB ′ =ΣyxB ′Tp BΣxxB ′ is diagonal ===========⇒ ABΣxxB ′Tp =ΣyxB
′Tp =⇒ ABΣxxB ′ =ΣyxB ′.
Hence, the second condition, eq. (49) is also satisfied. Therefore, any critical point of L(A,B) is a critical point of L̃(A,B).

heading: A.6 PROOF OF LEMMA 1
text: Proof. We have
L(A,B) = p∑ i=1 ‖Y −AIi;pBX‖2F = p∑ i=1 〈Y −AIi;pBX,Y −AIi;pBX〉F
= p∑ i=1 (〈Y ,Y 〉F + 〈Y ,−AIi;pBX〉F + 〈−AIi;pBX,Y 〉F +〈−AIi;pBX,−AIi;pBX〉F )
= p〈Y ,Y 〉F − 2〈Y ,A (
p∑ i=1 Ii;p
) BX〉F +
p∑ i=1 〈AIi;pBX,AIi;pBX〉F eq. (25) ===⇒
= pTr(Y Y ′)− 2 Tr (ATpBXY ′) + p∑ i=1 Tr (X ′B′Ii;pA ′AIi;pBX)
= pTr(Σyy)− 2 Tr (ATpBΣxy) + Tr ( XX ′B′
p∑ i=1 (Ii;pA ′AIi;p)B
) eq. (26) ===⇒
= pTr(Σyy)− 2 Tr (ATpBΣxy) + Tr (B′ (Sp ◦ (A′A))BΣxx) , which is eq. (17).

heading: A.7 PROOF OF THEOREM 2
text: Proof. The full rank matrices A∗ and B∗ given by eq. (18) and eq. (19) are clearly of the form given by Theorem 1 with Ip = Np := {1, 2, · · · , p}, and Πp = Ip. Hence, they define a critical point of L(A,B). We want to show that these are the only local minima, that is any other critical (A,B) is a saddle points. The proof is similar to the second partial derivative test. However, in this case the Hessian is a forth order tensor. Therefore, the second order Taylor approximation of the loss, derived in Lemma 4, is used directly. To prove the necessary condition, we show that at any other critical point (A,B), where the first order derivatives are zero, there exists infinitesimal direction along which the second derivative of loss is negative. Next, for the sufficient condition we show that the any critical point of the form (A∗,B∗) is a local and global minima.

heading: The necessary condition:
text: Recall that UIp is the matrix of eigenvectors indexed by the p−index set Ip and Π is a p × p permutation matrix. Since all the index sets Ir, r ≤ p are assumed to be ordered, the only way to have UNp = UIpΠ is by having Ip = Np and Π = Ip. Let A (with no zero column) and B define an arbitrary critical point of L(A,B). Then Based on the previous theorem, either A = UIrC with r < p or A = UIpΠD while in both cases B = B̂(A) given by eq. (6). If (A,B) is not of the form of (A∗,B∗) then there are three possibilities either 1) A = UIrCD with r < p, or 2)
A = UIpΠD with Ip 6= Np or 2) A = UNpΠD but Π 6= Ip. The first two cases corresponds to not having the “right” and/or “enough” eigenvectors, and the third corresponds to not having the “right” ordering. We introduce the following notation and investigate each case separately. Let ε > 0 and Ui;j ∈ Rn×p be a matrix of all zeros except the ith column, which contains uj ; the eigenvector of Σ corresponding to the jth largest eigenvalue. Therefore,
U ′i;jΣUi;j = U ′ i;jUΛU ′Ui;j = λjEi, (50)
where, Ei ∈ Rp×p is matrix of zeros except the ith diagonal element that contains 1. In what follows, for each case we define a encoder direction V ∈ Rn×p with ‖V ‖F = O(ε), and set the decoder direction W ∈ Rp×n as W = W̄ := (Sp ◦ (A′A))−1 TpV ′ΣyxΣ−1xx . Then we use eq. (30) and eq. (31) of Lemma 4, to show that the given direction (V ,W ) infinitesimally reduces the loss and hence, in every case the corresponding critical (A,B) is a saddle point.
1. For the case A = UIrCD, with r < p, note that based on the first item in Corollary 1, BΣxxB
′ is a p×p diagonal matrix of rank r so it has p−r zero diagonal elements. Pick an i ∈ Np such that (BΣxxB′)ii is zero and a j ∈ Np \ Ir. Set V = εUi;jD and W = W̄ . Clearly,
V ′A =εDU ′i;jUIrCD = 0, (51)
V ′V TpBΣxxB ′ =ε2DU ′i;jUi;j︸ ︷︷ ︸DTpBΣxxB′,
=ε2DEiDTpBΣxxB ′ = ε2D2TpEi (BΣxxB ′) = 0 and (52)
V ′ΣV =ε2DU ′i;jUΛU ′Ui;jD = ε 2λjD 2Ei. (53)
Notice, ‖V ‖F , ‖W ‖F = O(ε), so based on eq. (30) of Lemma 4, we have L(A + V ,B + W )− L(A,B) = Tr (V ′V TpBΣxxB ′)− Tr ( V ′ΣV Tp (Sp ◦ (A′A))−1 Tp
) + 2 Tr ( V ′A ( Sp ◦ ( BΣxyV Tp (Sp ◦ (A′A))−1 + (Sp ◦ (A′A))−1 TpV ′ΣyxB′
))) +O(ε3)
eq. (51) ====⇒ eq. (52)
L(A + V ,B + W )− L(A,B) = − Tr ( V ′ΣV Tp (Sp ◦ (A′A))−1 Tp ) +O(ε3)
eq. (53) =========⇒ A′A=DC′CD
L(A + V ,B + W )− L(A,B) = − ε2λj Tr ( D2EiD −1 (( T−1p SpT −1 p︸ ︷︷ ︸ ) ◦ (C ′C) )−1 D−1 ) +O(ε3) =
− ε2λj (( Ŝp ◦ (C ′C) )−1)
ii
+O(ε3).
Therefore, since ( Ŝp ◦ (C ′C) )−1 is a positive definite matrix, as ε→ 0, we have L(A +
V ,B + W ) ≤ L(A,B). Hence, any (A,B) = (UIrCD, B̂(UIrCD)) with r < p is a saddle point.
2. Next, consider the case where A = UIpΠD with Ip 6= Np. Then there exists at least one j ∈ Ip \ Np and i ∈ Np \ Ip such that i < j (so λi > λj). Let σ be the permutation corresponding to the permutation matrix Π. Also, let ε > 0 and Uσ(j);i ∈ Rn×p be a matrix of all zeros except the σ(j)th column, which contains ui; the eigenvector of Σ corresponding to the ith largest eigenvalue. Set V = εUσ(j);iD and W = W̄ . Then, since i /∈ Ip we have
V ′UIp =εDU ′ σ(j);iUIp = 0, (54)
V ′V =ε2DU ′σ(j);iUσ(j);iD = ε 2D2Eσ(j), and (55)
V ′ΣV =ε2DU ′σ(j);iUΛU ′Uσ(j);iD = ε 2λiD 2Eσ(j). (56)
Since ‖V ‖F , ‖W ‖F = O(ε), based on eq. (31) of Lemma 4, we have L(A + V ,B + W )− L(A,B) = Tr ( V ′V Π′ΛIpΠTpD −2)− Tr (V ′ΣV TpD−2) +2 Tr ( V ′UIpΠD ( Sp ◦ ( D−1Π′U ′IpΣV D −2 )))
+2 Tr ( V ′UIpΠD ( Sp ◦ ( D−2V ′ΣUIpΠD −1))) +O(ε3)
eq. (54) =========⇒ eq. (55),eq. (56)
L(A + V ,B + W )− L(A,B) = Tr ( ε2D2Eσ(j)Π ′ΛIpΠTpD −2)
−Tr ( ε2λiD 2Eσ(j)TpD −2)+O(ε3)
=ε2 Tr ( Eσ(j)Π ′ΛIpΠ︸ ︷︷ ︸Tp ) −ε2λi Tr ( Eσ(j)Tp ) +O(ε3)
=ε2λj Tr ( Eσ(j)Tp ) − ε2λi Tr ( Eσ(j)Tp ) +O(ε3)
=− ε2(p− σ(j) + 1) (λi − λj) +O(ε3). Note that in the above, the diagonal matrix Π′ΛIpΠ has the same diagonal elements as ΛIp but they are permuted by σ. So Eσ(j)Π ′ΛIpΠ selects σ(j) th diagonal element of Π′ΛIpΠ that is the j thdiagonal element of ΛIp , which is nothing but λj . Now, since i < j so λi > λj and σ(j) ≤ p, as ε→ 0, we have L(A+V ,B +W ) ≤ L(A,B). Hence, any (A,B) = (UIpΠD, B̂(UIpΠD)) is a saddle point.
3. Finally consider the case where A = UNpΠD with Π 6= Ip. Since Π 6= Ip, the permutation σ of the set Np, corresponding to the permutation matrix Π, has at least a cycle (i1i2 · · · ik), where 1 < i1 < i2 · · · < ik < p and 2 ≤ k ≤ p. Hence, Π can be decomposed as Π = Π(i1i2···ik)Π̂, where Π̂ is the permutation matrix corresponding to other cycles of σ. The cycle (i1i2 · · · ik) can be decomposed into transpositions as (i1i2 · · · ik) = (ikik−1) · · · (iki1), which in matrix form is Π(i1i2···ik) = Π(iki1)Π(iki2) · · ·Π(ikik−1). Therefore, Π can be decomposed as Π = Π(iki1)Π̃, where Π̃ = Π(iki2) · · ·Π(ikik−1)Π̂. Note that Π(iki1), the permutation matrix corresponding to transposition (iki1) is a symmetric involutory matrix, i.e. Π2(iki1) = Ip. Set V = ε(Ui1;i1−Uik;ik)Π̃D and W = W̄ . Again we replace V and W in eq. (31) of Lemma 4. There are some tedious steps to simplify the equation, which is given in appendix A.7.1. The final result is as follows. With the given V and W , the third and forth terms of the RHS of eq. (31) are canceled and the first two terms are simplified to
Tr ( V ′V Π′ΛNpΠTpD −2) =ε2λik(p− i1 + 1) + ε2λi1(p− im + 1), and (57) Tr ( V ′ΣV TpD
−2) =ε2λi1(p− i1 + 1) + ε2λik(p− im + 1), (58) in which, m = max{k − 1, 2}. This means that If the selected cycle is just a transposition (i1i2) then im = i2. But if for the selected cycle (i1i2 · · · ik), k is greater than 2 then im = ik−1. Using above equations, eq. (31) yields
L(A+V ,B+W )−L(A,B)=Tr ( V ′V Π′ΛIpΠTpD −2)−Tr (V ′ΣV TpD−2)+O(ε3) =ε2λik(p− i1 + 1) + ε2λi1(p− im + 1) −ε2λi1(p− i1 + 1)− ε2λik(p− im + 1) +O(ε3) =− ε2i1λik − ε2imλi1 + ε2i1λi1 + ε2imλik =− ε2 ((λi1 − λik)(im − i1)) +O(ε3). (59)
By the above definition of im, we have im − i1 > 0 and since i1 < ik, λi1 − λik > 0. Hence, the first term in the above equation is negative and as ε → 0, we have L(A + V ,B +W )−L(A,B) < 0. Therefore, any any (A,B) = (UIpΠD, B̂(UIpΠD)) with Π 6= Ip is a saddle point.

heading: The Sufficient condition:
text: From Lemma 1 we know that the loss L(A,B) can be written in the form of eq. (17). Use this equation to evaluate loss at (A∗,B∗) = ( UNpDp,D −1 p U ′ NpΣyxΣ −1 xx ) as follows
L(A∗,B∗) = pTr(Σyy)− 2 Tr (A∗TpB∗Σxy) + Tr ( B∗ ′ ( Sp ◦ ( A∗ ′ A∗ )) B∗Σxx ) =⇒
L(A∗,B∗) = pTr(Σyy)− 2 Tr ( UNpDpTpD −1 p U ′ Np ΣyxΣ −1 xxΣxy︸ ︷︷ ︸ ) + Tr (( Sp ◦ ( DpU
′ NpUNp︸ ︷︷ ︸Dp
)) D−1p U ′ Np ΣyxΣ −1 xxΣxxΣ −1 xxΣxy︸ ︷︷ ︸UNpD−1p ) =⇒
L(A∗,B∗) = pTr(Σyy)− 2 Tr ( TpDpD
−1 p︸ ︷︷ ︸U ′NpΣUNp︸ ︷︷ ︸ ) + Tr (( Sp ◦ (Ip)︸ ︷︷ ︸ ) DpD −1 p︸ ︷︷ ︸U ′NpΣUNp︸ ︷︷ ︸D−1p Dp︸ ︷︷ ︸ ) =⇒
L(A∗,B∗) = pTr(Σyy)− 2 Tr ( TpΛNp ) + Tr ( TpΛNp ) =⇒
L(A∗,B∗) = pTr(Σyy)− Tr ( TpΛNp ) = p Tr(Σyy)− p∑ i=1 (p− i+ 1)λi,
which is eq. (20), as claimed. Notice that the above value is independent of the diagonal matrix Dp. From the necessary condition we know that any critical point not in the form of (A∗,B∗) is a saddle point. Hence, due to the convexity of the loss at least one (A∗,B∗) is a global minimum but since the value of the loss at (A∗,B∗) is independent of Dp all these critical points yield the same value for the loss. Therefore, any critical point in the form of (A∗,B∗) is a local and global minima.

heading: A.7.1 SUPPLEMENTARY DETAILS OF THE PROOF OF THEOREM 2
text: To verify eq. (57), eq. (58), and eq. (59) in the proof of Theorem 2, we want to replace V and W in eq. (31) of Lemma 4 with V = ε(Ui1;i1 −Uik;ik)Π̃D and W = W̄ and simplify. eq. (31) is
L(A + V ,B + W )− L(A,B) = Tr ( V ′V Π′ΛIpΠTpD −2)− Tr (V ′ΣV TpD−2) +2 Tr ( V ′UIpΠD ( Sp ◦ ( D−1Π′U ′IpΣV D −2 )))
+2 Tr ( V ′UIpΠD ( Sp ◦ ( D−2V ′ΣUIpΠD −1))) +O(ε3).
We investigate each term on the RHS separately. but before note that
EiΠ̃TpΠ̃ ′ = ( Π̃TpΠ̃ ′ ) i,i Ei = (Tp)σ̃−1(i),σ̃−1(i) Ei = (p− σ̃−1(i) + 1)Ei, (60)
where, σ̃ and its function inverse σ̃−1 are permutations corresponding to Π̃ and Π̃′ respectively. Π̃TpΠ̃
′ is a diagonal matrix where diagonal elements of Tp are ordered based on σ̃−1. Moreover, recall that we decomposed the permutation matrix Π in A with a cycle (i1i2 · · · ik) as Π = Π(i1ik) Π(iki2) · · ·Π(ikik−1)Π̂︸ ︷︷ ︸ = Π(i1ik)Π̃, where i1, i2, · · · ik are fixed points of Π̂. Therefore, with σ̃ being the permutation corresponding to Π̃ we have
σ̃(i1) = i1 =⇒ σ̃−1(i1) = i1, and (61) σ̃(ik−1) = im =⇒ σ̃−1(ik) = im, (62)
where, m = max{k − 1, 2}. This means that If the selected cycle is just a transposition (i1i2) then im = i2. But if for the selected cycle (i1i2 · · · ik), k is greater than 2 then im = ik−1. For the first term we have
V ′V =ε2DΠ̃′(U ′i1;i1 −U ′ik;ik)(Ui1;i1 −Uik;ik)Π̃D U ′i1;i1 Uik;ik=0 =========⇒
V ′V =ε2DΠ̃′(U ′i1;i1Ui1;i1 + U ′ ik;ik Uik;ik)Π̃D U ′i1;i1 Ui1;i1=Ei1 ===========⇒ U ′ik;ik Uik;ik=Eik
V ′V =ε2DΠ̃′(Ei1 + Eik)Π̃D Π̃′(Ei1+Eik )Π̃ is diagonal===============⇒ V ′V =ε2Π̃′(Ei1 + Eik)Π̃D 2 =⇒
Tr ( V ′V Π′ΛNpΠTpD −2) = Tr(︷ ︸︸ ︷V ′V D−2Π̃′Π(i1ik)ΛNpΠ(i1ik)Π̃Tp)
= Tr ε2Π̃′(Ei1 + Eik) Π̃D2D−2Π̃′︸ ︷︷ ︸ Ip Π(i1ik)ΛNpΠ(i1ik)Π̃Tp  =ε2 Tr ( (Ei1 + Eik)Π(i1ik)ΛNpΠ(i1ik)Π̃TpΠ̃ ′ )
=ε2 Tr ( λikEi1Π̃TpΠ̃ ′ + λi1EikΠ̃TpΠ̃ ′ ) eq. (60) ===⇒
Tr ( V ′V Π′ΛNpΠTpD −2) =ε2λik(p− σ̃−1(i1) + 1)Ei1 + ε2λi1(p− σ̃−1(ik) + 1)Eik eq. (61)===⇒eq. (62) Tr ( V ′V Π′ΛNpΠTpD
−2) =ε2λik(p− i1 + 1)Ei1 + ε2λi1(p− im + 1)Eik , which is eq. (57) as claimed.
For the second term we have
V ′ΣV =ε2DΠ̃′(U ′i1;i1 −U ′ik;ik)UΛU ′(Ui1;i1 −Uik;ik)Π̃D =ε2DΠ̃′(U ′i1;i1UΛU
′Ui1;i1︸ ︷︷ ︸ λi1Ei1 −U ′i1;i1UΛU ′Uik;ik︸ ︷︷ ︸ 0
−U ′ik;ikUΛU ′Ui1;i1︸ ︷︷ ︸ 0 +U ′ik;ikUΛU ′Uik;ik︸ ︷︷ ︸
λikEik
)Π̃D
=ε2Π̃′(λi1Ei1 + λikEik)Π̃D 2 =⇒ Tr ( V ′ΣV TpD −2) = Tr(ε2Π̃′(λi1Ei1 + λikEik)Π̃D2TpD−2) =ε2 Tr ( λi1Ei1Π̃TpΠ̃ ′ + λikEikΠ̃TpΠ̃ ′ ) eq. (60) ===⇒
Tr ( V ′ΣV TpD −2) =ε2λi1(p− σ̃−1(i1) + 1) + ε2λik(p− σ̃−1(ik) + 1) eq. (61)===⇒eq. (62) Tr ( V ′ΣV TpD
−2) =ε2λi1(p− i1 + 1) + ε2λik(p− im + 1), which is eq. (58) as claimed.
Finally, we have to show that the third and the forth terms of the eq. (31) are canceled. First, observe that
Tr ( V ′UNpΠD ( Sp ◦ ( D−1Π′U ′NpΣV D −2 ))) =
Tr ( εDΠ̃′(U ′i1;i1 −U ′ik;ik)UNpΠ ( Sp ◦ ( Π′U ′NpΣV D −2 ))) =
εTr ( Π̃′(Ei1 −Eik)Π ( Sp ◦ ( Π′U ′NpΣV D −2 )) D ) =
ε2 Tr ( Π̃′(Ei1 −Eik)Π ( Sp ◦ ( Π′(λi1Ei1 − λikEik)Π̃ ))) =
ε2 Tr ( (Ei1 −Eik) (( ΠSpΠ̃ ′ ) ◦ ( ΠΠ′(λi1Ei1 − λikEik)Π̃Π̃′ ))) =
ε2 Tr (( ΠSpΠ̃ ′ ) ◦ ((Ei1 −Eik) (λi1Ei1 − λikEik)) ) =
ε2 Tr (( ΠSpΠ̃ ′ ) ◦ (λi1Ei1 + λikEik) ) , and
Tr ( V ′UNpΠD ( Sp ◦ ( D−2V ′ΣUNpΠD −1))) =
Tr ( εDΠ̃′(U ′i1;i1 −U ′ik;ik)UNpΠ ( Sp ◦ ( D−1V ′ΣUNpΠD −1))) = εTr ( Π̃′(Ei1 −Eik)Π ( Sp ◦ ( D−1V ′ΣUNpΠ ))) =
ε2 Tr ( (Ei1 −Eik)Π ( Sp ◦ ( Π̃′(λi1Ei1 − λikEik)Π )) Π̃′ ) =
ε2 Tr ( (Ei1 −Eik) (( ΠSpΠ̃ ′ ) ◦ ( ΠΠ̃′(λi1Ei1 − λikEik)ΠΠ̃′ ))) =
ε2 Tr ( (Ei1 −Eik) (( ΠSpΠ̃ ′ ) ◦ ( Π(i1ik)(λi1Ei1 − λikEik)Π(i1ik) ))) =
ε2 Tr ( (Ei1 −Eik) (( ΠSpΠ̃ ′ ) ◦ ((λi1Eik − λikEi1)) )) =
ε2 Tr (( ΠSpΠ̃ ′ ) ◦ ((Ei1 −Eik)(λi1Eik − λikEi1)) ) =
−ε2 Tr (( ΠSpΠ̃ ′ ) ◦ (λi1Eik + λikEi1) ) =
−ε2 Tr (( ΠSpΠ̃ ′ ) ◦ (λi1Eik + λikEi1) ) .
Now, note that in both cases the matrices that are multiplied elementwise with ΠSpΠ̃′ are diagonal and hence, we only need to look at diagonal elements of ΠSpΠ̃′. Moreover,
ΠSpΠ̃ ′ = Π(i1ik)Π(iki2) · · ·Π(ikik−1)Π̂SpΠ̂′Π(ikik−1) · · ·Π(iki2),
where, i1 · · · ik are fixed points of permutation corresponding to Π̂ so Π̂SpΠ̂′ has the same values at diagonal positions i1 and ik as the original matrix Sp. The only permutation that is only on the left side is Π(i1ik) which exchanges the i1 and ik rows of Sp. Since Sp is such that the elements at each row before the diagonal element are the same and ik > i1, we have the i1 and ik diagonal elements of ΠSpΠ̃′ have the same value. Let that value be denoted as s. Then the sum of the above two equations yields m(λi1 + λik)−m(λi1 + λik) = 0, as claimed.

heading: B DERIVATIVES OF THE LOSS FUNCTION

heading: B.1 FIRST AND SECOND ORDER FRÉCHET DERIVATIVE
text: In order to derive and analyze the critical points of the cost function which is a real-valued function of matrices we use the first and second order Fréchet derivatives as described in chapter 4 of Zeidler (1995). For a function f : Rn×m → R the first order Fréchet derivative at the point A ∈ Rn×m is a linear functional df(A) : Rn×m → R such that
lim V→0 |f(A + V )− f(A)− df(A)V | ‖V ‖F = 0,
where we used the shorthand df(A)V ≡ (df(A))(V ). Similarly, the 2nd derivative is a bilinear functional d2f(A) : Rn×m × Rn×m → R such that
lim V→0 |df(A + V )K − df(A)K − d2f(A)V K| ‖V ‖F = 0,
for all ‖K‖F ≤ 1, where again d2f(A)V K ≡ (d2f(A))(V ,K). The generalized Taylor formula then becomes:
f(A + V ) = f(A) + df(A)V + 1
2 d2f(A)V 2 + o(‖V ‖2),
Moreover, we derive functions ∇f : Rn×m → Rn×m and H(A) : Rn×m → Rn×m such that df(A)V = 〈∇f(A),V 〉F and d2f(A)V 2 = 〈H(A)V ,V 〉F , where again H(A)V ≡ H(A)(V ). Then clearly, A ∈ Rn×mis a critical point of f iff ∇f(A) = 0 and for such As the sign of the bilinear form 〈H(A)V ,V 〉over directions V determines the type of the critical point.
Extending the generalized Taylor theorem of Zeidler (1995), the second order Taylor expansion for the loss L(A,B) is then given by
L(A + V ,B + W )− L(A,B) =dAL(A,B)V + dBL(A,B)W + 1
2 d2AL(A,B)V 2
+dABL(A,B)V W + 1
2 d2BL(A,B)W 2 +RV ,W (A,B),
(63)
where, if ‖V ‖F , ‖W ‖F = O(ε) then ‖R(V ,W )‖ = O(ε3). Clearly, as at critical points where dAL(A,B)V +dBL(A,B)W = 0, as ε→ 0 we haveRV ,W (A,B)→ 0 and the sign of the sum of the second order partial Fréchet derivatives determines the type of the critical point very much similar to second partial derivative test for two variable functions. However, here for local minima we have to show the sign is positive in all directions and for saddle points have to show the sign is positive in some directions and negative at least in on direction. Finally, note that the smoothness of the loss entails that Fréchet derivative and directional derivative (Gateaux) both exist and (foregoing some subtleties in definition) are the same.

heading: B.2 FIRST AND SECOND ORDER DERIVATIVE OF THE LOSS WRT TO B
text: Lemma 5. The first and second (partial Fréchet ) derivative of the loss L(A,B) wrt to B is derived as follows.
dBL(A,B)W= −2 Tr (W ′ (TpA′Σyx − (Sp ◦ (A′A))BΣxx)) (64) = −2〈TpA′Σyx − (Sp ◦ (A′A))BΣxx,W 〉F . (65)
d2B2L(A,B)W 2 = 2〈(Sp ◦ (A′A))WΣxx,W 〉F = 2 Tr (W ′ (Sp ◦ (A′A))WΣxx) . (66)
Proof. Directly compute
L(A,B + W ) = p∑ i=1 ‖Y −AIi;p(B + W )X‖2F
= p∑ i=1 〈Y −AIi;p(B + W )X,Y −AIi;p(B + W )X〉F
= p∑ i=1 〈Y −AIi;pBX,Y −AIi;pBX〉F + p∑ i=1 〈Y −AIi;pBX,−AIi;pWX〉F
+ p∑ i=1 〈−AIi;pWX,Y −AIi;pBX〉F + p∑ i=1 〈−AIi;pWX,−AIi;pWX〉F
= L(A,B)− p∑ i=1 2〈Y −AIi;pBX,AIi;pWX〉+O(‖W ‖2F ) =⇒
L(A,B + W )− L(A,B) = −2 p∑ i=1 〈Y −AIi;pBX,AIi;pWX〉F +O(‖W ‖2F ) W→0 =⇒
dBL(A,B)W = −2 p∑ i=1 Tr(X ′W ′Ii;pA ′(Y −AIi;pBX))
= −2 Tr ( W ′ (( p∑ i=1 Ii;p ) A′Y X ′− ( p∑ i=1 Ii;pA ′AIi;p ) BXX ′ )) = −2 Tr (W ′ (TpA′Y X ′ − (Sp ◦ (A′A))BXX ′)) ,
which can be written as the given form. For the second derivative wrt B we have
dBL(A,B)W = −2〈TpA′Σyx − (Sp ◦ (A′A))BΣxx,W 〉F =⇒ dBL(A,B + W̄ )W = −2〈TpA′Σyx − (Sp ◦ (A′A)) (B + W̄ )Σxx,W 〉F
= −2〈TpA′Σyx − (Sp ◦ (A′A))BΣxx,W 〉F + 2〈(Sp ◦ (A′A)) W̄Σxx,W 〉F =⇒
dBL(A,B + W̄ )W − dBL(A,B)W = 2〈(Sp ◦ (A′A)) W̄Σxx,W 〉F , which by having W̄ → 0 results in the second order partial derivative.

heading: B.3 FIRST AND SECOND ORDER DERIVATIVE OF THE LOSS WRT TO A
text: Lemma 6. The first and second (partial Fréchet ) derivative of the loss L(A,B) wrt to A is derived as follows.
dAL(A,B)V = −2〈ΣyxB′Tp −A (Sp ◦ (BΣxxB′)) ,V 〉F , (67) d2ABL(A,B)V W= −2〈ΣyxW ′Tp −A (Sp ◦ (BΣxxW ′))−A (Sp ◦ (WΣxxB′)) ,V 〉F ,
(68)
d2A2L(A,B)V 2= 2〈V (Sp ◦ (BΣxxB′)) ,V 〉F . (69)
Proof. Directly compute
L(A + V ,B) = p∑ i=1 〈Y − (A + V )Ii;pBX,Y − (A + V )Ii;pBX〉F
= p∑ i=1 〈Y −AIi;pBX,Y −AIi;pBX〉F − p∑ i=1 〈Y −AIi;pBX,V Ii;pBX〉F
+ p∑ i=1 〈−V Ii;pBX,Y −AIi;pBX〉F + p∑ i=1 〈−V Ii;pBX,−V Ii;pBX〉F
= L(A,B)− p∑ i=1 2〈Y −AIi;pBX,V Ii;pBX〉F + p∑ i=1 〈V Ii;pBX,V Ii;pBX〉F
L(A + V ,B)− L(A,B) = − p∑ i=1 2〈Y −AIi;pBX,V Ii;pBX〉F +O(‖V ‖2F ) V→0 =⇒
dAL(A,B)V = − p∑ i=1 2〈Y −AIi;pBX,V Ii;pBX〉F
= −2 Tr(V ′(ΣyxB′ p∑ i=1 Ii;p −A p∑ i=1 Ii;pBΣxxB ′Ii;p)) =⇒
dAL(A,B)V = −2〈ΣyxB′Tp −A (Sp ◦ (BΣxxB′)) ,V 〉F =⇒ dAL(A + V̄ ,B)V = −2〈ΣyxB′Tp − (A + V̄ ) (Sp ◦ (BΣxxB′)) ,V 〉F
dAL(A + V̄ ,B)V − dAL(A,B)V = 2〈V̄ (Sp ◦ (BΣxxB′)) ,V 〉F V̄→0=⇒ d2A2L(A,B)(V , V̄ ) = 2〈V̄ (Sp ◦ (BΣxxB′)) ,V 〉F =⇒
d2A2L(A,B)V 2 = 2〈V (Sp ◦ (BΣxxB′)) ,V 〉F
dAL(A,B + W )V =− 2〈Σyx(B + W )′Tp,V 〉F −2〈−A (Sp ◦ ((B + W )Σxx(B + W )′)) ,V 〉F −2〈ΣyxB′Tp −A (Sp ◦ (BΣxxB′)) ,V 〉F =dAL(A,B)V − 2〈ΣyxW ′Tp,V 〉F −2〈−A (Sp ◦ (BΣxxW ′))−A (Sp ◦ (WΣxxB′)) ,V 〉F +O(‖W ‖2F ) =⇒
dAL(A,B + W )V − dAL(A,B)V = −2〈ΣyxW ′Tp,V 〉F − 2〈−A (Sp ◦ (BΣxxW ′))−A (Sp ◦ (WΣxxB′)) ,V 〉F +O(‖W ‖2F ) W→0 =⇒
d2ABL(A,B)V W = −2〈ΣyxW ′Tp −A (Sp ◦ (BΣxxW ′))−A (Sp ◦ (WΣxxB′)) ,V 〉F .

Review summaries:
This paper proposes a new loss function to compute the exact ordered eigenvectors of a dataset. The loss is motivated from the idea of computing the eigenvectors sequentially. However doing so would be computationally expensive, and the authors show that the loss function they propose (sum of sequential losses) has the same order (constant less than 7) of computational complexity as using the squared loss. A proof of the correctness of the algorithm is given, along with experiments to verify its performance.

This paper proposes a new loss function for performing principal component analysis (PCA) using linear autoencoders (LAEs). With this new loss function, the decoder weights of LAEs can eventually converge to the exact ordered unnormalized eigenvectors of the sample covariance matrix. The main contribution is to add the identifiability of principal components in PCA using LAEs and. Two empirical experiments were done to show the effectiveness of proposed loss function on one synthetic dataset and the MNIST dataset.