  The canonical polyadic decomposition (CPD) is a compact decomposition which
expresses a tensor as a sum of its rank-1 components. A common step in the
computation of a CPD is computing a generalized eigenvalue decomposition (GEVD)
of the tensor. A GEVD provides an algebraic approximation of the CPD which can
then be used as an initialization in optimization routines.
  While in the noiseless setting GEVD exactly recovers the CPD, it has recently
been shown that pencil-based computations such as GEVD are not stable. In this
article we present an algebraic method for approximation of a CPD which greatly
improves on the accuracy of GEVD. Our method is still fundamentally
pencil-based; however, rather than using a single pencil and computing all of
its generalized eigenvectors, we use many different pencils and in each pencil
compute generalized eigenspaces corresponding to sufficiently well-separated
generalized eigenvalues. The resulting "generalized eigenspace decomposition"
is significantly more robust to noise than the classical GEVD.
  Accuracy of the generalized eigenspace decomposition is examined both
empirically and theoretically. In particular, we provide a deterministic
perturbation theoretic bound which is predictive of error in the computed
factorization.
