Suppose that $A$ is a symmetric non-random matrix and $X\sim N(\mu,\Sigma)$ and $b \in R^n$ is a non-random vector. Then what is the distribution of $$X^tAX+b^tX \quad ?$$
The distribution without the linear term is solved in the answer here(Transformation of multivariate normal sum of chi-squared).
In the case of an invertible $A$ we can write $X^tAX+b^tX=(X-h)^tA(X-h)+k$ where $h=-\frac{1}{2}A^{-1}b$ and $k=-\frac{1}{4}b^tA^{-1}b$. However, the case where $A$ is not invertible is also of interest as it arises in practice. In the case where $n=1$ this corresponds to $A=0$ and thus the distribution above is simply a normal with mean $\mu_0 = b^t\mu$ and $\sigma_0 = b^t\Sigma b$. Is there a reducible solution in the more general case of $n>1$ for arbitrary non-invertible $A$? Perhaps an appropriate transformation can disentangle the quadratic form and the linear so that we have in some basis independent sum a normal and a linear combination of scaled non-central chi-squared with 1 degree of freedom.
Attempt: Write $A=P\Lambda P^t$, where $\Lambda$ is a diagonal of the eigenvalues and $P$ has the corresponding unit eigenvectors in it's rows, $PP^t = I$. As $A$ is not invertible there are $r>0$ zero eigenvalues. Then we can write $$X^tAX+b^tX = X^tP\Lambda P^tX+b^tPP^tX $$ Set $Y=P^tX$, then $$ X^tAX+b^tX = Y^t \Lambda Y + b^t PY$$ Assume w.l.o.g. that the eigenvalues which are $0$ and corresponding eigenvectors in $P$ are the last $r$ rows/columns in $\Lambda$. Then $$X^tAX+b^tX = (Y^\star)^t \Lambda^\star Y^\star + (b^\star)^tY^\star+(b')^tY' $$ $Y^\star$ are the $n-r$ first entries of $Y$ and $Y'$ the rest, $b^\star$ is the $n-r$ first entries of $P^tb$ and $b'$ the $r$ last. $\Lambda^\star \in R^{(n-r)\times (n-r)}$. Now the first two terms can be used to fill the square and the last term involves parts of $Y$ not involved in the first two terms but it is not independent of them as the covariance matrix of $Y$, $Cov(Y) = P^t\Sigma P$ is not diagonal. The goal is to identify the distribution of $X^tAX+b^tX $.
