Zagier's one-sentence proof of a theorem of Fermat Zagier has a very short proof (MR1041893, JSTOR) for the fact that every prime number $p$ of the form $4k+1$ is the sum of two squares. 
The proof defines an involution of the set $S= \lbrace (x,y,z) \in N^3: x^2+4yz=p \rbrace $ which is easily seen to have exactly one fixed point. This shows that the involution that swaps $y$ and $z$ has a fixed point too, implying the theorem.
This simple proof has always been quite mysterious to me. Looking at a precursor of this proof by Heath-Brown did not make it easier to see what, if anything, is going behind the scenes. 
There are similar proofs for the representation of primes using some other quadratic forms, with much more involved involutions.
Now, my question is: is there any way to see where these involutions come from and to what extent they can be used to prove similar statements?
 A: This paper by Christian Elsholtz seems to be exactly what you're looking for.  It motivates the Zagier/Liouville/Heath-Brown proof and uses the method to prove some other similar statements.  Here is a German version, with slightly different content.
Essentially, Elsholtz takes the idea of using a group action and examining orbits as given (and why not -- it's relatively common) and writes down the axioms such a group action would have to fulfill to be useful in a proof of the two-squares theorem.  He then algorithmically determines that there is a unique group action satisfying his axioms -- that is, the one in the Zagier proof.  The important thing is that having written down these (fairly natural) axioms, there's no cleverness required; finding the involution in Zagier's proof boils down to solving a system of equations.
A: Zagier's proof is non-constructive. Part II of this paper

*

*Aimeric Malter, Dierk Schleicher and Don Zagier, New Looks at Old Number Theory, AMM 120, No. 3 (2013), pp. 243–264 doi:10.4169/amer.math.monthly.120.03.243, JSTOR
gives a refined version with effective, although not very efficient, algorithm to find the decomposition of a prime number $p=4k+1$ into two squares. See also:

*

*Bhaskar Bagchi, Fermat's two squares theorem revisited, Journal of Science Education 4 (7) (1999) 59–67, doi:10.1007/BF02839014 (repository version).

Impressively fast algorithm for finding $n$ and $m$ in
$p=4k+1=n^2+m^2$ was given in:

*

*John Brillhart, Note on Representing a Prime as a Sum of Two Squares, Mathematics of Computation 26 No. 120 (1972) 1011–1013, doi:10.1090/S0025-5718-1972-0314745-6,

as an improvement of methods of Serret and Hermite. For example, the
Brillhart method gives
$$
10^{50}+577=7611065343808354245450401^2+6486268906873921642245424^2.
$$
Both Serret and Hermite ideas are similar to the H.J.S. Smith's method
who in 1855 gave an elementary proof of Fermat's two square theorem
(using the notion of palindromic continuants); see

*

*F. W. Clarke, W. N. Everitt, L. L. Littlejohn and S. J. R. Vorster , H. J. S. Smith and the Fermat Two Squares Theorem, AMM 106 No. 7 (1999) 652–665, doi:10.2307/2589495, JSTOR
from where the above example was taken.
A: A recent article gives another very elegant proof of the two-square theorem in the Zagier style, but is easier.  See

*

*Stan Dolan, A very simple proof of the two-squares theorem,
The Mathematical Gazette, Volume 105, Issue 564 (2021) p. 511. doi:10.1017/mag.2021.120 (the whole article is visible in the free preview)

The proof is not one sentence but it is still well worth reading.  One starts with the set $S = \{(x,y,z) \mid (x+y+z)^2 = p + 4xz, \text{ and } 
x,y,z \in \mathbb{Z}_{>0} \}$.  Short, clever arguments show that $S$ is finite,
$|S|$ is odd, and that there is a positive integer solution with $y=z$.
A: As the answers above linked to an old paper of mine (in German, and a somewhat
different English preprint), some readers might like to know that an updated version is to appear very soon and is now linked on my webpage:
http://www.math.tugraz.at/~elsholtz/WWW/papers/papers30nathanson-new-address3.pdf
In addition to the motivation of the Heath-Brown/Zagier proof it contains for example
a) a discussion of a lattice point proof (section 1.6)
b) much more historical information and links to other work
c) an alternative motivation of the Heath-Brown-Zagier proof, due to
Dijkstra (section 2.3)
A: [This is more of a comment to Franz Lemmermeyer's answer.]
Let $F(x,y)=ax^2+bxy+cy^2=[a,b,c]$ be an integral ($a,b,c\in\mathbb{Z}$) binary quadratic form of discriminant $\Delta=b^2-4ac$, and let $\mathcal{Q}=\mathcal{Q}(\Delta)$ be the collection of all such forms.  $GL_2(\mathbb{Z})\times\langle i\rangle$ acts on $\mathcal{Q}$ by change of variable, preserving the discriminant
$$
g=\left(\begin{array}{cc}\alpha&\beta\\\gamma&\delta\end{array}\right), \ F^g=F(\alpha x+\beta y,\gamma x+\delta y).
$$
where
$$
i=\left(\begin{array}{cc}\sqrt{-1}&0\\0&\sqrt{-1}\\\end{array}\right)
$$
(which has determinant $-1$ and preserves integrality $[a,b,c]^i=[-a,-b,-c]$).
Consider the actions of
$$
s=\left(\begin{array}{cc}0&-\sqrt{-1}\\\sqrt{-1}&0\\\end{array}\right), \ [a,b,c]^s=[-c,b,-a], \ s^2=1
$$
and
$$
f=\left(\begin{array}{cc}1&1\\0&-1\\\end{array}\right), \ [a,b,c]^{f}=[a,2a-b,a-b+c], \ f^2=1.
$$
If $\mathcal{Q}$ contains a fixed point for the action of $s$, (i.e. $[a,b,-a]$), then $\Delta=b^2+(2a)^2$ is a sum of two squares, and if $\mathcal{Q}$ contains a fixed point for the action of $f$, (i.e. $[a,a,c]$), then $\Delta=a(a-4c)$ factors.
One interesting aspect of these involutions is that they go outside $SL_2(\mathbb{Z})$ to the broader orthogonal group $O(\Delta;\mathbb{Z})$ preserving the discriminant,
$$
\Delta(a,b,c)=\left(a \ b \ c\right)
\left(
\begin{array}{ccc}
0&0&-2\\
0&1&0\\
-2&0&0\\
\end{array}
\right)
\left(
\begin{array}{c}
a\\
b\\
c\\
\end{array}
\right).
$$
Other than that, they are among the ``simplest'' involutions available.
Now, if $\Delta=p>0$ is prime, then $f$ has fixed points only when $p\equiv 1\bmod 4$, namely
$$
\pm\left[1,1,\frac{1-p}{4}\right], \ \pm\left[p,p,\frac{p-1}{4}\right].
$$
One might hope to finagle a finite set of ``reduced'' forms that contains only one of these fixed points and is closed under both involutions, but this seems somewhat complicated (or impossible?).  [For instance, looking at $s$ with uniqueness in mind, one might consider those forms with $a>0>c$. Considering what happens with these under $f$, we get the condition $a+c<b$.  Looking at this condition under $s$, we get $|a+c|<b$ (so $b>0$) and $b>0$ implies $b<2a$ for stability under $f$, etc.]
The other two elements from $O(\Delta)$ (acting from the left on the column $(a,b,c)$ from my own predjudices) in the ``one-sentence'' proof aren't involutions; they're of order 6 and inverse to one another:
$$
\left(
\begin{array}{ccc}
0&0&-1\\
0&1&-2\\
-1&1&-1\\
\end{array}
\right)\leftrightarrow
\left(\begin{array}{cc}0&\sqrt{-1}\\-\sqrt{-1}&-\sqrt{-1}\\\end{array}\right)=t,
$$
$$
\left(
\begin{array}{ccc}
-1&1&-1\\
-2&1&0\\
-1&0&0\\
\end{array}
\right)\leftrightarrow
\left(\begin{array}{cc}\sqrt{-1}&-\sqrt{-1}\\-\sqrt{-1}&0\\\end{array}\right)=t^{-1}
$$
i.e.
\begin{align*}
[a,b,c]^t&=[-c,b-2c,-a+b-c],\\
[a,b,c]^{t^{-1}}&=[-a+b-c,-2a+b,-a],
\end{align*}
under the association
$$
F^g\longleftrightarrow
\left(
\begin{array}{ccc}
\alpha^2&\alpha\gamma&\gamma^2\\
2\alpha\beta&\alpha\delta+\beta\gamma&2\gamma\delta\\
\beta^2&\beta\delta&\delta^2\\
\end{array}
\right)
\left(
\begin{array}{c}
a\\
b\\
c\\
\end{array}
\right).
$$
Again, these are some of the simpler finite order $O(\Delta)$ elements (from the $2\times2$ matrix standpoint).
Now, what happens to the finite set of ``reduced'' forms
$$
\mathcal{Q}_0=\{[a,b,c]\in\mathcal{Q} : a>0>c, b>0\}
$$
under $t$ and $t^{-1}$?  We have
\begin{align*}
[a,b,c]^t&\in\mathcal{Q} \Longleftrightarrow b<a+c,\\
[a,b,c]^{t^{-1}}&\in\mathcal{Q} \Longleftrightarrow b>2a \text{ and } b>a+c.
\end{align*}
These conditions seem to start matching up with the the earlier attempt to find reduced forms stabalized by both $s$ and $f$.  So much so that the map
$$
[a,b,c]\mapsto\left\{
\begin{array}{cc}
[a,b,c]^t & b<a+c\\
[a,b,c]^f& a+c<b<2a\\
[a,b,c]^{t^{-1}} & b>2a\\
\end{array}
\right.
$$
is an involution on $\mathcal{Q}_0$ giving us what we want!
A: It's been a while since I read Elsholtz's article, but after doing so I felt none the wiser. Below I have translated Heath-Brown's proof into the language of binary quadratic forms; Zagier's proof looks more interesting from this point of view (the connections to Gauss reduction are much closer), but when working out the details I got stuck in the middle. 
One essential ingredient for the proofs by Heath-Brown and Zagier was pointed out already by Frick in 1918, who showed that if $p = a^2 + 4b^2$ is an odd prime number, then the indefinite binary quadratic form $Q = (-b,a,b)$ with discriminant $p$ is Gauss reduced and is contained in the principal cycle.
For proving that such a form exists without assuming that $p$ is a sum of two squares, we consider all forms $(A,B,C)$ with discriminant $p$ such that $A < 0$ and $C > 0$. From $p = B^2 - 4AC$ it then follows that the set
$$ S = \{(A,B,C): B^2 - 4AC = p, A < 0, C > 0\} $$
is finite. The obvious map
$$ \mu: S \to S, \quad (A,B,C) \to (-C,B,-A) $$
is an involution; if $S$ had odd cardinality, it would follow that $\mu$ has a fixed point, say $(A,B,-A)$, from which we would get $p = B^2 + 4A^2$. Unfortunately, $S$ has even cardinality since the involution
$$ \nu: S \to S, \quad (A,B,C) \to (A,-B,C) $$
has no fixed points: this is because $B = 0$ implies $p = 4AC$, which is impossible for prime numbers $p$.
We now would like to find a subset $U \subset S$ of $S$ with odd cardinality on which $\mu$ is still defined. The most natural idea would be considering the forms with $B > 0$. For showing that this set of forms has odd cardinality, we have to define an involution $(A,B,C) \to (A',B',C')$ on this subset that has exactly one fixed point. To find such an involution, we start with $(A,B,C) \to (A,-B,C)$ and then apply reduction by changing the middle coefficient modulo $2A$ and then adjusting the last coefficient so that the discriminant is $p$. This gives 
$$ (A,-B,C) \to (A',B',C') = (A,-2A-B,A+B+C). $$ 
Now we are facing the problem that it is not clear at all that $B' = -2A-B > 0$, or that $C' = A+B+C > 0$. But if we set
$$ U = \{(A,B,C) \in S: A+B+C > 0 \}, $$
then the map
$$ \gamma: (A,B,C) \to (A,-2A-B,A+B+C) $$ 
actually is an involution on $U$. Moreover, $(A,B,C)$ is a fixed point if and only if $-2A-B = B$ and $A+B+C = C$, which is equivalent to $A = -B$. Since $p = B^2 - 4AC = B^2 + 4BC = B(B+4C)$ is prime, we must have $|A| = |B| = 1$. Since $A < 0$, this implies that the fixed point is $(-1,1,\frac{p-1}4)$; this form is equivalent to the principal form $(1,1,\frac{p-1}4)$.
The involution $\gamma$ on $U$ shows that $U$ has odd cardinality; the map
$$ (A,B,C) \to (-C,-B,-A) $$
is an involution on $S$ sending $U$ to $S \setminus U$, which impliesthat $|S| = 2 |U|$. The involution $\nu$ on $S$ sends elements with $B > 0$ to elements with $B < 0$, hence 
$$ T = \{(A,B,C) \in S: B > 0\} $$ 
has the same number of elements as $U$, and in particular, it has odd cardinality. Finally, $\mu$ is an involution on $T$, and now the Two-Squares Theorem follows.
References


*

*H. Frick,
 Über den Zusammenhang der Perioden quadratischer Formen
 positiver Determinante mit der Zerlegung einer Zahl in die
 Summe zweier Quadrate, Diss. ETH Zürich, 1918

A: It is interesting that not only Zagier took this proof from Heath-Brown. Heath-Brown (according to his own words) took this proof from Uspensky. 
This trick has different applications, see articles of Bykovskii On the arithmetic nature of some identities of the elliptic functions theory  and The arithmetic nature of the triple and quintuple product identities .
A: Let me answer your question "where do these involution come from" with an elementary geometric explanation. You can skip ahead to the pictures, which are somewhat self-explanatory, I hope. 
The elements in the the $S$, i.e. triplets $(x,y,z)$ such that $x^2+4yz=p$, can be visualized as a square of side length $x$ together with $4$ rectangles  of size $y\times z$, where we place the rectangles such that for each corner of the square, a corner of a rectangle coincides with the corner of the square and a side of length $y$ coincides with the side of the square clockwise adjacent to the corner, and the interior of the rectangles and the square do not intersect. See the pictures below for examples. 
While it might seem trivial to display the number $x^2+4yz$ in this fashion, it very nicely illustrate Zagier's involution
 $$(x, y, z)\mapsto\begin{cases}(x+2z, z, y-x-z) &\text{if }  x < y-z\\(2y-x, y, x-y+z)&\text{if } y-z < x < 2y\\(x-2y, x-y+z, y)&\text{if } 2y < x
 \end{cases} .$$ 
In fact, an element of $S$, visualized in the way described above covers an area of a square (of side length $x+2z$) with $4$ rectangles removed. However, given such a square with oblong rectangles removed, there are precisely two ways of cutting it into a smaller configuration of a square and $4$ rectangles. Interchanging the two possible cuttings is precisely what Zagier's involution does on the non fixed points. The only fixpoint occurs if the area covered is a square with $4$ squares removed. For a prime number $p=1+4k$, this happens presicely once, namely for the configuration associated to $(x,y,z)=(1,1,k)$. 
We provide a complete example for $p=41$ (so $k=10$). 
Each picture shows two elements of $S$, which are mapped to each other under Zagier's involution and their common shape (in grey). 











The last image dispays the unique fixpoint of the involution. 
The other involution $(x,y,z)\mapsto (x,z,y)$ can of course also be visualized; the rectangles are simply rotated. 
If you have this simple illustration of Zagier's involution in mind, it is easy to recover the formula of it, which I find hard to remember otherwise. In this sense, I hope this answers where this involution comes from; although I am not sure, that Zagier was thinking in these images, when he wrote his one-liner. 
I heard about this illustration from Günter M. Ziegler, and as I understand it, Aigner and Ziegler plan to include it in the next edition of their "Proofs from the book". 
The first source seems to be some notes from a A. Spivak, which can be found here:
A. Spivak : Крылатые квадраты (Winged squares), Lecture notes for the mathematical circle at Moscow State University, 15th lecture 2007
I would be interested in knowing where this illustration appeared first, since it does not (yet) appear to be well known. 
