My guess is, in the hypotesis that $b_k$ are large and hence the "non-linear" term $ \sum_{p=1}^{M}\frac{a_k\,a_q\,a^2_p}{b_k\,b_q}\,x_{k,p}\,x_{q,p}$ is neglectable, to use an iterative alogrithm. I start with a guess where the non-linear term is put to zero:
$$
x^{\left(0\right)}_{k,q} = \frac{1}{\left(\frac{a_k\,a_{q}^2}{b_k}+\frac{a_q\,a_{k}^2}{b_q}\right)}\,\left(c_{k,q}-a^2_k\,\delta_{k,q}\right)
$$