Tensor decompositions [36] provide powerful tools to represent and compress higher-order arrays. In this work, we focus on the symmetric rank-r Tucker decomposition. Applied to a d-th order symmetric tensor of size n → . . . → n, it expresses the input using a symmetric core tensor C C C ∈ R r→•••→r (r n) and a matrix Q ∈ R n→r with orthonormal columns. See Figure 1 for an illustration. Symmetric tensors abound in statistics as a means of representing higher-order moments and cumulants of a random vector. They are commonly used to analyze dependences in multivariate distributions [12,30,14], learn latent variable models via the method-of-moments [4,24], solve important signal processing problems [9,10], etc.
In most cases, the exact moment of the ground-truth distribution is not available. The idea of moment estimation is to use the d-th order sample moment tensor, denoted M M M d , instead. This is formed from observations {x 1 , . . . , x p } ⊆ R n of a random variable x as
There is a need for economical representations of M M M d , which would otherwise take O(n d ) space to store. Using the symmetric Tucker decomposition, one approximates M M M d solely by determining the rank-r factor matrix Q, which takes just rn memory space. (The core C C C is a byproduct once Q is identified; see Remark 1.) However, decomposing the sample moment still has upfront costs of O(pn d ) and O(n d ) to initially construct and store M M M d . For large data volumes p and high dimensionalities n, we turn to implicit (or tensor-free) and streaming strategies to escape such prohibitive overhead. The symmetric Tucker decomposition of moment tensors has other important motivations, besides compression. One major application is to extract the key features Q ∈ R n→r from moment tensors. Without the need for the core tensor C C C ∈ S d (R r ), this usage further reduces the output's storage. In some scientific domains, data has high volatility, e.g., in combustion simulations [1], outlier detection [25,55] and asset allocation [49,8,6]. Higher-order moments and cumulants (d ≥ 3) are sensitive to heavy tails and anomalies. Thus, finding the key subspace Q of variability in higher-order statistics is important for non-Gaussian (also known as non-normal) data analysis. This is akin to the role of principal component analysis (PCA) and the covariance matrix for Gaussian data.
In other settings, data observations are assumed to lie in a low-dimensional subspace and follow the low-rank factor model [61]. For example, this holds in the Fama-French model [22] of finance and also in psychometrics [48]. The resulting moment tensors naturally admit a Tucker format. Thus the decomposition enforces the low Tucker rank structure and improves moment estimation accuracy.
Due to the NP-hard nature of tensor decompositions [29], an approximate symmetric Tucker decomposition is computed by higher-order eigenvector decomposition (HOEVD) [15]. It is the eigendecomposition of the flattened tensor. It can be viewed as the symmetric version of higher-order singular value decomposition (HOSVD) [14] and efficient variants [64,43,62,45,11,2]. Despite their popularity, HOEVD o"ers no guarantee of returning the exact critical points for the problem.
We instead turn to iterative search methods to achieve the best possible solution, using HOEVD as an initialization. Newton [21,33], quasi-Newton [57] and trust-region [32] methods have been considered for general Tucker decomposition. These second-order optimizations rely on Riemannian Hessian updates on a the manifold of orthonormal matrices. For symmetric tensors, higher-order orthogonal iterations (HOOI) [16], Jacobi rotations [31,41] and power method [54] were proposed. Regalia briefly mentions a shifted variant of the power method in [54], which is akin to a gradient descent method, although strong convergence analysis is omitted in the work.
Our work is primarily motivated by moment tensor decompositions, which are fundamental objects in statistics and other disciplines. Neither the existing HOEVD nor iterative algorithms for Symmetric Tucker take advantage of the sample moment tensor structure (1). We show that streaming and implicit techniques can be applied to address the computational burden in the data-rich and high-dimensional regime. For covariance matrices which are centered second-order moments, this is relevant to the streaming PCA scheme, also known as Oja's method [51], and its adaptations [46,67,28].
Recently, a similar strategy was applied to optimization with moment tensor CP decomposition for learning Gaussian mixture models, see [59,52]. Building on techniques from [59], the present work also passes over only a subset of samples at a time and drives the optimization without forming gigantic moment tensors. Despite using similar approaches for increasing efficiency, we emphasize that CP and Tucker decompositions are useful for their own independent purposes. The motivations we list for feature extraction and data factor analysis specifically apply to Tucker decompositions. We refer readers to Section 4.5 for a concrete comparison.
In this paper, our major contribution is developing scalable algorithms for symmetric Tucker decomposition on moment tensors.
To start with, we propose the projected gradient descent (PGD) framework to the general symmetric Tucker tensor decomposition problem, see Algorithm 1. The approach is led by a first-order gradient update governed by the Tucker approximation cost. It is free of non-trivial second-order updates on manifolds.
Driven by their importance in data science, our attention then focuses on decomposing sample moment tensors (1). By exploiting the special input structure, we develop scalable PGD (SPGD) in Algorithm 2.
There are two main highlights. First, we adopt the implicit technique and operate on the input sample vectors directly. This much cheaper in memory and computation than constructing the sample moment tensor. Proposed in [59,52], the implicit strategy greatly reduces the e"ort for gradient computation. Second, SPGD is built on a streaming model. By passing over small batches of data observations, both storage and computational complexity are further reduced. Similar ideas have been used in CP moment tensor decomposition [59] for separate purposes.
We also design a scalable implementation of HOEVD (SHOEVD) for sample moment tensors (1), see Algorithm 3. This method is akin to the streaming PCA idea. We use the SHOEVD approximation (Algorithm 3) as the initialization to SPGD (Algorithm 2). Moreover, SHOEVD may be of independent interest due to its broad usages, e.g. [13,3]. The savings from the streaming and implicit strategies allow SHOEVD and SPGD to handle large-scale computations. The gains are summarized in Table 1. Here p is the total sample number and b is the batch size at each iteration. Usually, the user sets b p. These efficiency advantages are corroborated by extensive numerical experiments.
We demonstrate the key motivations of sample moment decompositions, including feature extraction and moment tensor low-rank estimation. Our proposed SPGD is respectively applied to real datasets in anomaly detection and portfolio allocation.
In terms of theoretical guarantee, the basic PGD (Algorithm 1) convergence result is provided. Employing the theories from manifold optimization and dynamical systems, we prove that for any symmetric tensor, PGD converges monotonically to a first-order critical point given arbitrary initialization. Furthermore, it converges to a second-order critical point almost surely. Scalars, vectors, matrices and tensors are represented as lowercase letters x, lowercase bold letters x, uppercase bold letters X, and calligraphic bold letters X X X, respectively. The (i 1 , . . . , i d )-th element of a d-way tensor is denoted by x i1,...,i d . The j-th column of a matrix X is x j . The spectral norm of a matrix is • 2 . The norm of a tensor is • , see Definition 5. The identity matrix is I n ∈ R n→n . The Moore-Penrose pseudoinverse of a matrix is denoted using †. The matrix symmetrization operator is written sym : R n→n ⇐ S 2 (R n ) and given by sym(X) =foot_0 2 (X + X ) for X ∈ R n→n . We write orth for the function that selects the orthonormal factor Q ∈ R n→r from the QR decomposition of a full-rank n → r matrix. 1 The notation X [d] is the elementwise exponentiation, where each coordinate of X is raised to the d-th power. The set {1, . . . , n} is denoted by [n] when n is a positive integer.
Definition 1 (symmetric tensors). A tensor X X X ∈ R n→•••→n is symmetric if its entries are invariant under any permutation of indices, i.e.
where Π(d) is the permutation group on [d]. Here, d and n are called the order and dimension of X X X. The space of all real-valued symmetric tensors of order d and dimension n is denoted by S d (R n ).
Definition 2 (symmetric outer product). The outer product of a vector x ∈ R n with itself d times is denoted as x ⊗d ∈ S d (R n ). The (i 1 , . . . , i d )-th entry of this product is
Definition 3 (symmetric Tucker product). For a symmetric tensor X X X ∈ S d (R n ) and a matrix Y ∈ R n→m , their symmetric Tucker product is defined as X X X with Y multiplied along each mode, written X X X • (Y, . . . , Y) ∈ S d (R m ). The (j 1 , . . . , j d )-th entry of the product is
, their inner product in all modes except the first one [21,17] is defined as:
where z j1,j 1 =
Let X X X ∈ S d (R n ) be a symmetric tensor. Fix a user-specified rank r (typically r n). Our goal is to find the best approximation of X X X represented by the symmetric Tucker product of a core tensor C C C ∈ S d (R r ) and a rank-r orthonormal matrix Q ∈ R n→r . The problem is formulated as follows:
In fact, the minimization (2) is equivalent to solving for a rank-r basis Q that maximizes the following cost function: max
Remark 1. The problem formulation transforms from (2) to ( * ) by solving the unconstrained linear least squares in (2) for the core, i.e.,
As X X X 2 is a constant in Q, minimizing the above quantity is inded the same as maximizing X X X • (Q, . . . , Q) 2 .
Remark 2. The full Tucker decomposition (2) takes O(nr + r d ) in memory to store the basis Q ∈ R n→r and the core C C C ∈ S d (R r ). Though still exponential in d, this cost savings can be substantial compared to the original X X X ∈ S d (R n ) requiring O(n d ) space if the rank r is much less than the ambient dimension n. Further, in applications such as tensor PCA [47] and anomaly detection Section 5.1, one only desires the feature subspace Q, and the core C C C need note be computed or stored. This reduces the cost to O(nr).
The feasible domain in ( * ) for Q is the Stiefel manifold:
which is a compact Riemannian submanifold of R n→r with dimension r(2n-r-1)
. Further, the cost function F is rotationally invariant. In fact, it only depends on the orthogonal projector P = QQ ∈ S 2 (R n ), since X X X • (Q, . . . , Q) = X X X • (P, . . . , P) . So, the optimization ( * ) naturally runs on the Grassmannian manifold, which consists of all orthogonal projectors P:
Here, Gr(n, r) is a compact Riemannian manifold of S 2 (R n ) with dimension r(nr).
Remark 3. No simple characterization of the critical points of ( * ) is known when d ≥ 3. This is in contrast to case of d = 2, where the classical Eckart-Young theorem [20] holds. Various possible extensions of Eckart-Young to tensors have been shown to fail, e.g., see [35].
This paper develops a projected gradient descent (PGD) method to solve the constrained optimization problem ( * ). The PGD framework works as follows: given an initial point, we move the current iterate Q t-1 along the negative gradient direction by a step size of ε t > 0 where the gradient is:
Then we project the update back to the Stiefel manifold (4) by selecting the Q factor from the QR decomposition. The procedure is described in Algorithm 1. See Section 6 for a convergence analysis of PGD.
An approximate solution to the best symmetric Tucker decomposition ( * ) is the HOEVD [15] result:
given by the leading r eigenvectors of mat(X X X) mat(X X X) .
HOEVD is widely used for Tucker approximation, beause it can be computed via a direct matrix eigendecomposition.
Generally speaking, the HOEVD solution is not a critical point of ( * ) [15]. Its approximation error to a d-th order tensor X X X is bounded [26, Theorems 6.9-10] as
Here, X X X hoevd = X X X • (Q hoevd Q hoevd , . . . , Q hoevd Q hoevd ) as in Remark 1 and X X X best is the optimum of (2). We observe the suboptimality of the HOEVD solution in Section 4.4. However according to [14, section 3.4], the HOEVD solution usually belongs to an "attractive" region around the desired local optima in numerical experiments. So it is a good candidate for initialization.
The main interest in this work is the decomposition of the sample moment tensor M M M d ∈ S d (R n ) (1). We design scalable algorithms for the symmetric Tucker decomposition ( * ), and the HOEVD approximation (7), of M M M d . An average of symmetric outer products, the sample moment structure (1) is the key to the scalability of the decomposition computation. Tensor operations with outer products can be reduced to much cheaper computations involving only their factor vectors. Specifically we use the following identities: for vectors x, y ∈ R n , z ∈ R m and a matrix Y ∈ R n→m , it holds [26] x
x ⊗d , y ⊗d ⇔ = x, y⇔ d ∈ R, and
In Algorithm 2, we introduce the scalable PGD (SPGD) when the input tensor is a sample moment (1). The symmetric Tucker cost ( * ) now depends solely on the data observations X = [x 1 , . . . , x p ] ∈ R n→p . Indeed due to (8),
We actually do not need to create or store the sample moment M M M d . Such an idea was called the implicit technique in [59,52].
Commonly, there is a huge amount of data in exceedingly high dimensions. We hence further turn to the streaming strategy to speed up the computations. Here data is revealed as a sequence of sample batches {X 1 , . . . , X t , . . . } ⊆ R n→b with batch size b (b p). We read each batch as it arrives but discard it afterwards. This streaming model is called the turnstile model [50,42,63].
Therefore, as the t-th (t ≥ 1) batch X t = [x t,1 , . . . , x t,b ] ∈ R n→b arrives, the cost function (11) is recast stochastically:
In light of (10), we can compute the stochastic gradient in a tensor-free manner:
Here, ( • ) [d-1] means raising each element of the input matrix to the (d -1)-th power. Computing the gradient implicitly via online samples reduces the cost to O(rnb + rb 2 ) flops per iteration. This improves upon O(rnp + rp 2 ) using the full amount of samples, and a cost that is exponential in d if we operate on the explicit tensor in S d (R n ). Moreover, a scalable HOEVD solver, shoevd in line 1 of Algorithm 2, is utilized for good initialization; see Section 3.2 for details. We also use AdaGrad [19] to automatically tune step sizes. We implement the column-wise extension of AdaGrad proposed in [28], since it is natural to independently update the orthonormal basis vectors in Q. The step size (ε ε ε) t becomes a row vector in R 1→r , rather than a scalar. See lines 4 and 5, where ./ denotes columnwise division.
SPGD is a two-phased iterative algorithm. It consists of
• Phase I -SHOEVD approximation (line 1),
• Phase II -symmetric Tucker decomposition (lines 2-8).
Algorithm 2 SPGD for sample moment Tucker decomposition Input:
step size 5:
We present the scalable HOEVD (SHOEVD) method in Algorithm 3. This serves as the Phase I initialization in SPGD (Algorithm 2), but it may be of independent interest. As HOEVD is essentially a matrix eigendecomposition, we follow the streaming PCA idea. The solver uses the PGD framework, but is applied to a di"erent cost function than (11). Specifically, we solve HOEVD iteratively finding Q ∈ St(n, r) that maximizes the following stochastic cost formed by a data batch
Using (10), the gradient can be calculated implicitly:
It costs O(rnb 1 + nb 2 1 + rb 2 1 ) flops to compute the gradient. Remark 4. SPGD and SHOEVD (Algorithms 2 and 3) obviate the upfront cost to form the tensor, which would require O(pn d ) flops. Each iteration temporarily stores an updated matrix Q t ∈ R n→r and an online sample batch X t ∈ R n→b , taking only nr + nb memory. This explains the complexity counts shown in Table 1.
for t = 1, . . . , T 1 do 5:
:
t + colnorm(G t ) [2] [1/2] 7:
end for 3) can be evaluated as follows. After obtaining Q ∈ St(n, r) by SPGD or SHOEVD (Algorithms 2 and 3), we draw fresh data X c = [(x c ) 1 , . . . , (x c ) pc ] ∈ R n→pc and then set:
The core computation costs O(r d ) in storage and O(p c rn + p c r d ) in flops.
In this section, we utilize synthetic tests to study the proposed algorithms for sample moment tensor decomposition. And real-world datasets are considered in Section 5. In this section, we show:
• the scalability of implicit compared to explicit methods in Figure 2,
• the speedup of streaming compared to full data usage in Figure 3,
• the robustness of the HOEVD approximation and the symmetric Tucker decomposition to noise and target rank choice in Figures 4 and 5, and
• a comparison between symmetric Tucker and CP in Table 2.
In these tests, we generate moment tensors with approximate low-rank symmetric Tucker structure (3). One way is to assume data instances x ∈ R n follow a linear factor model:
Here we fix the factor loading matrix B ∈ R n→rtrue . The latent factor f ∈ R rtrue is drawn from a coordinateindependent, zero-mean and unit variance normal-inverse Gaussian (NIG) distribution [5,6]. This distribution is skewed and heavy-tailed, aligning with some motivations discussed in Section 1. The noise
Note that in (15), the d-th moment of x roughly admits a low-rank Tucker format:
We define the inverse signal-to-noise ratio
The above quantity (rather than the deviation σ) is to measure the model's noise level.
Remark 6. Our goal is to simply perform the symmetric Tucker decomposition of the sample moment tensor M M M d ∈ S d (R n ). This is different from inference procedures on the linear factor model, e.g., independent component analysis (ICA) [14].
Remark 7. Besides the model (15) or other heavy-tailed datasets, we emphasize that our proposed algorithms are able to efficiently decompose M M M d when the sample vectors come from any probabilistic model.
All numerical tests are implemented using MATLAB R2022a on a MacBook Pro with 16GB RAM memory. Each reported plot is averaged over 5 independent simulations with independent, random initializations. The matrix B in ( 15) is randomly generated with i.i.d. standard normal entries and then fixed for all runs. We generate new f and ε ε ε ε ε ε ε ε ε to form the data samples in each stochastic simulation. The parameters of each experiment, namely n, r true , r, p, d, SNR -1 , are specified below.
We illustrate the scalability of implicit and explicit implementations in Figure 2. The experiments are nonstreaming. The implicit methods run on a fixed set of observations X ∈ R n→p . The explicit approaches require constructing the sample moment tensor M M M d ∈ S d (R n ) as overhead, and then operating on it. Due to resource limitations, the maximal tensor dimension n for explicit methods is smaller than maximal dimension for the implicit methods. We set the sample size to be p = 5n and choose the rank as r = r true = 5. We consider non-streaming methods, so b 1 = b 2 = p. The sample moment formation uses the symktensor package in the Tensor Toolbox for MATLAB [37]. The noise level is SNR -1 = 0.5, see (16). We conduct the experiments in two phases for order-3 and 4 sample moments M M M d . For Phase I (initialization), we compare the efficiency of the implicit HOEVD result,
and its explicit counterpart (7). The equivalence of the two methods is by (10). Then in Phase II (iterations), we run the SPGD Algorithm 2 and various explicit algorithms, namely the explicit PGD Algorithm 1, HOOI [15], quasi-Newton [57] and trust-region [32] methods. 2For the last three algorithms, we use the code tucker als.m in Tensor Toolbox [37], symalgQNlc.m in package [56] and lmlra3 rtr.m in TensorLab [65], respectively. HOOI and trust-region apply to general asymmetric tensors. We omit the trust-region method in d = 4 case as its code does not apply to the fourth order. The step size for Algorithm 2 is c 2 = 1000. Note that the implicit Algorithm 2 matches the explicit Algorithm 1 iteration-by-iteration, because there is no streaming. However, they calculate the updates di"erently because Algorithm 2 is implicit.
Figure 2 shows the running time of the two phases as the dimension n increases. The time on the vertical axes shows when the algorithm reaches a relative gradient to be less than 10 -12 . This measure of convergence precision is used in [21,57,32] and defined as:
Here grad F (Q) = (I n -QQ )∇F (Q) ∈ R n→r denotes the Riemannian gradient on the Stiefel manifold (4), where ∇F (Q) ∈ R n→r is the Euclidean gradient (6). We use the Riemannian gradient grad F (Q) because the possible update direction on the Stiefel manifold is constrained to the tangent space at Q. The time for checking the relative gradient is included in the runtime. It is evident from Figure 2 that the implicit HOEVD (17) and SPGD (Algorithm 2) are much faster than the explicit algorithms. The advantage of being tensor-free allows the implicit implementation to handle computation in higher dimensions. Meanwhile, the prohibitive upfront and storage costs in the explicit methods restrict their scalability.
We stress the gains provided by the streaming technique. This is demonstrated in stochastic experiments for (a) SHOEVD (Algorithm 3) and (b) SPGD (Phase II of Algorithm 2) in Figure 3. The vertical axes represent the approximation accuracy evaluated by the metric:
The measure is rotationally invariant. We plot the subspace error against runtime.
To compare streaming and full data methods, we first generate a pool of observations X ∈ R n→p . Each streaming method selects a subset of samples X(:, b(t-1)+1 : bt) from the pool at step t. We ensure the total sample number p is just enough for the streaming case with the largest batch size b max . The non-streaming counterparts use all the samples X at once. We set the number of iterations T = 250, and the total sample size p = T b max = 25, 000.
In Phase I (a), the full HOEVD refers to the implicit way (17), which is already more efficient than explicit HOEVD (7). The step size for all streaming cases of Algorithm 3 is c 1 = 0.05. In Phase II (b), all plots start at the same initialization by Algorithm 3 with T 1 = 20, b 1 = 50 and c 1 = 1. These parameters for Phase I are di"erent from the (a) setting. They do not a"ect the efficiency of (b). The streaming and full cases for the actual Phase II use step size c 2 = 1 and T 2 = 10 respectively.
It is seen in Figure 3 that the streaming Algorithms 2 and 3 have limited loss of accuracy compared to the methods that use the full dataset. On the other hand, the runtimes for all streaming cases are nearly negligible compared to timings for the full methods.
We test how the noise level in the observations and the choice of target rank a"ect HOEVD and symmetric Tucker decomposition in Figures 4 and 5. In the two experiments below, we plot the evolution of both phases of Algorithm 2, see Figures 4 and 5. As the HOEVD phase (in dashed lines) plateaus, its last iterate gives the initialization for Phase II (in solid lines). The number of iterations are shown on the horizontal axes. The parameters in Algorithm 2 are set to be
We first consider the algorithms' sensitivity to noise in the data samples. To test this, we generate a sequence of clean samples and then plant Gaussian noise with di"erent levels of SNR -1 . The algorithms' results are shown in Figure 4. In particular, note the improvement in Phase II over the HOEVD phase. The discrepancy is enlarged in the high noise regime (SNR -1 ≥ 1), when the HOEVD solver provides a poor solution. In contrast, relatively good approximation is obtained using the symmetric Tucker decomposition (Phase II) when SNR -1 is no bigger than 1.
Next, we test various rank choices in Figure 5. All cases pass over the same set of data observations. The algorithms' recovery qualities are reported in terms of the objective value F ( * ), and compared to the "ground-truth" value F (Q * ). As there is no deterministic moment M M M d in the streaming case, we estimate F using an additional set of testing samples X test ∈ R n→ptest drawn from (15). We calculate F implicitly using (9):
When r ≥ r true , the solutions obtained by Algorithm 2 fully approximate the moment tensor in the sense that their objective values match the ground truth. When r < r true , we can only partially recover the tensor. The less the rank is, the lower the objective value. Also note that as r decreases, HOEVD's initialization moves further away from the ground truth and the gap between the phases is more pronounced.
Figures 4 and 5 show that the robustness of the symmetric Tucker decomposition is stronger than that of the HOEVD approximation. For "hard" problem instances where noise is high or the rank is under-specified, the HOEVD phase (in dashed lines) is clearly improved upon during Phase II (in solid lines).
This subsection shows the situations where symmetric Tucker decomposition is more appropriate to use than the symmetric CP counterpart. We provide numerical evidence in the moment tensor approximation quality for low-rank factor model, while maintaining a low storage budget.
We employ the relative error as the accuracy metric:
The testing sample moment M M M d test ∈ S d (R n ) serves as the reference solution and is built from the dataset X test ∈ R n→ptest via (1). The approximation tensor M M M d ∈ S d (R n ) is obtained via either the symmetric Tucker or the symmetric CP [59] decompositions. These solutions admit the formats:
In (18), we denote the Tucker and CP rank by r and r respectively. The Tucker approximation is determined by the basis Q ∈ R n→r and its core C C C ∈ S d (R r ) is via (3). The CP format consists of the coefficients {λ i } r i=1 ⊂ R and the vector components {v i } r i=1 ⊂ R n . For the Tucker and CP solutions, we respectively apply implicit stochastic methods: our Algorithm 2 and Algorithm cp isym.m [59] with Adam optimization from Tensor Toolbox [37].
The data observations are drawn from the linear factor model (15). Here, the latent factor f ∈ R rtrue follows the standard Gaussian distribution with zero mean and identity covariance. The noise level is SNR -1 = 0.05. In one simulation, we first draw the full data set X test ∈ R n→ptest with dimension n = 500 and sample size p test = 10000. We form order d = 4 moment tensors. At each iteration of both algorithms, we compute the gradient based on a sample batch X t ∈ R n→b with batch size b = 50. Their objectives are evaluated based on the full testing tensor M M M d test . The algorithms take a complete pass over X test in one simulation (for Tucker) and one epoch (for CP). All epochs of one CP simulation operate on the same data X test .
In Table 2, we compare the accuracies of the symmetric Tucker versus CP decompositions, in reference to M M M d test . Regarding the target rank, the Tucker method sets it equal to the true rank: r true . There are two choices for the CP method: r true and r d-1 true . For each table entry, we report the lowest error from 5 independent simulations. We see that the Tucker estimation achieves the greatest accuracy on a low rank choice. On the other hand, CP decomposition su"ers substantially larger error (10000 times worse than Tucker) when its target rank is small. CP manages to fit moment tensors with high precision similar to Tucker's level for r = r d-1 true . However, comparing the storage costs, Tucker is much more economical in representing the moment tensor. From (18), Tucker costs rn + r+d-1 d ≈ r true n to store Q and C C C, whereas the CP decomposition costs r (n + 1) ≈ r d-1
true n in storage to achieve a similar approximation. We conclude that there are important situations where the Tucker structure provides a more natural and more accurate format for moment tensors than CP does. This holds especially when the data samples lie close to a low-rank subspace.
This section illustrates the applicability of moment Tucker decompositions to real datasets from anomaly detection and portfolio allocation.
The computations are not wholly streaming, as we need to store and preprocess the data first. However, only data batches are used while running Algorithms 2 and 3. The task results are based on one simulation.
Outliers can be discovered using higher-order statistics, when the mean and covariance information are not sufficient [18,1]. In hyperspectral imaging (HSI) [53], the detection goal is to find a low-rank feature subspace, denoted by Q ∈ St(n, r). After transforming HSI pixels by Q, unusual targets in the image are revealed [25]. Intuitively if the anomalies are captured by colspan(Q), they should stand out as bright spots in the transformed image.
This section tests three methods to extract the feature subspace from HSI data: Tucker factorization of the skewness (order-3 standardized moment) [44], HOEVD applied to the skewness, and eigendecompostion of the covariance. Here, skewness and covariance are formed from HSI pixels. The results are qualitatively and quantitatively evaluated according to their performance in capturing anomalies, as in [25].
Figure 6 shows the data structure of a typical HSI image. It is a 3-way hypercube X X X ∈ R n→p1→p2 . The data consists of n band images in R p1→p2 and p = p 1 p 2 HSI pixels in R n . We use a real hyperspectral image from the Airport-Beach-Urban dataset. 3 The scene was taken at Gulfport, Florida by the Airborne band
Visible/Infrared Imaging Spectrometer (AVIRIS) sensor. The top two images of Figure 7 are the false-color composite image and the detection reference image (with white spots depicting the anomalous aircraft). The number of bands is n = 191 and the size of each band is p 1 → p 2 = 100 → 100. We use scalable algorithms for the moment decomposition task.
First, the HSI cube X X X ∈ R n→p1→p2 is reshaped into an n → p matrix. This is whitened [23] to give X ∈ R n→p . We choose rank r = 4, and seek the leading rank-r subspaces Q ∈ R n→r from the skewness tensor M M M 3 ∈ S 3 (R n ) and the covariance matrix M 2 ∈ S 2 (R n ), both formed using the columns in X. The subspaces are found using SPGD (Algorithm 2), SHOEVD (Algorithm 3), and PCA. The parameters for Algorithms 2 and 3 are b 1 = b 2 = 100, c 1 = 0.35, c 2 = 0.5 and T 1 = 500, T 2 = 1500. For each subspace Q obtained, X is transformed into Q X ∈ R r→p . Each row of Q X is reshaped to a p 1 → p 2 matrix. The resulting bands 1 -4 are shown in Figure 7.
Since the reference map in Figure 7 is sparse, we further pursue a sparsest detection image by optimizing over weights on the bands. This is achieved by a constrained ℓ 1 -minimization:
This is solved by the MATLAB Optimization Toolbox function fmincon. Then we reshape X Qw * into a band image. The resulting sparse images are in the first column of Figure 7. Inspecting Figure 7, we see that SPGD and SHOEVD based on skewness distinguish between abnormal and normal pixels well. The subspace provided by SPGD enables especially clear detection, since the targets are the most pronounced. It outlines both aircraft individually in the sparse image. The covariance information, however, fails to identify anomalies.
false-color image reference map sparsest SPGD (Alg. 2) (d=3) band 1 band 2 band 3 band 4 SHOEVD (Alg. 3) (d=3) PCA (d=2) In Figure 8, we show a quantitative evaluation in anomaly detection, via the receiver operating characteristic-area under curve (ROC-AUC). This quantifies how accurate a method is in classifying normal and abnormal objects, where a higher AUC value indicates better detection. To draw the ROC curve, we whiten the matrix X Q and assign the anomaly score for each pixel as its ℓ 2 -norm. The dashed curve indicates a random classifier. Note that the ROC of skewness-based detection schemes lie strictly above that of PCA. The AUC values corresponding to SPGD, SHOEVD and PCA are 0.9947, 0.9900 and 0.8893.
10 -4 10 -3 10 -2 10 -1 10 0 0 0.2 0.4 0.6 0.8 1 false alarm rate detection probability SPGD (Alg. 2) SHOEVD (Alg. 3) PCA random
Moment estimation is sometimes improved by enforcing a low-rank structure [6]. In portfolio allocation, we seek a distribution over investments to maximize future returns. Judicious predictions can be obtained by including higher-order moment estimations in allocation optimization, as financial data is heavy-tailed [34].
Higher moments of financial data should roughly admit a low-rank symmetric Tucker format, because of the Fama-French model [22]. In this section, we test two ways to calculate moments: the sample moment M M M d (1) and the low-rank symmetric Tucker approximation to M M M d . We evaluate the moments by deriving portfolios based on them and comparing the portfolio performances, as in [6].
We use the real HFR hedge fund dataset. 4 It consists of daily return rates (in percentage) of n = 30 HFRX indices from 01/02/2009 to 12/31/2019. The time series is split into in-sample data (2009 -2016) X in ∈ R n→pin (p in = 2016) and out-of-sample data (2017 -2019) X out ∈ R n→pout (p out = 754). The insample data is centered. We calculate the averaged sample moments (1):
Then their low-rank Tucker approximations are computed by SPGD (Algorithm 2). The rank is set to r = 15, and (3) is used to solve for the cores once the subspaces Q are determined. We run Phase II of Algorithm 2 twice. Starting from a random initialization, we obtain Q for M 2 . This is used an initialization for decomposing M M M 3 . Finally we tackle M M M 4 , starting from the Q for M M M 3 . The parameters are set to c 2 = 1, 10 -2 , 10 -5 ; T 2 = 1000, 1000, 200; and b 2 = 500, 500, 500 for the three runs respectively. This produces
Given moment estimates, we select the allocation weights that maximize the expected utility (EU) [6, Supplementary Material, Eqn. ( 17)]:
where the risk aversion parameter µ is set to be 1. We estimate E (xx) ⊗d ≈ M M M d in (19), and then apply the MATLAB Optimization Toolbox function fmincon to output the optimal weights. Meanwhile, estimating E (xx) ⊗d ≈ M M M d Tucker results in another weight vector.
The two portfolios are evaluated on the out-sample data. The acquired daily returns {y t } ⊆ R are collected in a vector y = X out w * ∈ R pout . In Figure 9, we plot the cumulative daily return rates ( t t =1 (1 + y t 100 ) -1) → 100 at days t = 1, . . . , p out . We observe that the low-rank moment estimates produce a better return rate than the plain sample moments do. Despite mild fluctuations, the two strategies always produce positive return rates. A naive strategy, the equally-weighted portfolio, is the most unstable. It sometimes yields negative returns.
01/03/17 01/02/18 01/02/19 01/02/20 -4 -2 0 2 4 6 8 time return rate(%) cumulative daily return (2017-2019) low-rank moment sample moment equally weighted Other statistical measures are shown in Table 3. The annual geometric mean of returns over the years 2017 -2019 is
(yt-y) d pout . Large geometric mean and low absolute values of moments are desirable allocation outcomes. In each row, we highlight the preferred score in bold. They are all achieved by the portfolio based on low-rank Tucker approximations to the moments.
This section presents a convergence theory for the basic PGD method (Algorithm 1) applied to symmetric Tucker decomposition ( * ). We note several points in the analysis: 1) The tensor X X X needs not arise as a moment tensor; even when it does, we do not analyze the e"ects of streaming. 2) Proving convergence of the non-streaming algorithm to critical points is needed first, because ( * ) is non-convex, with a complex landscape when d ≥ 3 (recall Remark 3). 3) Our results do not immediately follow from known general properties of plain GD or projected GD. We use the specific properties of the Tucker objective function and QR retraction.
Recall the set-up in Section 2.2. Because the Tucker factorization cost ( * ) is non-unique as
the convergence of {Q t } is not relevant. Instead we focus on the convergence of the projectors, P t := Q t Q t ∈ Gr(n, r). Thus, while Algorithm 1 runs on the Stiefel manifold and generates the sequence of bases {Q t } ⊆ St(n, r), we analyze the sequence of subspaces {P t } ⊆ Gr(n, r) on the Grassmannian manifold. Given scalar step sizes {ε t } ⊆ R ≥0 , the update scheme starts from the PGD iteration on Q:
The corresponding update on P is defined as
St(n, r)
Gr(n, r)
The sequence {P t } is assessed in terms of the cost function
As noted in Section 2.2, it is equivalent to the cost in ( * ), as f (P) = F (Q) when P = QQ . We prove that the iterates {P t } converge to a critical point of f on the Grassmannian. This is stronger than just showing convergence of the objective values {f (P t )} in R.
The first and second-order criticality conditions are stated as follows in the manifold setting.
Here, grad f (P) ∈ T P Gr(n, r) ⊆ S 2 (R n ) denotes the Riemannian gradient, which is defined as the Euclidean gradient ∇f (P) ∈ S 2 (R n ) projected to T P Gr(n, r). The tangent space at P is given by
Here, Hess f (P) :
The first theorem is a uniform guarantee: Algorithm 1 always converges to a first-order critical point regardless where the iteration begins. Theorem 6.1. Let X X X ∈ S d (R n ) be any symmetric tensor, and consider the objective function f defined by (22). There exists a constant Γ * (X X X, r) > 0 such that if the step sizes {ε t } ⊆ R ≥0 satisfy sup t ε t ∞ Γ * and inf t ε t > 0, the following holds. For all initializations Q 0 ∈ St(n, r), the sequence {P t } ⊆ Gr(n, r) generated via (20)-( 21) by Algorithm 1 converges monotonically to a first-order critical point (23) of f . The convergence is at no less than an algebraic rate.
Next, in fact, Algorithm 1 to a second-order critical point, for almost all initializations. Theorem 6.2. Let X X X ∈ S d (R n ) be any symmetric tensor, and consider the cost f defined by (22). Use a constant step size ε for all iterates. There exist a constant Γ * * (X X X, r) > 0 and a measure-zero subset I ⊆ St(n, r) such that if 0 < ε ∞ Γ * * , then for all initializations Q 0 ∈ St(n, r) \ I, the sequence {P t } ⊆ Gr(n, r) generated via (20)-( 21) converges to a second-order critical point (25) of f . We build the proofs of Theorems 6.1 and 6.2 respectively in Sections 6.2 and 6.3. The main technical tools are a convergence theorem for real-analytic cost functions (based on the ' Lojasiewicz inequality) and the center-stable manifold theorem. These are stated in Appendix A. Detailed proofs for supporting propositions and lemmas are in the Supplementary Materials. They involve non-trivial calculations for the QR retraction and the symmetric Tucker cost function (22).
Remark 8. The convergence of PGD in Gr(n, r) is not for free. Previously, Regalia found example tensors X X X where our Algorithm 2 with a step size approaching infinity shows an oscillatory cycle, rather than convergence. See [54, Figures 1 and 4].
To establish Theorem 6.1, we apply the guarantee stated in Theorem A.1, originally from [58,Theorem 2.3]. It requires us to prove {P t } meets C1 -C4 conditions in Theorem A.1.
We first compute the Riemannian gradient of f (22) on Gr(n, r). Although Theorem 6.1 is in terms of P ∈ Gr(n, r), our proof uses Q ∈ St(n, r) due to the explicit update rule in (20). Thus, we state the relation between grad f (P) and Q, ∇F (Q). As short-hand notation, define w : Gr(n, r) ⇐ S 2 (R n ) by
Note that w(P) is a positive semidefinite (PSD) matrix.
Proposition 6.3. For Q ∈ St(n, r) and P = QQ ∈ Gr(n, r), the following holds grad f (P) = sym 2 (I n -P) w(P)P = sym
Proposition 6.3 is proven in Appendix B. Next, we analyze the update rule on P ∈ Gr(n, r). For t ≥ 0 and step size ε t > 0, define α t : Gr(n, r) ⇐ S 2 (R n ) by α t (P) = P + 2ε t Pw(P)P ∈ S 2 (R n ).
The expression α t (P) and the gradient grad f (P) help to represent the update on P. Lemma 6.4. Let a step t ≥ 0. The difference between the current iterate P t and the next P t+1 (20)-( 21) is approximated as follows:
The proof of Lemma 6.4 is in Appendix C. It relies on the Taylor expansion of (20).
The following lemma gives properties of α t (P) and shows that when the step size is sufficiently small, grad f (P t )α t (P) † and grad f (P t ) are on the same order. Lemma 6.5. Let t ≥ 0 and P ∈ Gr(n, r). The pseudoinverse α t (P) † ∈ S 2 (R n ) satisfies α t (P) † α t (P) = P. Moreover, there exists a constant Γ 1 (X X X, r) > 0 such that if
Please refer to Appendix D for the proof of Lemma 6.5. It utilizes the fact that α t is a bounded continuous function of P.
Proof Theorem 6.1. With the above useful lemmas and proposition in place, in this proof, we show that our update scheme for P ∈ Gr(n, r) (20)-( 21) satisfies the conditions C1 -C4 in Theorem A.1. Fix a step t ≥ 0. Assume 0 < ε t ∞ Γ 1 where Γ 1 is as in Lemma 6.5.
(C1) We apply (30) that O grad f (P t )α t (P t ) † has the same scale as O ( grad f (P t ) ). By Lemma 6.4, it follows
It implies that there exists a constant C 1 (X X X, r) > 0 such that
By (59) in Appendix C, grad f (P t )α t (P t ) † = (I n -P t )w(P t )P t α t (P t ) † = (I n -P t )w(P t )α t (P t ) † , and so grad f (P t )α t (P t ) † , α t (P t ) † grad f (P t ) = trace (I n -P t )w(P t )α t (P t ) † (I n -P t )w(P t )α t (P t ) † = 0. It implies sym grad f (P t )α t (P t ) † = ∈ 2 2 grad f (P t )α t (P t ) † .
By triangle inequality,
Here, the second inequality is due to the first half of (30). We further constrain the step size as follows
C 1 sup P∈Gr(n,r) grad f (P) ,
so that
Note that sup P∈Gr(n,r) grad f (P) is finite because grad f is continuous on Gr(n, r) and the Grassmann is compact. (We also exclude the trivial case when grad f (P) ≡ 0.) Continuing from (32) and substituting (34), we have
)ε t grad f (P t ) = 2 3 ε t grad f (P t ) .
Letting κ = inf t→∞ 2 3 ε t > 0 ensures that C1 (46) holds. In preparation for C2, we derive an upper bound for P t+1 -P t . The idea is similar to (32):
The second inequality is obtained from the latter half of ( 30) and (34).
(C2) Keep the assumption (33) ε t . Consider the first-order Taylor's expansion for f on Gr(n, r) [7, section 10.7]. Inserting (31), the di"erence of the current and next f is:
grad f (P t ), P t+1 -P t ⇔ + O P t+1 -P t 2 = grad f (P t ), 4ε t sym grad f (P t )α t (P t ) † + O ε 2 t max grad f (P t ) 3 , grad f (P t ) 2 . (36) Since grad f (P t ) is symmetric, the first term in the result of (36) is grad f (P t ), 4ε t sym grad f (P t )α t (P t ) † = 4ε t grad f (P t ), grad f (P t )α t (P t ) † = 4ε t grad f (P t )P t , grad f (P t )α t (P t ) † P t = 4ε t grad f (P t )α t (P t ) † α t (P t ), grad f (P t )α t (P t ) † P t . The second equation is due to α t (P t ) † P t = α t (P t ) † . The definition (28) gives α t (P t ) = P t (I n + 2ε t w(P t ))P t . Continuing from above, grad f (P t ), 4ε t sym grad f (P t )α t (P t ) † = 4ε t grad f (P t )α t (P t ) † P t (I n + 2w(P t )), grad f (P t )α t (P t ) † P t ≥ 4ε t grad f (P t )α t (P t ) † P t 2 = 4ε t grad f (P t )α t (P t ) † 2 ≥ 8 9 ε t grad f (P t ) 2 .
(
Here we used that I n + 2ε t w(P t ) I n for ε t ≥ 0, since w(P t ) ∈ S 2 (R n ) is PSD. The second inequality is derived from (30).
Considering the remainder term in (36), there is a constant C 2 (X X X, r) > 0 such that
If the step size ε t is small enough to satisfy
where Γ 2 > 0 is defined in (33), then
Inserting into (38),
The final inequality is by (35). Thus C2 holds if we set σ = 1
Proof of Theorem 6.2. Fix step size 0 < ε ∞ Γ * * . From Theorem 6.1, the sequence {P t } always converges to a first-order critical point of f (23), thus a fixed point of φ (21) by C3 in Section 6.2. In this proof, we consider the case when the limit P ∈ Gr(n, r) is a "bad" critical point, meaning it is a first-order critical point but fails the second-order criticality (25). The idea is to apply Theorem A.2 to show that the set of Q 0 that converge to bad critical points has zero measure in St(n, r).
Because φ is a local di"eomorphism by Lemma 6.8, we are eligible to apply Theorem A.2. It gives an open neighborhood B P ⊆ Gr(n, r) around each bad critical point P satisfying the conditions in Theorem A.2. As Gr(n, r) is second-countable, it has the Lindelöf property, i.e. every open cover has a countable subcover. So, there exists a countable set C ⊆ Gr(n, r) consisting of some of the bad critical points such that P∈Gr(n,r) is a bad critical point
Define the subset
We claim that P 0 ∈ Gr(n, r) : φ t (P 0 ) converges to a bad critical point as t ⇐ ∀ ⊆ J .
To see this, assume lim t→∞ φ t (P 0 ) = P is bad. By (43), P ∈ B P for some P ∈ C. Since B P is open, there is t 0 ≥ 0 such that φ t (P 0 ) ∈ B P for all t ≥ t 0 . As P is a fixed point of φ, by (50) in Theorem A.2, we obtain φ t0 (P 0 ) ∈ W P . In other words, P 0 ∈ φ -t0 (W P ). So, P 0 ∈ J as claimed.
We claim that J is measure-zero in Gr(n, r). For each P ∈ C, at least one eigenvalue of Hess f (P) is greater than 0. By Lemma 6.7, we know at least one eigenvalue of Dφ(P) is greater than 1. By the first result of Theorem A.2, the dimension of W P is strictly less than the dimension of Gr(n, r). So W P is measure-zero in Gr(n, r) by [40,Corollary 6.12]. Given that the step size satisfies 0 < ε ∞ Γ * * , the map φ is a local di"eomorphism by Lemma 6.8, and so is φ t for all t ≥ 0. Thus their pre-images of measure-zero sets have zero measure. As (φ t ) -1 ≡ φ -t , the set φ -t (W P ) has zero measure for each t ≥ 0. Then J is a countable union of measure-zero sets, so has zero measure.
To conclude, define I ⊆ St(n, r) as the pre-image of J ⊆ Gr(n, r) under the map Q ⇐ QQ . Since the map is a submersion from the Stiefel manifold onto the Grassmannian, I has zero measure in St(n, r) as J does in Gr(n, r). For all initializations Q 0 in St(n, r) but outside the null set I, the sequence {φ t (P 0 )} = {P t } converges to a second-order critical point of f if 0 < ε ∞ Γ * * . The proof of Theorem 6.2 is complete.
In this paper, we developed the PGD method (Algorithm 1) for symmetric Tucker tensor decomposition. It has a simple update scheme compared to other iterative solvers running on Riemannian manifolds. We further designed the SPGD method (Algorithm 2) and SHOEVD solver (Algorithm 3) to decompose the sample moment tensor (1). The algorithms are free of constructing and storing moment tensors, and their iterations are conducted by streaming only a small subset of samples. Tremendous computational savings were seen in numerical experiments with minimal loss of accuracy. We also showed the applicability of the Tucker decomposition of moment tensors in comparison to the CP format, as well as to detect anomalies and allocate assets on real datasets. Finally, utilizing Riemannian manifold optimization theory, we derived theoretical guarantees for the PGD sequence on the Grassmannian. The sequence always converges to a first-order critical point, and almost surely to a second-order critical point.
There are a number of directions that would be interesting to study in the future.
1. Implementation for cumulants. Cumulant tensors are central to non-Gaussian data analysis. They are nonlinear combinations of moments. It would be useful to derive scalable implementations for sample cumulants.
2. Extending the analysis. It is certainly worth investigating convergence guarantees in the streaming setting with adaptive step sizes of Algorithms 2 and 3. We also want to explore under what conditions is the limiting point guaranteed to be a global optimizer.
The gradient of f (22) is ∇f (P) = d X X X • (I n , P, . . . , P) , X X X • (P, . . . , P)⇔ -1 + d X X X • (P, . . . , P) , X X X • (I n , P, . . . , P)⇔ -1 = w(P)P + Pw(P).
Next, we relate the gradients of f (P) and F (Q). The Euclidean gradient of F is
Therefore, I n -QQ ∇F (Q)Q = 2 (I n -P) w(P)QQ = 2 (I n -P) w(P)P.
By (52), we have grad f
The proof of Proposition 6.3 is complete.
C Proof of Lemma 6.4
Proof. We derive the update on P t (21) based on the iterate Q t (20). We decompose the update direction ε t ∇F (Q t ) into two parts respectively in colspan(Q t ) and colspan(Q
denotes the QR factorization. (Note this quantity is full-rank according to C3 in Section 6.2 and orth : St(n, r) ⇐ St(n, r) is C ∞ at full-rank inputs [66]. ) By the derivative formula of the QR factorization given in [68, Algorithm 1] (originally [66]),
To see the second equality, I n -Q t Q t = I n -Q t Q t , because Q t and Q t are both the orthonormal bases for colspan(Q t ). Moreover, given the linear operator s : R r→r ⇐ R r→r , 5 the term
Here the second line refers to the symmetric term
. We transform the right-hand side of ( 56) from an expression in terms of Q to one with P.
For the first component of ( 56), we separate
Also,
by the fact that pseudoinversion satisfies X † Y † = (YX) † if X and Y respectively have orthonormal rows and columns. Inserting (54) gives
where the last equality is by (53) and the definition of α t (28). Moreover,
The formula of grad f (52) is used in the last equality. Therefore, the first component of (56) reads
As for second component of (56), note
= 4 grad f (P t )α t (P t ) † 2 . 5 The operator s is a composition of element-wise multiplication and skew-symmetrization. See [68, Algorithm 1] for details.
The of Q t , which preserves norms, gives the first equality. The second is from (59). The third component of ( 56) is dealt with directly by (27),
Thus, (56) consists of the main component (60) and the residuals O ε 2 t grad f (P t )α t (P t ) † 2 , O ε 2 t grad f (P t ) 2 as needed in (29). This completes the proof of Lemma 6.4.
Proof. Consider the step size ε t > 0 and projector P ∈ Gr(n, r). Definition (28) reads α t (P) = P + 2ε t Pw(P)P = P (I n + 2ε t w(P)) P.
(
So colspan(α t (P)) ⊆ colspan(P). Also for any vector x ∈ colspan(P),
x α t (P)x = x P(I n + 2ε t w(P))Px = x (I n + 2ε t w(P))x ≥ x 2 , as Px = x and I n + 2ε t w(P) I n . It follows colspan(α t (P)) = colspan(P), and α t (P) has only r non-zero singular values, each greater than 1. By properties of pseudoinverses, α t (P) † α t (P) = Proj range(αt(P)) = P,
and α t (P) † has only r non-zero singular values, each less than 1. Hence
Constrain the step size by 0 < ε t ∞ Γ 1 := 1 4 ∈ r sup P∈Gr(n,r) Pw(P)P ,
where sup P∈Gr(n,r) Pw(P)P is finite because the function P ⇐ Pw(P)P is continuous and Gr(n, r) is compact. (We exclude the trivial case Pw(P)P ≡ 0 when X X X = 0.) From ( 62), Pα t (P) † = α t (P) † α t (P)α t (P) † P = -α t (P) † (Pα t (P)).
Then by ( 63) and ( 64), Pα t (P) † =α t (P) † (Pα t (P)) ∞ α t (P) † Pα t (P)
We now compare grad f (P) and grad f (P)α t (P) † . By (27), grad f (P)P = (I n -P)w(P)P.
Inserting the above equation into (27), grad f (P) = sym 2(I n -P)w(P)P = sym 2 grad f (P)P .
Then due to triangle inequality and (65), grad f (P) ∞ ∈ 2 grad f (P)α t (P) † + ∈ 2 grad f (P) Pα t (P) † = ∈ 2 grad f (P)α t (P) † + ∈ 2 grad f (P)α t (P) † Pα t (P) † ∞ ∈ 2 grad f (P)α t (P) † 1 + Pα t (P)
grad f (P)α t (P) † .
Inserting (66) shows the last inequality. Thus, grad f (P)α t (P) † ≥ ∈ 2 3 grad f (P) .
On the other hand, as α t (P) † = Pα t (P) † , it holds grad f (P)α t (P) † ∞ grad f (P)P α t (P) † ∞ ∈ 2r 2 grad f (P) .
The results ( 67) and ( 63) give the last derivation. Combining the above two inequalities completes the proof of Lemma 6.5.
E Proof of Proposition 6.6
Proof. Let P ∈ Gr(n, r) and %P ∈ T P Gr(n, r). We first show the computation of the Riemannian Hessian. Applying
Recall (51) that ∇f (P) = w(P)P + Pw(P). We di"erentiate ∇f (P) to obtain the Euclidean Hessian of f :
∇
where the di"erential is
= 2d(d -1) sym X X X • (I n , %P, P, . . . , P) , X X X • (I n , P, P, . . . , P)⇔ -1 = d(d -1) X X X • (I n , (%P)P + P%P, P, . . . , P) , X X X • (I n , I n , P, . . . , P)⇔ -1 = d(d -1) X X X • (I n , %P, P, . . . , P) , X X X • (I n , I n , P, . . . , P)⇔ -1 .
(70)
The last equality follows from (%P)P + P%P = %P in (24). (71)
We compute the matrix M. Due to the rank constraint on P, we utilize its block structure to facilitate the computation. W.l.o.g., we suppose the rank-r orthogonal projection P is: P = I r 0 r→(n-r) 0 (n-r)→r 0 (n-r)→(n-r)
∈ Gr(n, r), by choosing an appropriate for R n . From (52), grad f (P) = sym 2(I n -P)w(P)P = 0, and ( 24), %P = (%P) P + P%P for all %P ∈ T P Gr(n, r), we obtain block matrices for w(P) and {%P i }. Namely there are PSD matrices A ∈ S 2 (R r ), C ∈ S 2 (R n-r ) and matrices H i ∈ R (n-r)→r for i ∈ [r(nr)] such that
Since %P i are eigenmatrices of Hess f (P), by (41) we see Here we note that since A is PSD and ε > 0, the inverse (I r + 2εA) -1 exists and is symmetric. Thus the left component in the inner product (80) is 2 sym
For ease of analysis, we ensure the uniqueness of the QR factorization by using the one where all diagonal entries of the R factor are positive. But none of our algorithms require this choice.
Similar to our algorithm scheme, the comparable algorithms are mostly gradient and Hessian-driven methods. Since we could not identify an implementation for the Jacobi rotation method[31,41], it is not included in the comparison.
The dataset is available at http://xudongkang.weebly.com/data-sets.html.
The HFRX indices dataset can be found at: https://www.hfr.com/family-indices/hfrx .
√ r in(47). (C3) "⇒": From grad f (P) = 0 and (27), we knowgrad f (P)Q = 0 = grad F (Q) = (I n -QQ )∇F (Q),