Is high-dimensional Pauli decomposition impossible? I've been trying to decompose a $3 \times 3$ density matrix with 3-dimensional Pauli matrices but it doesn't work for all matrices.
For example, the density matrix of the state $|0\rangle + |1\rangle + |2\rangle$ can be decomposed by obtaining the coefficients of the equation $\rho = a_{X}X + a_{Y}Y + a_{Z}Z + a_{V}V + a_{X^2}X^2 + a_{Y^2}Y^2 + a_{Z^2}Z^2 + a_{V^2}V^2 + a_{I}I $ using trace, e.g. $a_X = Tr(\rho X)/3 $. Here, $Y=XZ, V=XZ^2$, and X and Z are $3 \times 3$ Pauli matrices shown in https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices.
On the other hand, when I tried to decompose the state $|0\rangle + |2\rangle$ using the above equation, the result was not the same as $|0\rangle + |2\rangle$.
Is it impossible to decompose a high-dimensional matrix using Pauli matrices?

Cross-posted on qc.SE
 A: For a general $3\times 3$ matrix, which contains 9 elements, you would need 9 standard matrices. For a density matrix, which is normalized, one degree of freedom is removed by the fact that the trace of the matrix is equal to 1. In that case you need 8 matrices. The Gell-mann matrices is a good choice for this purpose. Since there are only 3 Pauli matrices, they cannot represent an arbitrary normalized $3\times 3$ matrix.
A: For a $3$-level system, the proper generalization of the Pauli matrices are the Gell-Mann or the clock-and-shift matrices.  (see also Patera J, Zassenhaus H. The Pauli matrices in n dimensions and finest gradings of simple Lie algebras of type $A_{n-1}$, Journal of Mathematical Physics. 1988 Mar;29(3):665-73.) Explicitly the latter are
\begin{align}
A_0&=\left(
\begin{array}{ccc}
 0 & 1 & 0 \\
 0 & 0 & 1 \\
 1 & 0 & 0 \\
\end{array}\right)\, ,\quad 
A_0^2\, ,\quad B_0=\left(
\begin{array}{ccc}
 0 & \omega  & 0 \\
 0 & 0 & \omega ^2 \\
 1 & 0 & 0 \\
\end{array}
\right)\, ,\quad 
B_0^2\, \\
C_0&=\left(
\begin{array}{ccc}
 0 & \omega ^2 & 0 \\
 0 & 0 & \omega  \\
 1 & 0 & 0 \\
\end{array}
\right)\, ,\quad 
C_0^2\, ,\quad 
D_0=\left(
\begin{array}{ccc}
 1 & 0 & 0 \\
 0 & \omega  & 0 \\
 0 & 0 & \omega ^2 \\
\end{array}
\right)\, ,\quad
D_0^2
\end{align}
and the unit matrix, where $\omega^3=1$.  (The phases are not unique in the definition of each of these matrices.)  As defined they satisfy Tr$[T_i^\dagger T_j]/3=\delta_{ij}$.
If your state is $\vert\psi\rangle=\frac{1}{\sqrt{3}}(1,1,1)^\top$, then
$$
\rho=\frac{1}{3}\left(
\begin{array}{ccc}
 1 & 1 & 1 \\
 1 & 1 & 1 \\
 1 & 1 & 1 \\
\end{array}
\right)=\frac{1}{3} \left( A_0+A_0^2+\hat{\mathbb{1}}\right)\, .
$$
If your state is $\vert\psi\rangle=\frac{1}{\sqrt{2}}(1,1,0)^\top$ then
the expansion coefficients in the order given above are
$$
\frac{1}{6}\left(1,1,\omega ^2,\omega ,\omega ,\omega ^2,\omega ^2+1,\omega +1,2\right)\, .
$$
The $3\times 3$ angular momentum matrices for $L_x,L_y$ and $L_z$ (if this is what you mean by higher-dimensional Pauli matrices) do not form a basis for $3\times 3$ Hermitian matrices.
There is a set of mutually unbiased bases in prime dimension of when the dimension is a power of a prime.  This would include dimension $3$ and $8$ (i.e. $3\times 3$ matrices and $8\times 8$ matrices).  The matrices above (excluding the unit) do decompose into 4 sets of 2 commuting matrices: $\{A_0, A_0^-\}$ etc so they form a full set of mutually unbiased matrices.
In power-of-prime dimension (like $8=2^3$) one needs to use finite field methods: see

Gibbons KS, Hoffman MJ, Wootters WK. Discrete phase space based on finite fields. Physical Review A. 2004 Dec 3;70(6):062101.

For $8=2^3$ you would need a cubic extension of $\mathbb{Z}_2$ to proceed.  As there are multiple choices of irreducible polynomials, it's not so easy to check if different sets are equivalent.
There is no known complete set of mutually unbiased bases in composite dimension.  The case $d=6$ was studied extensively as a prototype by Brierley and Weigert, but no complete set was found. See for instance

Brierley S, Weigert S. Maximal sets of mutually unbiased quantum states in dimension 6. Physical Review A. 2008 Oct 13;78(4):042312.

A: A more general perspective:
If you want to know whether a set of vectors $\vec v_i$ spans the full space, you can compute the Gram matrix
$$
G_{ij} = (\vec v_i,\vec v_j)
$$
and determine its rank: It equals the dimension of the space spanned.
In the specific example of a space of matrices, the scalar product $(\cdot,\cdot)$ is given by $\mathrm{tr}(A^\dagger B)$.
