In order to make the notation shorter for the remaining calculations in my answer I will denote the matrix with
$$
H = \frac{1}{2}
\begin{bmatrix}
-q^\top \\ q_0\,I_3 + q^\times
\end{bmatrix}.
$$
In this case you have four equations but only three unknowns. So there might not be a solution that exactly matches all equations, especially when the time derivative of the quaternion is perturbed, i.e. by numerical rounding errors. However, there is an elegant least squares solution that minimizes $\|\dot{\textbf{q}} - H\,\omega\|$. Here I use $\textbf{q}$ to denote the entire quaternion $\begin{bmatrix}q_0 & q^\top\end{bmatrix}^\top$. This least squares solution can shown to be
$$
\omega = \left(H^\top H\right)^{-1} H^\top \dot{\textbf{q}}.
$$
It can be noted that if the quaternion is of unit length then the matrix, whose inverse is taken, simplifies to
$$
\left(H^\top H\right)^{-1} = 4\,I_3.
$$
Thus the solution for the angular velocity in a reduced form can also be written as
$$
\omega = 4\,H^\top \dot{\textbf{q}} = 2
\begin{bmatrix}
-q & q_0\,I_3 \!-\! q^\times
\end{bmatrix}
\begin{bmatrix} \dot{q}_0 \\ \dot{q}
\end{bmatrix}.
$$
