You are given a set of vectors, and you must find the best choice in this set. It is unlikely that an analytical method will give you the answer. Instead, I suggest you find a way to quickly discard vectors that are bad candidates.

Here is what might be the basis for a practical algorithm.

Typically, covariance matrices can be approximated using a low-rank matrix. That is, $A$ is diagonalizable and most of its eigenvalues are near zero relative to the highest eigenvalues.

For simplicity, let me assume that we can approximate $A$ using a one-dimensional projection. That is, we have that $A \approx \lambda y y^{\top}$. Then we can estimate $(x-u)^T A (x-u)$ as $((x-u) \cdot y)^2$. This last expression can be computed much faster. Naturally, you can extend this analysis with more eigenvalues for more accuracy. Let me write $P(x)$ this estimate of $(x-u)^T A (x-u)$. The important thing is that $P(x)$ can be computed much faster than $(x-u)^T A (x-u)$ and is somewhat accurate. Find $\kappa$ such that $| (x-u)^T A (x-u)- ((x-u) \cdot y)^2 |\leq \kappa$.

This can be used to quickly prune out $x$'s that cannot possibly be the best choice. That is, if $\max_x P(x) = M$ then any $x$ such that $ P(x) < M-\kappa$ can be rejected.

This leaves you with a smaller sets of vectors $x$ over which you can do the full computation: $\arg \max (x-u)^T A (x-u)$.

You are given a set of vectors, and you must find the best choice in this set. It is unlikely that an analytical method will give you the answer. Instead, I suggest you find a way to quickly discard vectors that are bad candidates.

Here is what might be the basis for a practical algorithm.

Typically, covariance matrices can be approximated using a low-rank matrix. That is, $A$ is diagonalizable and most of its eigenvalues are near zero relative to the highest eigenvalues.

For simplicity, let me assume that we can approximate $A$ using a one-dimensional projection. That is, we have that $A \approx \lambda y y^{\top}$. Then we can estimate $(x-u)^T A (x-u)$ as $((x-u) \cdot y)^2$. This last expression can be computed much faster. Naturally, you can extend this analysis with more eigenvalues for more accuracy. Let me write $P(x)$ this estimate of $(x-u)^T A (x-u)$. The important thing is that $P(x)$ can be computed much faster than $(x-u)^T A (x-u)$ and is somewhat accurate. Quickly find $\kappa$ such that $| (x-u)^T A (x-u)- ((x-u) \cdot y)^2 |\leq \kappa$. You might be able to compute $\kappa$ from the eigenvalues of $A$.

This fast estimation be used to quickly prune out $x$'s that cannot possibly be the best choice. That is, if $\max_x P(x) = M$ then any $x$ such that $ P(x) < M-\kappa$ can be rejected.

This leaves you with a smaller sets of vectors $x$ over which you can do the full computation: $\arg \max (x-u)^T A (x-u)$.