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).
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.)