Any Pauli string has exactly one non-zero element in each row and column. Moreover, if you switch between $I$ and $Z$, and between $X$ and $Y$ in a Pauli string then the pattern of zeros will be the same. I'll call it a type of a Pauli string. Also note that non-zero elements patterns never intersect for different Pauli string types.
The position of non-zero element tells you what type a Pauli string has. For example, in $I\otimes X \otimes X \otimes I \otimes I$ the 0-based index of the unique non-zero element in the first column is $(01100)_2$. The Pauli string $Z\otimes Y \otimes X \otimes I \otimes Z$ has the same pattern of zeros.
So, by looking at indices (in binary notation) of non-zero elements in the first column of $M$ you immediately know what types of Pauli strings are in the decomposition of $M$. And you immediately know how matrix $M$ decomposes into a sum $\sum_t M_t$, where $M_t$ is a weighted sum of Pauli strings of the same type $t$.
Thus, the problem reduces to finding a decomposition of $M_t$ into an exact weighted sum of Pauli strings of type $t$.
There are $2^n$ Pauli strings in each type. And matrix $M_t$ has $2^n$ non-zero elements. In essence, we have a system of linear equations here. There is no way to solve it fast in general. But if we know that $k$ is small and coefficients are restricted, then we can try some specific methods, I guess. Though, it still looks like a knapsack problem.
One can try to solve it just for the type $I\otimes I \otimes \dots \otimes I$. In this case $M_t$ is a diagonal matrix and we have to find its decomposition into a weighted sum of $Z^{b_1}\otimes Z^{b_2} \otimes \dots \otimes Z^{b_n}$ where $b_i$ are bits.
Update
Actually, it might be solvable. I'll explain just for the case where $M_t$ is diagonal. Let's take first $m$ bits and group strings
$$
Z^{b_1}\otimes Z^{b_2} \otimes \dots \otimes Z^{b_{m}} \otimes Z^{c_{m+1}} \otimes \dots \otimes Z^{c_{n}}
$$
together, where $b_i$ are fixed and $c_j$ vary (so there are $2^m$ groups).
A linear combination of matrices in a group is some matrix
$$
Z^{b_1}\otimes Z^{b_2} \otimes \dots \otimes Z^{b_{m}} \otimes D_{b_1\dots b_m}
$$
where $D_{b_1\dots b_m}$ can be any diagonal matrix of size $2^{n-m} \times 2^{n-m}$. It's true for each group.
What we want to achieve is to decompose $M_t$ into "groups" at first. That is, to find decomposition
$$
M_t = \sum_{b_1\dots b_m} Z^{b_1}\otimes Z^{b_2} \otimes \dots \otimes Z^{b_{m}} \otimes D_{b_1\dots b_m},
$$
for some $D_{b_1\dots b_m}$. Note that if $k$ is small ($k<2^m$), then some groups will not be present in the decomposition, i.e. $D_{b_1\dots b_m}=0$.
Now, consider the first element on the diagonal in $D_{b_1\dots b_m}$. In $Z^{b_1}\otimes Z^{b_2} \otimes \dots \otimes Z^{b_{m}} \otimes D_{b_1\dots b_m}$ it will appear on positions with indices $x2^{n-m}$, $x=0,\dots,2^{m}-1$, with some $\pm$ sign. Same is true for each combination $b_1\dots b_m$.
Thus, if we consider elements of $M_t$ on those positions $x2^{n-m}$, we get a system of linear equations for those $2^m$ first elements of $D_{b_1\dots b_m}$. Which we just solve. If in the solution the first element of $D_{b_1\dots b_m}$ is not zero, then some Pauli strings in $M_t$ decomposition must be in this group. If it's zero, then the group is either missing, or has at least two Pauli strings which are parts of $M_t$. This already gives us a lot of information.
To refine and compose it we can do the same for the second, third, and other elements of each $D_{b_1\dots b_m}$. Or try a different split of qubits into two parts. It seems we only need a small number of runs.
