$
\def\l{\left}
\def\r{\right}
\def\o{{\tt1}}
\def\p{{\partial}}
\def\g#1#2{\frac{\p #1}{\p #2}}
\def\c#1{\color{red}{#1}}
\def\B{{\mathbb B}}
\def\R{{\mathbb R}^p}
\def\E{{\cal E}}
$Let $\{e_k\}$ denote the standard basis vectors for $\R$ and $M$ the gradient of $q,\;$ then
$$M = \g{q}{x} \quad\implies\quad dq = M\,dx$$
The gradient of a matrix $(G)$ with respect to a vector $(x)$ is a third-order tensor $(\B)$ with components
$$\eqalign{
\B_{ijk} &= \g{G_{ij}}{x_k} \\
}$$
Using this tensor we can write the other terms in the problem as
$$\eqalign{
B_k &= \B\cdot e_k \;\;\implies\;\; \B = \sum_{k=1}^p B_k\star e_k \\
\B\cdot x &= \sum_{k=1}^p B_k\star\l(e_k\cdot x\r) = \sum_{k=1}^p B_k x_k \\
G &= A + \B\cdot x \\
}$$
where $(\star)$ denotes the tensor/dyadic product,
and $(\cdot)$ the dot standard product.
Let's also define a new vector $(w)$ with components
$$\eqalign{
w_k = \sum_{k=1}^p G^{-1}qq^TG^{-1}:B_k \\
}$$
where $(:)$ denotes the double-dot product and is an alternative
notation for the trace, i.e.
$$\eqalign{
A:B &= \sum_{i=1}^m\sum_{j=1}^n A_{ij}B_{ij}
\;=\; {\rm Tr}\!\l(AB^T\r) \\
A:A &= \big\|A\big\|_F^2 \\
}$$
Now we can write the function and calculate its differential and gradient.
$$\eqalign{
f &= q^TG^{-1}q \\
&= G^{-1}:qq^T \\
df&= G^{-1}:\l(dq\,q^T+q\,dq^T\r) + qq^T:dG^{-1} \\
&= \l(G^{-1}+G^{-T}\r):dq\,q^T - qq^T:G^{-1}dG\,G^{-1} \\
&= 2G^{-1}:dq\,q^T - G^{-1}qq^TG^{-1}:dG \\
&= 2G^{-1}q:M\,dx - G^{-1}qq^TG^{-1}:\B\cdot dx \\
&= \l(2MG^{-1}q - G^{-1}qq^TG^{-1}:\B\r)\cdot dx \\
&= \l(2MG^{-1}q - w\r)\cdot dx \\
\g{f}{x} &= 2MG^{-1}q - w \\
}$$