Consider a $2^n\times 2^n$ Hermitian matrix $M$ containing up to two non-zero elements, which are $1$ (so, either $M_{ii}=1$ for some $i$, or $M_{ij}=M_{ji} = 1$ for some $i$ and $j$). Each such matrix can be expressed a linear combination of Pauli operators. For example: $$ \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} = \frac{1}{2} (XX+YY) \ ,\quad \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} = \frac{1}{4} (1 - Z_1 + Z_2 -Z_1 Z_2) \ . $$
What would be the general expression for such an expansion, for an arbitrary $n$ and $i,j=1\ldots2^n$?
CLARIFICATION
I'm looking for a solution different from the one which involves taking trace with all the possible Pauli operators. The problem with that solution is that it requires $O(4^n)$ operations, since this is the total number of Paulis on $n$ qubits.
The motivation for my question, of course, comes from expanding an arbitrary Hermitian $H$ in terms of Pauli operators. The number of parameters in such a matrix is $O(2^n)$, which sets the lower bound on the solution complexity — which, I believe, is possible to achieve. One way to do this would be to decompose $M$ from my question in $O(\operatorname{poly}(n)) $ steps. This would allow to expand $H$ in $O(2^n)$ instead of $O(4^n)$ steps, via expressing it in terms of Paulis entry by entry.