I have seen in some ML tutorial that functions of the form $f(\vec x) = \vec x^T A \vec x$ ($\vec x \in \mathbb{R}^n$ and $A$ is an $n\times n$ real matrix) can be PAC learned.

Can anyone point me to a reference to this fact?

One approach is to use linear regression.

The matrix $A$ has $n^2$ entries; treat each of those entries as an unknown, so you have $n^2$ unknowns. Now when $x$ is known, $f(x)$ is a linear function of those unknowns. If you have a set of equations $f(x_i)=y_i$ that are all known to be exactly correct, then you can solve for the unknowns by using linear algebra. If they are known to be approximately correct (e.g., $||f(x_i)-y_i||_2$ is small), you can use least squares to find the solution for the unknowns that minimizes the sum of the squared errors (that minimizes $\sum_i ||f(x_i)-y_i||_2$). If a $1-\epsilon$ fraction of them are approximately correct, where $\epsilon \ll 1/n$, then you can pick a random subset of $n$ of them, guess that all of the ones in the subset are correct, use least squares to solve for the unknowns, and check whether that solution seems valid given the remaining equations, and repeat until you find a good solution.

Perhaps to complement D.W.'s answer: given an unknown matrix $A$, we are interested in finding a consistent algorithm. That is, given a list of $m$ vectors and their values $(\vec x_1,v_1),\dots,(\vec x_m,v_m)$, find some matrix $A^*$ such that for every $\vec x_j$ we have $\vec x_j^TA^* \vec x_j = v_j$ for all $j \in [m]$. Since $\vec x_j$'s are given, we just need to solve the following linear system: find $(a_{qr})_{q,r \in [n]}$ such that
$$\forall j \in [m]: \sum_{q = 1}^n\sum_{r = 1}^na_{qr}x_{j}(q)x_{j}(r) = v_j$$ (here $x_j(q)$ is the $q$-th coordinate of $\vec x_j$). Since we know that

1. There exists some solution (namely $A$ that we do not know)
2. This is a linear system with $\mathcal O(n^2)$ variables.

We can find a solution $A^*$ to the above system in poly time. In particular, this implies that there is a consistent algorithm. Now, suppose that we limit ourselves to matrices that take integer values in $[-M,M]$; then there are $(2M)^{n^2}$ possible hypotheses. In other words, if our hypothesis class $\mathcal C$ is the set of all quadratic functions in the range $[-M,M]$, the log size of $\mathcal C$ is $\mathcal O(n^2\log M)$.

The existence of a consistent algorithm, along with the above bound on the log size of $\mathcal C$ implies that quadratic functions are $(\varepsilon, \delta)$-PAC learnable in time polynomial in $\frac 1\varepsilon,\log \frac1\delta, n^2,$ and $\log M$.

However, I am still looking for a reference to this (seems like a simple enough fact to be cited somewhere).