Invertibility is not really a big deal because almost any Gaussian, non-invertible MA$(q)$ model can be changed to an invertible MA$(q)$ model representing the same process by changing the parameter values. This is mentioned in most textbooks for the MA(1) model but it is true more generally.
As an example, consider the MA(2) model
$$
z_t = (1-0.2B)(1-2B)w_t, \tag{1}
$$
where $w_t$ is white noise with variance $\sigma_w^2$. This is not an invertible model because $\theta(B)$ has one root equal to 0.5 inside the unit circle. However, consider the alternative MA(2) model obtained by changing this root to its reciprocal value of 2 such that model takes the form
$$
z_t = (1-0.2B)(1-0.5 B)w_t' \tag{2}
$$
where $w_t'$ has variance $\sigma_w'^2=4\sigma_w^2$. You can easily verify that models (1) and (2) both have the same autocovariance functions and hence specify the same distribution for the data if the process is Gaussian.
To make the model identifiable such that there is a one-to-one mapping from $\theta_1,\theta_2,\dots,\theta_q,\sigma_w^2$ to the distribution of the data, the parameter space is therefore by convention restricted to that of invertible models. This particular convention is preferred since the model then can be put directly in AR$(\infty)$ form with coefficients $\pi_1,\pi_2,\dots$ satisfying the simple difference equation $\theta(B)\pi_i=0$.
If we didn't impose this restriction on the parameter space, the likelihood function of an MA$(q)$ would in general have up to $2^q$ local optima (if the MA polynomial has $q$ distinct real roots) which is something we want to avoid.
You can always move roots from inside to outside the unit circle with a corresponding change in the white noise variance using the above technique, except in cases where the MA-polynomial has one or more roots exactly on the unit circle.
