This work considers the problem of computing the canonical polyadicdecomposition (CPD) of large tensors. Prior works mostly leverage data sparsityto handle this problem, which is not suitable for handling dense tensors thatoften arise in applications such as medical imaging, computer vision, andremote sensing. Stochastic optimization is known for its low memory cost andper-iteration complexity when handling dense data. However, exisitingstochastic CPD algorithms are not flexible enough to incorporate a variety ofconstraints/regularizations that are of interest in signal and data analytics.Convergence properties of many such algorithms are also unclear. In this work,we propose a stochastic optimization framework for large-scale CPD withconstraints/regularizations. The framework works under a doubly randomizedfashion, and can be regarded as a judicious combination of randomized blockcoordinate descent (BCD) and stochastic proximal gradient (SPG). The algorithmenjoys lightweight updates and small memory footprint. In addition, thisframework entails considerable flexibility---many frequently used regularizersand constraints can be readily handled under the proposed scheme. The approachis also supported by convergence analysis. Numerical results on large-scaledense tensors are employed to showcase the effectiveness of the proposedapproach.
Quick Read (beta)
Block-Randomized Stochastic Proximal Gradient for Low-Rank Tensor Factorization
This work considers the problem of computing the canonical polyadic decomposition (CPD) of large tensors. Prior works leverage data sparsity to handle this problem, which is not suitable for handling dense tensors that often arise in applications such as medical imaging, computer vision, and remote sensing. Stochastic optimization is known for its low memory cost and per-iteration complexity when handling dense data. However, existing stochastic CPD algorithms are not flexible to incorporate a variety of constraints/regularization terms that are of interest in signal and data analytics. Convergence properties of many such algorithms are also unclear. In this work, we propose a stochastic optimization framework for large-scale CPD with constraints/regularization terms. The framework works under a doubly randomized fashion, and can be regarded as a judicious combination of randomized block coordinate descent (BCD) and stochastic proximal gradient (SPG). The algorithm enjoys lightweight updates and small memory footprint. This framework entails considerable flexibility—many frequently used regularizers and constraints can be readily handled. The approach is supported by convergence analysis. Numerical results on large-scale dense tensors are presented to showcase the effectiveness of the proposed approach.
Canonical polyadic decomposition (CPD) [previously known as parallel factor analysis (PARAFAC)] [1, 2, 3] is arguably the most popular low-rank tensor decomposition model. CPD has found applications in various fields, such as analytical chemistry , social network mining , hyperspectral imaging , topic modeling , and time series analysis ; also see [9, 10, 11] for more applications in communications.
Computing the CPD of a tensor, however, is a challenging optimization problem . Many algorithms have been proposed through the years [13, 14, 15, 3]. To keep pace with the ever growing volume of available data, one pressing challenge is to compute CPD at scale. The classic alternating least squares (ALS) algorithm  has an elegant algorithmic structure, but suffers from a number of numerical issues [16, 17] and is hardly scalable. In recent years, many new CPD algorithms have appeared, triggered by the advances in big data analytics and first-order optimization [18, 19, 14, 13, 20]. Many of these algorithms leverage data sparsity to scale up CPD—by cleverly using the zero elements in huge tensors, computationally costly key operations in ALS (e.g., the matricized tensor times Khatri-Rao product (MTTKRP) operation) can be significantly simplified. Consequently, the classic ALS algorithm can be modified to handle CPD of huge and sparse tensors.
However, when the tensor to be factored is dense—i.e., when most entries of the tensor are nonzero—the sparsity-enabled efficient algorithms [13, 14, 18, 19, 6] are no longer applicable. Note that large and dense tensors arise in many timely and important applications such as medical imaging , hyperspectral imaging , and computer vision . In fact, since big dense tensors typically cost a lot of memory (e.g., a dense tensor with a size of occupies 57.52GB memory if saved as double-precision numbers), it is even hard to load them into the RAM of laptops, desktops, or servers.
Stochastic gradient (SG) method is a powerful tool for handling optimization problems involving a large amount of data, which is known for its low per-iteration memory and computational complexities . A number of stochastic optimization based CPD algorithms have been proposed in the literature [24, 25, 26]. Specifically, The works in [24, 25] work in an iterative manner. In each iteration, the algorithm samples a random subset of the tensor entries and update the corresponding parts of the latent factors using the sampled data. The algorithms have proven quite effective in practice, and feature distributed implementation . The challenge here is that every tensor entry only contains information of a certain row of the latent factors, and updating the entire latent factors may require a high complexity. This may lead to slow improvement of the latent factor estimation accuracy. More importantly, this update strategy loses the opportunity to incorporate constraints/regularization terms on the whole latent factors, since the sampled entries only contain partial information of them. This is undesired in practice, since prior information on the latent factors are critical for enhancing performance, especially in noisy cases.
Recently, a stochastic algorithm that ensures updating one entire latent factor in every iteration was proposed in . Instead of sampling tensor entries, the algorithm works via sampling tensor fibers that contain information of the whole latent factors. This algorithm exhibits very good empirical performance when the tensor rank is low. However, this algorithm works with at least as many fibers as the tensor rank, which in some cases gives rise to high per-iteration complexity. In addition, the algorithm in  did not explicitly offer implementations that take into considerations of constraints or regularization terms—although this can be fixed with some modifications. Lastly, convergence properties of many stochastic CPD algorithms such as those in [24, 26] are unclear.
Contributions In this work, we propose a new stochastic algorithmic framework for computing the CPD of large-scale dense tensors. Our contributions include:
A Doubly Randomized Computational Framework for Large-Scale CPD. We propose an efficient and flexible computational framework for CPD of large dense tensors. Our method is a judicious combination of randomized block coordinate descent (BCD) [27, 28] and stochastic proximal gradient (SPG) [29, 30]. In each iteration, our algorithm first samples a mode from all modes of the tensor. Then, the algorithm samples some fibers of this mode and updates the corresponding latent factor via stochastic proximal operations. Such a combination exhibits an array of attractive features: It admits smaller per-iteration memory and computational complexities relative to the existing fiber sampling based method in , particularly in high-rank cases. More importantly, it is very flexible in terms of incorporating regularization terms and constraints on the latent factors.
Convergence Analysis. Both BCD and SPG are well studied topics in the optimization literature [31, 28, 27]. However, convergence properties of the proposed framework is not immediately clear due to the nonconvex nature of CPD. The existing block-randomized SGD (BR-SGD) framework  only considers convex optimization. Another work  considers nonconvex optimization, which adopts a Gauss-Seidel type block updating strategy without randomization. The conditions for convergence in  are not easy to check or guarantee in the context of tensor factorization. In contrast, we offer tailored convergence analyses of the proposed algorithm leveraging block randomization, and show that the proposed optimization strategy features sub-sequence convergence to a stationary point—which is a necessary condition for attaining local or global optima.
Implementation-friendly Adaptive Stepsize Scheduling. One of the most challenging aspects in stochastic optimization is selecting a proper stepsize schedule. To make the proposed algorithms friendly to use by practitioners, we propose a practical and adaptive stepsize schedule based on the celebrated Adagrad algorithm . Adagrad is an adaptive stepsize selection method devised for single-block gradient descent. Nonetheless, we find through extensive simulations that it largely helps reducing the agonizing pain of tuning stepsize when implementing our multi-block algorithm for CPD. In addition, we show that the adaptive stepsize-based algorithm converges to a stationary point almost surely.
A quick demonstration of the effectiveness of the proposed algorithms is shown in Fig. 1, where the average mean squared error (MSE) of the estimated latent factors [cf. Eq. (22)] against the number of MTTKRP computed (which dominates the complexity) is plotted. One can see that the proposed algorithm largely outperforms a couple of state-of–the-art algorithms for constrained CPD. More thorough numerical results can be seen in Sec. VI.
Part of the work was submitted to ICASSP 2019 . In this new version, we have included detailed convergence proofs and the new adaptive stepsize based algorithm. More extensive simulations and real-data experiments are also included.
Notation. We follow the established conventions in signal processing. , , , and denote scalar, vector, matrix, and tensor, respectively; denotes the Euclidean norm, i.e., and , respectively; , , and denote outer product, Khatri-Rao product, and Hadamard product, respectively, unless otherwise specified; denotes the vectorization operator that concatenates the columns of ; means that all the entries of are nonnegative; and denote transpose and pseudo-inverse, respectively; denotes the cardinality of set ; denotes the largest eigenvalue of a matrix. Unless otherwise specified, we denote the total expectation by the subscript-less operator .
We first introduce some notions used in tensor algebra.
II-A Tensors and CPD
An th order tensor is an array whose entries are indexed by coordinates; i.e., denotes an element of the tensor with a size of . Like matrices, tensors can be represented as sum of rank-one components:
where “” denotes the outer product of vectors, and is an matrix that is often referred to as the mode- latent factor. When is the minimal integer that satisfies the expression in (1), the right hand side in (1) is called the canonical polyadic decomposition of the tensor . At the entry level, the CPD can be expressed as
for . The CPD of a tensor is essentially unique under mild conditions11 1 The latent factors ’s that constitute the data are unique up to some trivial ambiguities like column permutations and scalings .. The CPD of a tensor can be obtained via minimizing a certain criterion as follows:
In the sequel, we will often use the shorthand notation to denote , where
II-B Unfolding, ALS and MTTKRP
The matricization operation, or matrix unfolding of a tensor, has proven very useful in designing tensor factorization algorithms. The mode- unfolding of a tensor is a matrix where
where the matrix is defined as
The elegant form of the unfoldings has enabled the famous alternating least squares (ALS) algorithm  for handling Problem (3) with the LS objective. Specifically, ALS solves the following cyclically for :
Problem (5) is nothing but a least squares problem that admits the following closed-form solution:
if . Note that is not difficult to compute by exploiting the Khatri-Rao structure of [13, 2, 1]. However, when the problem dimension is large (which often happens in applications such as medical imaging, remote sensing, and computer vision), solving the seemingly simple problem in (5) can be computationally prohibitive. The reason is that both and can be very large matrices. In particular, the so-called matricized tensor times Khatri-Rao product (MTTKRP) operation, i.e.,
that happens in every iteration of ALS costs flops (or, if ). This is quite costly even if is moderately large. Many works have considered fast algorithms for computing MTTKRP, but these methods are mainly enabled by judiciously exploiting sparsity of the tensor data [39, 18]. Computing MTTKRP for dense tensors has also been considered. Nonetheless, these works are often concerned with practical implementation schemes such as parallelization and memory-efficient computation strategies, but the number of computational flops required is naturally high for the dense tensor case; see, e.g., [40, 41].
In a lot of applications, some prior knowledge on the latent factors is known—e.g., in image processing, ’s are normally assumed to be nonnegative ; in statistical machine learning, sometimes the columns of are assumed to be constrained within the probability simplex [36, 42]; i.e.,
In those cases, the following criterion is often of interest:
Compared to the unconstrained version, Problem (7) is even harder to handle. Some recent methods combine first-order constrained optimization and ALS [14, 13] to make the tensor factorization algorithms more flexible in handling constraints and regularization terms—but the complexity orders of those algorithms often scale similarly as that of ALS, since these algorithms do not avoid computing that is the bottleneck for computing CPD.
II-C Stochastic Optimization
When the tensor is large and dense, working with the entire dataset could be computationally and memory-wise expensive. A popular workaround is to apply stochastic optimization—i.e., sampling parts of the data at random and use the sampled piece to update the latent factors. Using Eq. (2), Problem (3) with the LS objective is equivalent to the following:
The objective function in (8) can be understood as an empirical risk . Using this observation, the algorithms in [24, 25] randomly sample a subset of entries indexed by and update the pertinent parts of the latent factors (note that the th entry of tensor contains the information of for ) using the sampled entries of the tensor. For example,  uses a stochastic gradient (SG) based approach and update the ’s that are associated with the sampled entries. The sampling method in  is similar, while the update is not gradient-based but Gauss-Newton or ALS applied to the sampled set of entries (or, sub-tensors, to be precise). The upshot of this line of work is that the per-iteration complexity can be quite low.
Despite of such favorable complexity savings, the approaches in [24, 25] have a couple of limitations. One challenge is that many useful prior information cannot be incorporated in the algorithm. The reason is that these algorithms update part of the rows of ’s, while many useful priors are defined w.r.t. the columns of the latent factors, e.g., the probability simplex constraint in (6) and the total variation constraint that is heavily used in image processing. For example, the algorithm in  samples subtensors (where ) to update the corresponding ’s. Under such a scheme, it is hard to handle constraints like that are critical in statistical learning [42, 43, 44, 45]. Third, convergence properties of these methods are often unclear.
An alternative  is to leverage the tensor data structure by considering randomly sampled fibers of tensors. Note that a mode- fiber of (cf. Fig. 2) is a row of the mode- unfolding . Now, assuming that one samples a set of mode- fibers indexed by , then can be updated by solving a ‘sketched version’ of Problem (5):
If , then the sketched system of linear equations can be over-determined. Hence, one can update by solving the dimensional linear system
Similar to the ALS algorithm, after updating , the algorithm moves to mode- fibers and repeats the same for updating . The downside of this method is that it needs to sample at least fibers for each update, and can be larger than in tensor decomposition. In addition, the work in  focused on unconstrained cases while did not offer implementations for constrained/regularized cases. This can be compensated by replacing (10) with a constrained least squares solver or certain constraint enforcing operations. However, convergence properties of doing so are also unclear.
III Proposed Algorithm
In this work, we propose a new stochastic optimization strategy for CPD. Our method combines the insights from ALS and fiber sampling, but allows . This is instrumental in practice, since it is the key for achieving low per-iteration complexity. The proposed algorithm can easily handle a variety of constraints and regularizations that are commonly used in signal processing and data analytics—which is reminiscent of stochastic proximal gradient (SPG) [30, 46]. In addition, we provide convergence analyses to back up the proposed approach.
III-A Basic Idea: Unconstrained Case
We first consider Problem (3). Our idea is to apply SA while exploiting the tensor fiber structure. Specifically, at each iteration, we sample a set of mode- fibers for a certain as the method in  does. However, instead of exactly solving the least squares subproblems (5) for all the modes following a Gauss-Seidel manner in each iteration, we update using a doubly stochastic procedure. To be more precise, at iteration , we first randomly sample a mode index . Then, we randomly sample a set of mode- fibers that is indexed by . Let such that
where we have
and we used the shorthand notations
The latent variables are updated by
Observe that is an estimate of the gradient applied to taken w.r.t. the mode- variable , and the update is an iteration of the SG algorithm with a minibatch size for solving the problem in (5).
The proposed update is very efficient, since the most resource-consuming update in algorithms such as those in [13, 14] is avoided. The corresponding part costs only flops—and is under our control. Note that the first step in this procedure is different from standard ALS-type algorithms that update the block variables cyclically instead of updating a randomly sampled block. As we will show, this modification greatly simplifies our convergence analysis.
III-B Constrained and Regularized Case
As mentioned, there are many cases in practice where considering regularizations or constraints on ’s can benefit the associated tasks. Since our framework updates an entire in each iteration, it is friendly for incorporating a large variety of commonly used constraints/regularizations—which is more flexible relative to the entry sampling based approaches in [24, 25]. The algorithm can be easily extended to handle the constrained/regularized case. Consider:
where denotes a structure-promoting regularizer on . Note that if is desired, we can write is defined as the indicator function of set :
Using the same fiber sampling strategy as in the previous subsection, we update by
If is a closed proper convex function, the update (14a) can be solved by applying the proximal operator of , which is often denoted as
Many ’s admit simple closed-form solutions for their respective proximal operators, e.g., when is the indicator function of the nonnegative orthant and ; see Table I and more details in [13, 47]. The complexity of computing (15) is similar to that of the plain update in (12), and thus is also computationally efficient. An overview of the proposed algorithm can be found in Algorithm 1, which we name Block-Randomized SGD for CPD (BrasCPD).
|randomized pivot search ||in expectation|
|monotonic||monotone regression |
|unimodal||unimodal regression |
In the table, is the number of optimization variables.
IV Convergence Properties
In this section, we offer tailored convergence analyses for BrasCPD. To this end, two most relevant works from the optimization literature are  and . The work in  considers block-randomized SGD, but only for the convex case—while our problem is nonconvex. The work in  considers the Gauss-Seidel type block SGD (i.e., cyclically updating the blocks), instead of the block-randomized version as BrasCPD uses. There, convergence is established using a number of assumptions that are not easy to check or guarantee in the context of CPD, e.g., that the bias of the stochastic oracle is bounded. We will show that, by using the block-randomization strategy and the proposed stochastic oracle construction, such an assumption can be circumvented22 2 We should note that the major motivation for using block randomization strategy is theoretical guarantees—since our goal is a convergence-guaranteed algorithmic framework; in practice, we observe cyclically updating the latent factors works as well..
To facilitate our discussions, let us define and as the random variables (r.v.) responsible for selecting the mode and fibers in iteration , respectively. These r.v.s are distributed as
where , such that is the set of all subsets of with size . We observe
Denote as the filtration generated by the r.v.s such that the th iterate is determined conditioned on . The stochastic gradient in (11) is an unbiased estimate for the full gradient w.r.t.
The proof of the above is straightforward and thus skipped. Fact 1 says that our block stochastic gradient is an unbiased estimation for the “block gradient” . This fact will prove quite handy in establishing convergence. The two-level sampling strategy (i.e., block sampling and fiber sampling, respectively) makes the gradient estimation w.r.t. unbiased up to a scaling factor (see Appendix A). This connection intuitively suggests that the proposed algorithm should behave similarly as an SG algorithm.
IV-A Unconstrained Case
We will use the following assumptions:
The stepsize schedule follows the Robbins-Monro rule :
The updates are bounded for all .
Assumption 1 is a principle for stepsize scheduling, which is commonly used in SG algorithms. Assumption 2 is considered a relatively strong assumption. In practice, there are several simple ways to make ’s bounded. One pragmatic modification is to change the objective to . Another method is as mentioned in [31, 13]. At iteration , one may add a proximal term to the cost function, which will effectively prevent from being unbounded. Following both ways, the updates can still be handled by simple proximal operations for the functions in Table I.
There are an array of problem structures that are useful for studying convergence of the algorithm.
For any and mode , there exists a constant such that
where and are extracted/constructed from following the respective definitions.
Eq. (2) holds because the objective function w.r.t. is a plain least squares fitting criterion, which is known to have a Lipschitz continuous gradient—and the smallest Lipschitz constant is .
We have the following convergence property for BrasCPD in the unconstrained case:
The proof is relegated to Appendix B. The above proposition implies that there exists a subsequence of the solution sequence that converges to a stationary point in expectation. The use of the expectation notion is due to the randomness in the algorithm. We should mention that the SG update and the block sampling step are essential for establishing convergence—and using the exact solution to (9) as in  may not have such convergence properties.
IV-B Constrained/Regularized Case
When , the gradient of the objective function of (13) may be undefined. In this case, a solution is stationary if , , where
There exists a sequence such that
where is the stepsize sequence following Assumption 1.
The BrasCPD produces a convergent solution sub-sequence:
The convergence result in Proposition 2 inherits a drawback from single-block stochastic proximal gradient algorithms for nonsmooth nonconvex optimization. To be specific, the relatively strong assumption 3 is required to ensure convergence. Assumption 3 implies that the variance of the gradient estimation error converges to zero. This is not entirely trivial. One way to fulfill this assumption is to increase the minibatch size along the iterations, e.g., by setting [30, 33]:
Another popular way for achieving (19) is to use some advanced variance reduction techniques such as SVRG —which may go beyond the scope of this paper and thus is left out of the discussion. Also notice that as the convergence analysis is pessimistic, in practice constant minibatch size works fairly well—as we will see soon.
V An Adaptive Stepsize Scheme
One may have noticed that the convergence theories in Propositions 1-2 do not specify the sequence except for the two constraints in Assumption 1. This often gives rise to agonizing tuning experience for practitioners when implementing stochastic algorithms.
Recently, a series of algorithms were proposed in the machine learning community for adaptive stepsize scheduling when training deep neural networks [52, 53, 54]. Most of these works are variants of the Adagrad algorithm . The insight of Adagrad can be understood as follows: If one optimization variable has been heavily updated before, then it is given a smaller stepsize for the current iteration (and a larger stepsize otherwise). This way, all the optimization variables can be updated in a balanced manner. Adagrad was proposed for single-block algorithms, and this simple strategy admits many provable benefits under the context of convex optimization . For our multi-block nonconvex problem, we extend the idea and propose the following updating rule: In iteration , if , then, for all and all , we have
where . We note that , are technical conditions used for establishing theoretical convergence. In practice, setting does not hurt the performance and we also observe a slight gain in runtime performance when . The Adagrad version of block-randomized CPD algorithm is very simple to implement. The algorithm is summarized in Algorithm 2, which is named AdaCPD.
As one will soon see, such a simple stepsize strategy is very robust to a large number of scenarios under test—i.e., in most of the cases, AdaCPD performs well without tuning the stepsize schedule. In addition, the AdaCPD algorithm works well for both the constrained and unconstrained case.
Proposition 3 asserts that the algorithm converges almost surely. The proof is relegated to Appendix D. Our proof extends the idea from a recent paper  that focuses on using Adagrad for solving single-block nonconvex problems. As mentioned, our two-level sampling strategy makes our algorithm very similar to single-block SGD with a scaled gradient estimation (cf. Appendix A), and thus with careful modifications, the key proof techniques in  goes through. Nevertheless, we detail the proof for being self-containing.
VI Numerical Results
In this section, we use simulations and real-data experiments to showcase the effectiveness of the proposed algorithm.
VI-A Synthetic Data Simulations
VI-A1 Data Generation
Throughout this subsection, we use synthetic third-order tensors (i.e., ) whose latent factors are drawn from i.i.d. uniform distribution between and —unless otherwise specified. This way, large and dense tensors can be created. For simplicity, we set for all and test the algorithms on tensors having different ’s and ’s. In some simulations, we also consider CPD for noisy tensors, i.e., factoring data tensors that have the following signal model:
where is the noiseless low-rank tensor and denotes the additive noise. We use zero-mean i.i.d. Gaussian noise with variance in our simulations, and the signal-to-noise ratio (SNR) (in dB) is defined as .
A number of baseline algorithms are employed as benchmarks. Specifically, we mainly use the AO-ADMM algorithm  and the APG algorithm  as our baselines since they are the most flexible algorithms with the ability of handling many different regularizations and constraints. We also present the results output by the CPRAND algorithm . Note that we are preliminarily interested in constrained/regularized CPD. Because CPRAND operates without constraints, the comparison is not entirely fair (e.g., CPRAND can potentially attain smaller cost values since it has a much larger feasible set if other algorithms operate with constraints). Nevertheless, we employ it as a benchmark since it uses the same fiber sampling strategy as ours. To make the comparisons more comprehensive, we also offer a simple modification for CPRAND to incorporate constraints/regularization terms. Specifically, we apply the proximal operators associated with the constraint/regularization terms to the original CPRAND updates, and we denote this baseline as CPRAND-Prox. All the algorithms are initialized with the same random initialization; i.e., ’s entries follow the uniform distribution between 0 and 1.
VI-A3 Parameter Setting
For BrasCPD, we set
where is the number of iterations, and typically takes a value in between 0.001 and 0.1, and we try multiple choices of in our simulations. Our experience is that, under such settings, is the main tuning parameter that affects the performance of BrasCPD. The batch size is typically set to be below 25 throughout this section. For AdaCPD, we fix , , and for all the simulations. For CPRAND, we follow the instruction in the original paper  and sample fibers for each update.
VI-A4 Performance Metrics
To measure the performance, we employ two metrics. We mainly use the estimation accuracy for the latent factors, for , as our performance indicator. The accuracy is measured by the mean squared error (MSE) which is as defined in [58, 59]:
where denotes the estimate of and ’s are under the constraint —which is used to fix the intrinsic column permutation in CPD. We also use the cost function, i.e., as a reference in some of the simulation tables.
Since the algorithms under test have very different operations and subproblem-solving strategies, it may be challenging to find an exactly unified complexity measure. In this section, we show the peformance of the algorithms against the number of MTTKRP operations used, since is the most costly step that dominates the complexity of all the algorithms under comparison. For the stochastic optimization/sketching based algorithms, we also test the MSE/cost value of the algorithms against runtime and the number of sampled entries for updating the latent factors. The latter is particularly meaningful under the stochastic settings, since it affects the communication overhead between the data storage units (e.g., the hard disks) and the computation units, e.g., CPUs and GPUs. All the simulations are conducted in Matlab. The results in this section are obtained from 50 trials with different randomly generated tensors.
Fig. 1 in Sec. I has shown the MSE performance of the algorithms in a relatively small-size example, where , and the nonnegativity constraints are used in the algorithms. In that simulation, we use so that every 500 iterations of the proposed algorithm compute a full MTTKRP. One can see that for this relatively easy case, all the algorithms can reach a good estimation accuracy for the latent factors. Nevertheless, the proposed methods exhibit remarkably higher efficiency.
Fig. 3 shows the MSEs of the estimated latent factors by the algorithms under and benchmarked by more baselines. The result is the median of 50 Monte Carlo trials; we use median here since mean is dominated by outlying trials, even if there is only one outlying trial. We set so that the proposed algorithms use 5,000 iterations to compute a full MTTKRP. All the algorithms use nonnegativity constraints except CPRAND. One can see that under this setting, the most competitive algorithms are CPRAND, CPRAND-Prox, and AdaCPD.
In Fig. 4, we increase the rank to be . There are several observations in order: First, the stochastic algorithms (i.e., BrasCPD, AdaCPD, CPRAND and CPRAND-Prox) are much more efficient relative to the deterministic algorithms (AO-ADMM and APG). After 60 MTTKRPs computed, the stochastic algorithms often have reached a reasonable level of MSE. This is indeed remarkable, since 60 MTTKRPs are roughly equivalent to 20 iterations of AO-ADMM and APG. Second, CPRAND and CPRAND-Prox are not as competitive in this high-rank regime. In particular, BrasCPD with gives the most promising performance. However, the performance of BrasCPD is affected a bit significantly by the parameters . One can see that using and , the algorithm does not give promising results under this setting. Third, AdaCPD yields the second lowest MSEs after 60 MTTKRPs—while not using any parameter tuning.
Figs 5-6 show the MSEs of the stochastic algorithms (i.e., BrasCPD, AdaCPD, CPRAND and CPRAND-Prox) against the number of sampled entries and runtime, respectively. Here, we set and . One can see that BrasCPD and AdaCPD do not use many data samples to reach a good MSE level. This indicates that the communication overhead performance of the proposed algorithms is promising. In terms of runtime, CPRAND and CPRAND-Prox start with quick decrease of the MSE—they reach MSE within 50 seconds, while the proposed approaches reach this level of accuracy after 100 seconds. Nevertheless, the MSEs of CPRAND and CPRADN-Prox are somehow stuck at this level, but the proposed algorithms can reach a much better accuracy for estimating the latent factors. We would like to remark that the runtime performance of stochastic algorithms are affected by the programming language used (i.e., Matlab in this case). Typically, interpreted languages (e.g., Matlab and Python) are not specialized for handling “for” loops, which is heavily used in stochastic algorithms (especially when is small). Hence, real-system implementations for these algorithms could be much faster.
Table II shows the mean and median of the MSEs and cost values output by the algorithms when the tensor rank varies under . All the algorithms are stopped after 60 full MTTKRPs are used. One can see that BrasCPD exhibits a quite competitive MSE performance if a proper is chosen, under the employed stepsize schedule in (21). However, one can see that when changes, there is a risk that BrasCPD runs into numerical issues and yields unbounded solutions. This suggests that BrasCPD may need extra care for tuning its stepsize. In principle, when the problem setting changes, the “best” of BrasCPD also changes. Our experience is that when increases, using a properly scaled up may help accelerate convergence. On the other hand, AdaCPD always outputs reasonably good results. More importantly, AdaCPD runs without tuning the stepsize parameters—which shows the power of the adaptive stepsize scheduling strategy. Another remark is that CPRAND and CPRAND-Prox work well when and . In particular, CPRAND-Prox largely outperforms CPRAND when , showing the benefit of explicitly considering constraints.
Table III shows the performance under different ’s when . In general, when increases, the performance of all the algorithms improves—with a fixed , a larger means more data and more “degrees of freedom” available, which normally leads to better performance. Again, BrasCPD with a proper and AdaCPD in general outperform the baselines. One particular observation is that, although the mean and median MSEs of AdaCPD are both low, the median is sometimes much better than the mean (cf. the case when ), which indicates that there exist outlying trials (i.e., trials where AdaCPD does not produce very low MSE results). The median-mean gap is less often observed for other algorithms, e.g., BrasCPD and CPRAND. Fig. 7 may better illustrate the situation, where the histograms of the MSE (in dB) of the algorithms are shown. One can see that, although the worst-case result of AdaCPD is still acceptable (with MSE), the MSE of AdaCPD clearly has a larger variance compared to BrasCPD with . This shows a trade-off between the easiness of stepsize scheduling and the risk of converging to less accurate solutions.
Table IV shows the median of the MSEs after 450 MTTKRPs of under same simulation setting (which is roughly equivalent to 150 iterations of AO-ADMM and APG). One can see that the performance of all the algorithms have improved. Note that due to resource limitations, we could not run 450 MTTKRPs for every table in this section and had to stop at 60 MTTKRPs. Hence, the capabilities of all the algorithms may not have fully shown up in those tables. Nonetheless, the 60-MTTKRP results can serve as a reference under limited time and computational resources.
Tables V-VI show the estimation accuracy of the latent factors by the algorithms under different SNRs. Here, we use a low-rank case where and , under which CPRAND and CPRAND-Prox are more competitive. In a noisy environment, the ability of handling constraints/regularizations is essential for a CPD algorithm, since prior information on the latent factors can help improve estimation accuracy. Table V and Table VI test the cases where is elementwise nonnegative and the columns of reside in a scaled version of the probability simplex, respectively. One can see from the two tables that both BrasCPD (with a proper ) and AdaCPD work very well. In Table VI, one can see that BrasCPD again shows its sensitivity to the choice of , with and actually not working. One can also see that, in the low-SNR regime (SNR= and dB), BrasCPD with proper and AdaCPD outperform the baselines. CPRAND-Prox also works well for the nonnegativity constraint case (especially for SNR= and 40dB), but not as promising for the simplex constraints.
VI-C Real-Data Experiment
In this subsection, we test our algorithm on a constrained tensor decomposition problem; i.e., we apply the proposed BrasCPD and AdaCPD to factor hyperspectral images. Hyperspectral images (HSIs) are special images with pixels measured at a large number of wavelengths. Hence, an HSI is usually stored as a third-order tensor with two spatial coordinates and one spectral coordinate. HSIs are dense tensors and thus are suitable for testing the proposed algorithms. We use sub-images of the Indian Pines dataset that has a size of and the Pavia University dataset33 3 Both datasets are available online: http://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes that has a size of . We apply the stochastic algorithms to the datasets by fixing in this section.
Tables VII-VIII show the cost values of the nonnegativity constrained optimization algorithms under different ranks, after computing 120 MTTKRPs for all three modes, which corresponds to 120 iterations for AO-ADMM and APG (we use this “all-mode MTTKRP” in this section since the tensors are unsymmetrical and thus single-mode MTTKRPs cannot be directly translated to iterations in batch algorithms). One can see that the proposed algorithms show the same merits as we have seen in the simulations: BrasCPD can exhibit very competitive performance when is properly chosen (e.g., when and for the Indian Pines dataset); in addition, AdaCPD gives consistently good performance without tuning the stepsize manually. Particularly, on the Pavia University dataset, AdaCPD gives much lower cost values compared to other algorithms. We also present the results of CPRAND and CPRAND-Prox. Note that since CPRAND does not have constraints, its cost value is naturally lower than the other methods. Hence, this baseline is only for reference. Both CPRAND and CPRAND-Prox work reasonably well for the two datasets, especially on Indian Pines under low rank. However, we note that for Pavia University, CPRAND-Prox does not converge well for the case.
Fig. 8 shows how the cost values change along with the iterations on the Pavia University data using . One can see that BrasCPD () and AdaCPD reduce the cost value quickly in this case. After 120 iterations (equivalent to 120 all-mode full MTTKRPs), the batch algorithm APG eventually reaches the same cost value level of those of BrasCPD () and AdaCPD.
To conclude, we proposed a block-randomized stochastic proximal gradient based CPD algorithmic framework for large-scale dense tensors. The framework works under a doubly stochastic manner, which randomly selects a mode and then samples a set of fibers for updating the associated latent factor. The framework has a series of nice features including being able to quickly improve estimation accuracy of the latent factors, being flexible with incorporating constraints and regularizations, and having rigorous convergence guarantees. We also proposed a practical and effective adaptive stepsize scheduling method that is reminiscent of recent advances in neural network training algorithms. Simulations and real-data experiments show that the proposed algorithms outperform a number of state-of-art constrained CPD algorithms when dealing with large dense tensors.
Appendix A Connection between and
Let be any integer from , consider the following conditional expectation:
where is the Dirac function. In the above chain, (a) is due to Fact 1 and (b) is obtained by evaluating the expectation with respect to the possible modes . The last equality shows that is a scaled gradient of the objective function of (3) taken w.r.t. . The block sampling step together with fiber sampling entails us an easy way to estimate the full gradient w.r.t. all the latent factors in an unbiased manner.
Appendix B Proof of Proposition 1
Let and be two nonnegative sequences such that is bounded, converges and diverges, then we have
where we have denoted
Taking expectation conditioned on the filtration and the chosen mode index , we have
As is bounded for all , and all the are bounded under Assumption 2, we have for all and some . Taking the expectation w.r.t. yields
Note that . Taking the total expectation (w.r.t. all random variables in ) gives
Summing up (27) from to , we have
Taking , the above implies that
where we have used , and denotes the global optimal value. Note that the right hand side above is bounded from above because . Hence, using Lemma 1 we conclude:
Appendix C Proof of Proposition 2
For the constrained case, let us denote as the objective function. Unlike the unconstrained case where we measure convergence via observing if the gradient vanishes, the optimality condition of the constrained case is a bit more complicated. Consider the following optimization problem
where is continuously differentiable while is convex but possibly nonsmooth. The deterministic proximal gradient algorithm for handling this problem is as follows:
Define , the update can also be represented as which is analogous to the gradient descent algorithm. It can be shown that implies that the necessary optimality condition is satisfied, and thus can be considered as a “generalized gradient”. In the multi-block setting of (13), we define:
To show that the BrasCPD algorithm finds a stationary point, our goal is to show the subsequence convergence of to zero for all as .
Our update is equivalent to the following:
for a randomly selected , i.e., the above is the proximal operator. For a given , we have
by the optimality of for solving Problem (29).
By the block Lipschitz continuity of the smooth part (cf. Fact 2), we have
where denotes the smooth part in the objective function and
Combining the two inequalities, we have
where we have defined:
The inequality in (C-B) can be further written as
Again, taking expectation conditioning on the filtration and , we can upper bound by
where (a) is due to the Cauchy-Schwartz inequality, and (b) is a consequence of the non-expansiveness of the proximal operator of convex . Taking the total expectation, we have
Summing up the inequality from to ,
Since , we have , therefore,
such that . Taking , and by the assumption that , we can conclude that
using Lemma 1.
To complete the proof, we observe that
where the last inequality is obtained via applying the nonexpansive property again. Note that both terms on the right hand side converge to zero. Hence, this relationship implies that
Note that by our sampling strategy, we have
However, since is independent of the random seed , we have
This proves the proposition.
Appendix D Proof of Proposition 3
The insight of the proof largely follows the technique for single-block Adagrad , with some careful modifications to multiple block updates. One will see that the block sampling strategy and the block-wise unbiased gradient estimation are key to apply the proof techniques developed in  to our case. To show convergence, let us first consider:
 Let , , and . Then, we have
The proof is simple and elegant; see [55, Lemma 4].
 Consider a random variable . If , then
Let us consider the block-wise again:
Plugging in our update rule under AdaCPD, one can see that
Taking expectation w.r.t. (the random seed that is responsible for selecting fibers) conditioning on the filtration and the selected block , the middle term is zero—since the block stochastic gradient is unbiased [cf. Fact 1]. Hence, we have reached the following
Taking total expectation on both sides, we have
From the above inequality and the assumption that is bounded from above by , we can conclude that
by summing up all the inequalities in (40) from to .
Taking and observe that:
Note that we have exchanged the order of the limits and expectations, since the expectation is taking on nonnegative terms. Using Lemma 2, one can easily show the first term above is bounded from above by , where is a constant. To see the second term is bounded, observe
where means taking expectation w.r.t. all the random variables except for for , and the second inequality is due to the effect of the telescope summation. Since we have assumed that ’s are bounded, the right hand side is bounded from above. Therefore, we have reached the conclusion
Applying Lemma 3, one can see that
Since , one immediate result is that any appears infinitely many times in the sequence , according to the second Borel-Cantelli lemma. This leads to
holds for , where is the subsequence of such that block is sampled for updating. Hence, with probability one there exists a subsequence such that at the corresponding iterations block is sampled for updating. It is not hard to show that
by the assumption that are all bounded. This directly implies that . Together with Lemma 1, we have
-  T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM review, vol. 51, no. 3, pp. 455–500, 2009.
-  N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos, “Tensor decomposition for signal processing and machine learning,” IEEE Trans. Signal Process., vol. 65, no. 13, pp. 3551–3582, 2017.
-  J. D. Carroll and J.-J. Chang, “Analysis of individual differences in multidimensional scaling via an N-way generalization of “Eckart-Young” decomposition,” Psychometrika, vol. 35, no. 3, pp. 283–319, 1970.
-  B. Rao and K. Kreutz-Delgado, “An affine scaling methodology for best basis selection,” IEEE Trans. Signal Process., vol. 47, no. 1, pp. 187 –200, jan 1999.
-  J. Sun, D. Tao, and C. Faloutsos, “Beyond streams and graphs: dynamic tensor analysis,” in Proceedings of the 12th ACM SIGKDD international conference on Knowledge discovery and data mining. ACM, 2006, pp. 374–383.
-  C. I. Kanatsoulis, X. Fu, N. D. Sidiropoulos, and W.-K. Ma, “Hyperspectral super-resolution: A coupled tensor factorization approach,” IEEE Trans. Signal Process. to appear, 2018.
-  A. Anandkumar, R. Ge, D. Hsu, S. M. Kakade, and M. Telgarsky, “Tensor decompositions for learning latent variable models,” The Journal of Machine Learning Research, vol. 15, no. 1, pp. 2773–2832, 2014.
-  A. Anandkumar, D. Hsu, and S. M. Kakade, “A method of moments for mixture models and hidden markov models,” in Conference on Learning Theory, 2012, pp. 33–1.
-  N. D. Sidiropoulos, R. Bro, and G. B. Giannakis, “Parallel factor analysis in sensor array processing,” IEEE Trans. Signal Process., vol. 48, no. 8, pp. 2377–2388, Aug. 2000.
-  N. D. Sidiropoulos and X.-Q. Liu, “Identifiability results for blind beamforming in incoherent multipath with small delay spread,” IEEE Trans. Signal Process., vol. 49, no. 1, pp. 228–236, Jan. 2001.
-  X. Fu, N. D. Sidiropoulos, J. H. Tranter, and W.-K. Ma, “A factor analysis framework for power spectra separation and emitter localization,” IEEE Trans. Signal Process., vol. 63, no. 24, pp. 6581–6594, 2015.
-  C. J. Hillar and L.-H. Lim, “Most tensor problems are np-hard,” Journal of the ACM (JACM), vol. 60, no. 6, p. 45, 2013.
-  K. Huang, N. D. Sidiropoulos, and A. P. Liavas, “A flexible and efficient algorithmic framework for constrained matrix and tensor factorization,” IEEE Trans. Signal Process., vol. 64, no. 19, pp. 5052–5065, 2016.
-  Y. Xu and W. Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” SIAM Journal on imaging sciences, vol. 6, no. 3, pp. 1758–1789, 2013.
-  A. P. Liavas and N. D. Sidiropoulos, “Parallel algorithms for constrained tensor factorization via alternating direction method of multipliers,” IEEE Trans. Signal Process., vol. 63, no. 20, pp. 5450–5463, 2015.
-  C. Navasca, L. De Lathauwer, and S. Kindermann, “Swamp reducing technique for tensor decomposition.” in EUSIPCO, 2008, pp. 1–5.
-  P. Comon, X. Luciani, and A. L. De Almeida, “Tensor decompositions, alternating least squares and other tales,” Journal of Chemometrics: A Journal of the Chemometrics Society, vol. 23, no. 7-8, pp. 393–405, 2009.
-  U. Kang, E. Papalexakis, A. Harpale, and C. Faloutsos, “Gigatensor: scaling tensor analysis up by 100 times-algorithms and discoveries,” in Proc. ACM SIGKDD 2012, 2012, pp. 316–324.
-  E. E. Papalexakis, C. Faloutsos, and N. D. Sidiropoulos, “Parcube: Sparse parallelizable tensor decompositions,” in Joint European Conference on Machine Learning and Knowledge Discovery in Databases. Springer, 2012, pp. 521–536.
-  D. Mitchell, N. Ye, and H. De Sterck, “Nesterov acceleration of alternating least squares for canonical tensor decomposition,” arXiv preprint arXiv:1810.05846, 2018.
-  A. L. Alexander, J. E. Lee, M. Lazar, and A. S. Field, “Diffusion tensor imaging of the brain,” Neurotherapeutics, vol. 4, no. 3, pp. 316–329, 2007.
-  A. Shashua and T. Hazan, “Non-negative tensor factorization with applications to statistics and computer vision,” in Proceedings of the 22nd international conference on Machine learning. ACM, 2005, pp. 792–799.
-  L. Bottou, F. E. Curtis, and J. Nocedal, “Optimization methods for large-scale machine learning,” SIAM Review, vol. 60, no. 2, pp. 223–311, 2018.
-  N. Vervliet and L. De Lathauwer, “A randomized block sampling approach to canonical polyadic decomposition of large-scale tensors,” IEEE J. Sel. Topics Signal Process., vol. 10, no. 2, pp. 284–295, 2016.
-  A. Beutel, P. P. Talukdar, A. Kumar, C. Faloutsos, E. E. Papalexakis, and E. P. Xing, “Flexifact: Scalable flexible factorization of coupled tensors on hadoop,” in Proc. SIAM SDM 2014. SIAM, 2014, pp. 109–117.
-  C. Battaglino, G. Ballard, and T. G. Kolda, “A practical randomized CP tensor decomposition,” SIAM Journal on Matrix Analysis and Applications, vol. 39, no. 2, pp. 876–901, 2018.
-  A. Beck and L. Tetruashvili, “On the convergence of block coordinate descent type methods,” SIAM journal on Optimization, vol. 23, no. 4, pp. 2037–2060, 2013.
-  Y. Nesterov, “Efficiency of coordinate descent methods on huge-scale optimization problems,” SIAM Journal on Optimization, vol. 22, no. 2, pp. 341–362, 2012.
-  S. Ghadimi and G. Lan, “Accelerated gradient methods for nonconvex nonlinear and stochastic programming,” Mathematical Programming, vol. 156, no. 1-2, pp. 59–99, 2016.
-  ——, “Stochastic first-and zeroth-order methods for nonconvex stochastic programming,” SIAM Journal on Optimization, vol. 23, no. 4, pp. 2341–2368, 2013.
-  M. Razaviyayn, M. Hong, and Z.-Q. Luo, “A unified convergence analysis of block successive minimization methods for nonsmooth optimization,” SIAM Journal on Optimization, vol. 23, no. 2, pp. 1126–1153, 2013.
-  H. Wang and A. Banerjee, “Randomized block coordinate descent for online and stochastic optimization,” arXiv preprint arXiv:1407.0107, 2014.
-  Y. Xu and W. Yin, “Block stochastic gradient iteration for convex and nonconvex optimization,” SIAM Journal on Optimization, vol. 25, no. 3, pp. 1686–1716, 2015.
-  J. Duchi, E. Hazan, and Y. Singer, “Adaptive subgradient methods for online learning and stochastic optimization,” Journal of Machine Learning Research, vol. 12, no. Jul, pp. 2121–2159, 2011.
-  X. Fu, C. Gao, H.-T. Wai, and K. Huang, “Block-randomized stochastic proximal gradient for constrained low-rank tensor factorization,” in submitted to IEEE ICASSP 2019, 2019.
-  E. C. Chi and T. G. Kolda, “On tensors, sparsity, and nonnegative factorizations,” SIAM Journal on Matrix Analysis and Applications, vol. 33, no. 4, pp. 1272–1299, 2012.
-  S. A. Vorobyov, Y. Rong, N. D. Sidiropoulos, and A. B. Gershman, “Robust iterative fitting of multilinear models,” IEEE Trans. Signal Process., vol. 53, no. 8, pp. 2678–2689, Aug 2005.
-  X. Fu, K. Huang, W.-K. Ma, N. Sidiropoulos, and R. Bro, “Joint tensor factorization and outlying slab suppression with applications,” IEEE Trans. Signal Process., vol. 63, no. 23, pp. 6315–6328, 2015.
-  E. E. Papalexakis, U. Kang, C. Faloutsos, N. D. Sidiropoulos, and A. Harpale, “Large Scale Tensor Decompositions: Algorithmic Developments and Applications,” IEEE Data Engineering Bulletin, Special Issue on Social Media and Data Analysis, vol. 36, no. 3, pp. 59–66, Sep. 2013.
-  N. Ravindran, N. D. Sidiropoulos, S. Smith, and G. Karypis, “Memory-efficient parallel computation of tensor and matrix products for big tensor decomposition,” in 2014 48th Asilomar Conference on Signals, Systems and Computers, Nov 2014, pp. 581–585.
-  J. Li, C. Battaglino, I. Perros, J. Sun, and R. Vuduc, “An input-adaptive and in-place approach to dense tensor-times-matrix multiply,” in SC’15: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2015, pp. 1–12.
-  N. Kargas, N. D. Sidiropoulos, and X. Fu, “Tensors, learning, and’kolmogorov extension’for finite-alphabet random vectors,” arXiv preprint arXiv:1712.00205, 2017.
-  S. Ibrahim and X. Fu, “Stochastic optimization for coupled tensor decomposition with applications in statistical learning,” in Proc. IEEE DSW 2019, 2019.
-  N. Kargas and N. D. Sidiropoulos, “Learning mixtures of smooth product distributions: Identifiability and algorithm,” arXiv preprint arXiv:1904.01156, 2019.
-  P. A. Traganitis, A. Pagès-Zamora, and G. B. Giannakis, “Blind multiclass ensemble classification,” IEEE Trans. Signal Process., vol. 66, no. 18, pp. 4737–4752, Sep. 2018.
-  L. Xiao and T. Zhang, “A proximal stochastic gradient method with progressive variance reduction,” SIAM Journal on Optimization, vol. 24, no. 4, pp. 2057–2075, 2014.
-  N. Parikh and S. Boyd, “Proximal algorithms,” Foundations and Trends in optimization, vol. 1, no. 3, pp. 123–231, 2013.
-  J. Duchi, S. Shalev-Shwartz, Y. Singer, and T. Chandra, “Efficient projections onto the -ball for learning in high dimensions,” in Proc. the 25th international conference on Machine learning. ACM, 2008, pp. 272–279.
-  J. B. Kruskal, “Nonmetric multidimensional scaling: a numerical method,” Psychometrika, vol. 29, no. 2, pp. 115–129, 1964.
-  R. Bro and N. D. Sidiropoulos, “Least squares algorithms under unimodality and non-negativity constraints,” Journal of Chemometrics: A Journal of the Chemometrics Society, vol. 12, no. 4, pp. 223–247, 1998.
-  H. Robbins and S. Monro, “A stochastic approximation method,” in Herbert Robbins Selected Papers. Springer, 1985, pp. 102–109.
-  D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
-  M. D. Zeiler, “Adadelta: an adaptive learning rate method,” arXiv preprint arXiv:1212.5701, 2012.
-  T. Dozat, “Incorporating nesterov momentum into adam,” 2016.
-  X. Li and F. Orabona, “On the convergence of stochastic gradient descent with adaptive stepsizes,” arXiv preprint arXiv:1805.08114, 2018.
-  X. Chen, S. Liu, R. Sun, and M. Hong, “On the convergence of a class of adam-type algorithms for non-convex optimization,” arXiv preprint arXiv:1808.02941, 2018.
-  K. Huang, N. Sidiropoulos, E. Papalexakis, C. Faloutsos, P. Talukdar, and T. Mitchell, “Principled neuro-functional connectivity discovery,” in Proc. SIAM SDM 2015, 2015.
-  X. Fu, W.-K. Ma, K. Huang, and N. D. Sidiropoulos, “Blind separation of quasi-stationary sources: Exploiting convex geometry in covariance domain,” IEEE Trans. Signal Process., vol. 63, no. 9, pp. 2306–2320, May 2015.
-  L. D. Lathauwer and J. Castaing, “Blind identification of underdetermined mixtures by simultaneous matrix diagonalization,” IEEE Trans. Signal Process., vol. 56, no. 3, pp. 1096 –1105, Mar. 2008.
-  J. Mairal, F. Bach, J. Ponce, and G. Sapiro, “Online learning for matrix factorization and sparse coding,” The Journal of Machine Learning Research, vol. 11, pp. 19–60, 2010.
Xiao Fu (S’12-M’15) is an Assistant Professor in the School of Electrical Engineering and Computer Science, Oregon State University, Corvallis, Oregon, United States. He received his Ph.D. degree in Electronic Engineering from The Chinese University of Hong Kong (CUHK), Hong Kong, in 2014. He was a Postdoctoral Associate in the Department of Electrical and Computer Engineering, University of Minnesota, Minneapolis, MN, United States, from 2014 to 2017. His research interests include the broad area of signal processing and machine learning. He received a Best Student Paper Award at ICASSP 2014. Two of his co-authored papers received Best Student Paper Awards at IEEE CAMSAP 2015 and IEEE MLSP 2019, respectively. He serves as the treasurer of IEEE Signal Processing Society Oregon Chapter. He also serves as a member of the EURASIP Technical Area Committee in Signal Processing for Multisensor Systems for the term of 2020-2023. He is a member of IEEE.
Shahana Ibrahim received her B.Tech. degree in Electronics and Communication Engineering from National Institute of Technology, Calicut, India, in 2012. She had been working as System Validation Engineer at Texas Instruments, Bengaluru, India, from 2012 to 2017. She received her M.S. degree in Electrical Engineering from Oregon State University, Corvallis, Oregon, United States, in 2019. She is currently pursuing her PhD degree in Electrical Engineering at Oregon State University, Corvallis, Oregon, United States. Her research interests are in the broad areas of statistical machine learning and signal processing.
Hoi-To Wai (S’11–M’18) received his PhD degree from Arizona State University (ASU) in Electrical Engineering in Fall 2017, B. Eng. (with First Class Honor) and M. Phil. degrees in Electronic Engineering from The Chinese University of Hong Kong (CUHK) in 2010 and 2012, respectively. He is an Assistant Professor in the Department of Systems Engineering & Engineering Management at CUHK. He has held research positions at ASU, UC Davis, Telecom ParisTech, Ecole Polytechnique, LIDS, MIT. Hoi-To’s research interests are in the broad area of signal processing, machine learning and distributed optimization, with a focus of their applications to network science. His dissertation has received the 2017’s Dean’s Dissertation Award from the Ira A. Fulton Schools of Engineering of ASU and he is a recipient of a Best Student Paper Award at ICASSP 2018.
Cheng Gao is a Computer Science PhD student at University of Missouri, Columbia, Missouri, United States. He is currently focusing his research on sparse learning and bioinformatics. He earned his Master’s degree in Electrical and Computer Engineering from Oregon State University, Corvallis, Oregon, United States, in 2019 with research focus on large-scale tensor decomposition algorithms. He received his Bachelor’s degree in Automation from Wuhan University of Technology, China in 2016.
Kejun Huang received his B.Eng. degree in communication engineering from the Nanjing University of Information Science and Technology, China, in 2010 and his Ph.D. degree in electrical engineering from the University of Minnesota, Minneapolis, in 2016. He is an assistant professor in the Department of Computer and Information Science and Engineering at the University of Florida, Gainesville. He was a postdoctoral associate in the Department of Electrical and Computer Engineering at the University of Minnesota, Minneapolis, from 2016 to 2018. His research interests include machine learning, signal processing, optimization, and statistics. He is a Member of the IEEE.