Symmetric Rank-1 Decomposition for Density Matrices Let $(H,\langle\cdot,\cdot\rangle)$ be an $n$-dimensional complex Hilbert space. For concreteness, you can just take $H=\mathbb{C}^n$ with standard inner product. Note that we will use the physicist's convention of the inner product being complex-linear in the second entry.
Consider the symmetric $k$-fold tensor product product of $H$ with itself, which we will denote by $H^{\otimes_s^k}$. Similarly, consider the symmetric $k$-fold tensor product of the dual of $H$ with itself, which we will denote by $H^{*,\otimes_s^k}$. We are interested in elements of the tensor product $H^{\otimes _s^k} \otimes H^{*,\otimes_s^k}$. Of course using the Riesz representation theorem, we can identify $H^*$ with $H$ using the inner product and so $H^{\otimes _s^k} \otimes H^{*,\otimes_s^k}$ is the complex span of elements of the form
$$|f\rangle\langle g|, \qquad f,g\in H^{\otimes_s^k}, \tag{1}$$
where we have used Dirac's bra-ket notation.
We are not interested in the whole space $H^{\otimes _s^k} \otimes H^{*,\otimes_s^k}$, but only those elements which are self-adjoint, by for an element of the form (1) means $f=g$. My question is the following:

Question. Is it possible to write every self-adjoint element of $H^{\otimes_s^k} \otimes H^{*,\otimes_s^k}$ into a linear combination
  $$\sum_{j=1}^N a_j |f_j^{\otimes k}\rangle\langle f_j^{\otimes_k}|, \tag{2}$$
  where $N\in\mathbb{N}$, $a_j\in\mathbb{C}$, and $f_j\in H$?

I know that it is possible to write every element $f\in H^{\otimes_s^k}$ as
$$f=\sum_{j=1}^N a_j f_j^{\otimes k},$$
and consequently every self-adjoint of element of $H^{\otimes_s^k}\otimes H^{*,\otimes_s^k}$ can be written as
$$\sum_{j=1}^N a_j(|f_j^{\otimes k}\rangle\langle g_j^{\otimes k}| + |g_j^{\otimes k}\rangle\langle f_j^{\otimes k}|).$$
However, I do not know how to prove such a symmetric rank-1 type decomposition like (2).
 A: First, a little discussion on these self-adjoint elements.  While the answer to
your question (about the form (2)) is certainly Yes, it makes more sense to
restrict the coefficients $a_j$ to real numbers rather than complex; the
answer is then still Yes, but with real coefficients the individual terms in
the sum are themselves constrained to be self-adjoint.  The reason is that the
adjoint operation conjugates the coefficients, so that it is compatible with
the complex inner product we are using.
(In linear-algebra language, if that is familiar to you, we can identify the
kets as column vectors, the bras as row vectors, the bra-ket correspondence
and the adjoint operation as the conjugate transpose, and the standard bra-ket
notation as matrix multiplication (so that $\langle f|g\rangle$ is the scalar
product of row vector $\langle f|$ and column vector $|g\rangle$, while
$|f\rangle\langle g|$ is the outer product of column vector and row vector).
So these self-adjoint elements are just Hermitian matrices, equal to their
conjugate transposes.)
Second, as far as I can tell your question is unrelated to the tensor-product
structure of the Hilbert space: the property of self-adjointness is unrelated
to this structure, and so this only complicates the notation by adding
$\otimes k$ superscripts on everything.  I will drop these superscripts.
On to your question: I think you want to know whether paired ``off-diagonal''
terms $|f\rangle\langle g| + |g\rangle\langle f|$ can be written as a linear
combination of terms like $|f_j\rangle\langle f_j|$.  This is straightforward;
one way is to define $|h\rangle=|f\rangle+|g\rangle$; then
$$|h\rangle\langle h|=|f\rangle\langle f|+|g\rangle\langle g|+|f\rangle\langle g|+|g\rangle\langle f|$$
and so
$$|h\rangle\langle h|-|f\rangle\langle f|-|g\rangle\langle g|=|f\rangle\langle g|+|g\rangle\langle f|$$
which is of the form you want.
But!  This is not usually the most useful way to do such a decomposition.  The
problem is that the coefficients, while real, are not all positive.  The
density matrix $\rho$ is positive semidefinite, which means that we really
should strive for a decomposition like (2) with coefficients $a_j$ not just
real, but nonnegative.  The way to do this is with an eigendecomposition of
$\rho$:  Since $\rho$ is Hermitian, it has a complete set of orthonormal
eigenvectors, i.e. $\rho|e_j\rangle=\lambda_j|e_j\rangle$ with $\lambda_j\ge0$
and $\langle e_i|e_j\rangle=\delta_{ij}$.  The eigendecomposition of $\rho$ is
then
$$\rho=\sum_j \lambda_j|e_j\rangle\langle e_j|$$
of the form (2) with all coefficients nonnegative.  (This can be interpreted as
a probabilistic mixture in which the system is in state $|e_j\rangle$ with
probability $\lambda_j$.  But note that the eigendecomposition is not
unique if there are degenerate eigenvalues.)
