Is it possible to decompose a hermitian and unitrary matrix $A$ into the sum of the Pauli matrix Kronecker products?
For example, I have a matrix 16x16 and want it to be decomposed into something like $$A =\sum_{i,j,k,l}h_{ijkl}\sigma_i\otimes\sigma_j\otimes\sigma_k\otimes\sigma_l$$ I can't find any clear explanation how it can be done. I found only this for two-dimentional case: https://michaelgoerz.net/notes/decomposing-two-qubit-hamiltonians-into-pauli-matrices.html
For any matrix $A$ we can write $$ A =\sum_{i,j,k,l}h_{ijkl}\cdot \frac{1}{4}\sigma_i\otimes\sigma_j\otimes\sigma_k\otimes\sigma_l, $$ where $$ h_{ijkl} = \frac{1}{4}\text{Tr}\big((\sigma_i\otimes\sigma_j\otimes\sigma_k\otimes\sigma_l)^\dagger \cdot A\big) = \frac{1}{4}\text{Tr}\big((\sigma_i\otimes\sigma_j\otimes\sigma_k\otimes\sigma_l) \cdot A\big) $$ This is because 256 matrices $\frac{1}{4}\sigma_i\otimes\sigma_j\otimes\sigma_k\otimes\sigma_l$ form an orthonormal basis in the linear space of all $16 \times 16$ matrices under Hilbert-Schmidt inner product given by $( A,B) = \text{Tr}(B^\dagger A)$.
If you have an orthonormal basis $e_i$ then any element of a Hilbert space can be written as $v = \sum_i (v,e_i)e_i$. That's it.
If dimension is not power of 2, then it's impossible, of course. But there are generalizations of Pauli matrices https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices. So the similar formula can be written.
$\newcommand{\bs}[1]{{\boldsymbol #1}} \newcommand{\tildebssigma}{\tilde{\bs\sigma}} \newcommand{\bssigma}{{\bs\sigma}}$Yes, products of Pauli matrices form a basis for the set of Hermitian matrices (of dimensions that are powers of $2$).
More specifically, fix an integer $n$ and let $N\equiv 2^n$, define $\bssigma\equiv(\sigma_x,\sigma_y,\sigma_y)$, and $\tildebssigma\equiv (I,\sigma_x,\sigma_y,\sigma_y)$. Moreover, let $J\in\{0,1,2,3\}^{n}$ be a tuple of $n$ integers with each $J_i\in\{0,1,2,3\}$. Consider the matrices of the form $$ \tildebssigma_J\equiv\prod_{k=1}^n \tildebssigma^{(k)}_{J_k}, $$ where $\tildebssigma_j^{(k)}$ denotes the Pauli matrix $\tildebssigma_j$ applied on the $k$-th qubit. You can check that these are all Hermitians. Moreover, $\tildebssigma_J$ is traceless for all $J\neq(0,...,0)$ (for which $\tildebssigma_{(0,...,0)}=I$).
There are $4^n=N^2$ matrices of this form (one for each possible choice of $J$), and for any pair of tuples $J,K$ we have $\operatorname{Tr}(\tildebssigma_J\tildebssigma_K)=N\delta_{JK}$. Moreover, the space of Hermitian $N\times N$ matrices also has dimension $N^2$. It follows that $\{\tildebssigma_J\}_J$ is a basis for this space. Explicitly, you can decompose an arbitrary Hermitian matrix $H$ as $$H=\frac{1}{N}\sum_J \operatorname{Tr}(\tildebssigma_J H)\tildebssigma_J.$$ Note that the coefficients in any such expansion are always real. This is not by chance. Indeed, being more careful, we should state that the set of real linear combinations of products of Pauli matrices give the set of Hermitian matrices (notice that if $A$ is Hermitian then $\alpha A$ is Hermitian iff $\alpha\in\mathbb R$, so this is not surprising).
More general matrices can be generated if we allow for complex coefficients in the expansion. Indeed, note that $$\frac{I+Z}{2}=\begin{pmatrix}1&0 \\ 0 & 0\end{pmatrix}, \qquad \frac{I-Z}{2}=\begin{pmatrix}0&0 \\ 0 & 1\end{pmatrix}, \\ \frac{X+iY}{2}=\begin{pmatrix}0 & 1 \\ 0 & 0\end{pmatrix} \qquad \frac{X-iY}{2}=\begin{pmatrix}0 & 0 \\ 1 & 0\end{pmatrix}. $$ Therefore the complex span of Pauli matrices can be used to generate arbitrary $2\times 2$ matrices. This then translates into the same result for arbitrary $2^n$-dimensional spaces, as if $\mathcal V\equiv\{v_k\}$ is a basis for $V$, then the sets of tensor products of elements of $\mathcal V$ form a basis for $V^{\otimes n}$.
This is a special instance of the more general fact that any matrix can be decomposed uniquely as sum of a Hermitian and a skew-Hermitian matrix, as $$A=\frac{A+A^\dagger}{2}+i\frac{A-A^\dagger}{2i},$$ and the fact that products of Pauli matrices give you a basis for the set of Hermitian matrices.
