  We give efficient algorithms for finding power-sum decomposition of an input
polynomial $P(x)= \sum_{i\leq m} p_i(x)^d$ with component $p_i$s. The case of
linear $p_i$s is equivalent to the well-studied tensor decomposition problem
while the quadratic case occurs naturally in studying identifiability of
non-spherical Gaussian mixtures from low-order moments.
  Unlike tensor decomposition, both the unique identifiability and algorithms
for this problem are not well-understood. For the simplest setting of quadratic
$p_i$s and $d=3$, prior work of Ge, Huang and Kakade yields an algorithm only
when $m \leq \tilde{O}(\sqrt{n})$. On the other hand, the more general recent
result of Garg, Kayal and Saha builds an algebraic approach to handle any
$m=n^{O(1)}$ components but only when $d$ is large enough (while yielding no
bounds for $d=3$ or even $d=100$) and only handles an inverse exponential
noise.
  Our results obtain a substantial quantitative improvement on both the prior
works above even in the base case of $d=3$ and quadratic $p_i$s. Specifically,
our algorithm succeeds in decomposing a sum of $m \sim \tilde{O}(n)$ generic
quadratic $p_i$s for $d=3$ and more generally the $d$th power-sum of $m \sim
n^{2d/15}$ generic degree-$K$ polynomials for any $K \geq 2$. Our algorithm
relies only on basic numerical linear algebraic primitives, is exact (i.e.,
obtain arbitrarily tiny error up to numerical precision), and handles an
inverse polynomial noise when the $p_i$s have random Gaussian coefficients.
  Our main tool is a new method for extracting the linear span of $p_i$s by
studying the linear subspace of low-order partial derivatives of the input $P$.
For establishing polynomial stability of our algorithm in average-case, we
prove inverse polynomial bounds on the smallest singular value of certain
correlated random matrices with low-degree polynomial entries that arise in our
analyses.
