I am reading a book on Hilbert spaces and it mentions the following fact about quadratic forms and bilinear functionals:
If $\hat{\psi}$ is a the quadratic form induced by a bilinear functional $\psi$ on a complex vector space then we have that
$$\psi(x, y) = \hat{\psi}\bigg(\frac{1}{2}(x + y)\bigg) - \hat{\psi}\bigg(\frac{1}{2}(x - y)\bigg) + i\hat{\psi}\bigg(\frac{1}{2}(x + iy)\bigg) - i\hat{\psi}\bigg(\frac{1}{2}(x - iy)\bigg).$$
Where does this identity come from?
It is called Polarization Identity. The $1/2$ should be $1/4$. Also, the formula applies to sesquilinear forms rather than bilinear. It is usually written as $$ \psi(x,y)=\frac14\,\sum_{k=0}^3i^k\,\hat\psi(x+i^k y). $$
Not sure what you expect regarding where "it comes from". It is the identity that relates the inner product with the norm in a Hilbert space. From that point of view, it could be written as $$ \langle x,y\rangle=\sum_{k=0}^3i^k\,\|x+i^ky\|^2. $$ It is likely that the first mathematicians who considered complex Hilbert spaces figured it out. It is not a far-fetched generalization of the identity in the real case, $$ \langle x,y\rangle=\frac14\,(\|x+y\|^2-\|x-y\|^2). $$