Is there a "right" proof of Riemann's Theta Relation? Let $\theta$ denote the usual Jacobi Theta function (with auxiliary parameter $\tau = i$, for simplicity), i.e.
$$
\theta(z) = \sum_{n \in \mathbb{Z}} \exp(-\pi (a + n)^2 + 2 \pi i n z) \ .
$$
I'm interested in Riemann's quartic theta relation, which is (barring mistakes):
$$
\theta(x_1) \theta(y_1) \theta(u_1) \theta(v_1) = \frac{1}{2} \sum_{\eta \in \lbrace 0, 1/2, i/2, (1+i)/2 \rbrace } c_\eta \theta(x + \eta) \theta(y + \eta) \theta(u+\eta) \theta(v+\eta)
$$
where $c_\eta$ are some exponential factors, and $x_1$, $y_1$, $u_1$, $v_1$ are certain linear functions of the variables $x$, $y$, $u$, $v$ as follows:
$$
x_1 = (x + y +  u + v) / 2
$$
$$
y_1 = (x + y - u - v) / 2
$$
$$
u_1 = (x - y + u - v) / 2
$$
$$
v_1 = (x - y - u + v) / 2 \ .
$$
(This is taken from Mumford's "Tata lectures on theta I", chapter 1, section 5.)
There are many variations.  They're all fairly straightforward to prove with bare hands and the above formula; that is precisely the line taken in Mumford's book and every other reference I've seen.  However, I've been promised (somewhere) that every fact about special functions should have a "nice" interpretation coming from representation theory, and in particular from their interpretation as matrix coefficients (here, of the Heisenberg group).  Is there a clean proof along those lines?
To put it another way, if I only told you some properties of $\theta$ and its variants -- being an eigenfunction of certain operators, periodicity conditions and so on -- but not its formula, is there an enlightening reason to expect something like the Riemann theta formula to hold?
As usual, apologies if I've missed this in a standard reference.
 A: I think there is indeed a "right" proof of Riemann's quartic theta relation. Namely, as a corollary of Mumford's first and second "fundamental identities" (described in 6.4, 6.5 of Tata III). Mumford shows us that there is a simple machinery (thetas with quadratic forms) producing an infinite number of theta relations. This abundance of theta relations is not a totally obvious fact (at least not until one finds the fundamental identities).  
Riemann's quartic theta relation is apparently one very remarkable identity, or at least as remarkable as the $4\times 4$ rational symmetric positive definite and orthogonal matrix $$A=\frac{1}{2}\begin{pmatrix} 1 & 1 & 1 & 1 \cr 1 &1 &-1&-1 \cr 1&-1&1&-1 \cr 1&-1&-1&1 \end{pmatrix}.$$
Let us displace the question of "proving" riemann's quartic relation with confirming the identity, i.e. as a physical process of displacing and shuffling thetas. Then we should look maybe for a reason why $A$ is so remarkable. Or, at least see where else $A$ occurs.  
But let's agree on what a theta function is. This requires replacing the upper half plane $\mathbb{H}$ consisting of $z=x+iy$ with $y>0$ with Siegel's upper half space $\mathbb{H}_g$ consisting of symmetric $g\times g$ complex matrices $T$ whose imaginary part is positive definite. An element $T\in \mathbb{H}_g$ yields a hermitian form on $\mathbb{C}^g$, by defining $T[z]:={}^tzTz$. Likewise for any $g\times h$ matrix $N$ we set $T[N]={}^t N T N$. It's also useful for a ring $R$ to let $R(g,h)$ denote all $g\times h$ matrices over $R$. 
Computing from Mumford's second identity requires (explicit) coset representatives for the quotient $$\mathbb{Z}(g,4){}^t A / \mathbb{Z}(g,4){}^tA \cap \mathbb{Z}(g,4).$$ But the rational matrix $A$ is peculiar enough to have $A={}^tA=A^{-1}$ leading the coset representatives to take the form $(\eta, \eta, \eta, \eta)$, where $\eta$ varies over $\frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g$.
Identifying the coset representatives with the diagonal in $$\frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g \times \frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g  \times\frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g   \times\frac{1}{2}\mathbb{Z}^g / \mathbb{Z}^g  $$ is perhaps reason enough to affirm $A$'s singularity. The identification of this diagonal with these cosets may in fact have some other root system interpretation (in sense of Jeff Harvey's answer involving string things).   
But in a sense, after Mumford, it is the matrix $$A=\begin{pmatrix} 1 &1 \cr 1&-1 \cr \end{pmatrix}$$ which is responsible for riemann's quartic relation. This is Theorem 7.4 in Tata III. 
Now a comment regarding the OP's secondary question concerning, for lack of a better expression, a priori knowledge of holomorphic functions on the complex torus $\mathbb{C}^g / \mathbb{Z}^g + T \mathbb{Z}^g$, for $T \in \mathbb{H}_g$.
At the beginning of the theory, by whatever divine interference, we define a theta function $\theta:\mathbb{C}^g \times \mathbb{H}_g \to \mathbb{C}$ by the expression 
$$\theta(z,T)=\sum_{n \in \mathbb{Z}^g} \exp \pi i(T[n]+2{}^tnz).$$ In maintaining one's honesty (and i can only speak for myself) it is worth thinking about the necessity of our condition $im(T)>0$ in guaranteeing that $\theta(z,T)$ is a holomorphic function in $z$. One should also derive the periodicity and quasi-periodicity relations of $\theta(z,T)$ with respect to translations by  $\mathbb{Z}^g$ and $T\mathbb{Z}^g$ in $\mathbb{C}^g$.
From the definition we find ourselves with a holomorphic function on $\mathbb{C}^g$ periodic under translations by $\mathbb{Z}^g$ and satisfying the quasi-periodicity relation $$f(z+Tm) =\exp \pi i(-T[m] - 2{}^t mz) f(z).$$ The following observation is already contained on pp.121 of Tata I.
If $f:\mathbb{C}^g \to \mathbb{C}$ is any entire function satisfying the quasiperiodicity relations of $\theta(z,T)$, then $f= c\theta(z,T)$ for some constant $c\in \mathbb{C}$. Indeed $\mathbb{Z}^g$-periodicity permits a fourier expansion $$f(z)=\sum_{n\in \mathbb{Z}^g} c_n \exp 2\pi i {}^tn z.$$ Now quasiperiodicity permits us to derive recursive relations amongst the coefficients $c_n, \\ n\in \mathbb{Z}^g$. Explicitly, we can relate $c_{n}$ to $c_{n + \epsilon_k}$, where $\epsilon_k$ is the $k^{th}$ unit vector. It's amusing to think of how running around the lattice generated by $\mathbb{Z}^g T$ varies the cofficients in the fourier expanision.
So knowing that $f$ is a quasiperiodic holomorphic function already yields the formula for the theta function $\theta(z,T)$ -- whose usual formula we should recognize already as a fourier expansion. But I also know nothing about $\theta(z,T)$ as eigenfunction.
Now concerning "promises" relating the Heisenberg representations to facts on theta functions. From what I can tell, the only promise is Stone-von Neumann-Mackey's theorem on irreducible unitary representations (which are faithful on the centre) of the Heisenberg group: all such representations are equivalent. So given two representations (and there are many), there is an "intertwining" among them. But for some facts, one suffers. In 7.4 Tata III there is a description of riemann's quartic relation as arising from such an intertwining. However the Heisenberg group is defined over the finite adeles, i.e. $Heis(2g, \mathbb{A}_f)$. I would like to know where triality, spinors, or $so(8)$ interact in $\mathbb{A}_f$.
A: Let me consider one of the Riemann theta relations in Mumford, since as you say, once you know one of them it is straightforward to derive the rest. Relation (R5) is
$$ \theta_{00}(x) \theta_{00}(y) \theta_{00}(u) \theta_{00}(v) - \theta_{01}(x) \theta_{01}(y) \theta_{01}(u) \theta_{01}(v) - \theta_{10}(x) \theta_{10}(y) \theta_{10}(u) \theta_{10}(v) $$
$$+ \theta_{11}(x) \theta_{11}(y) \theta_{11}(u) \theta_{11}(v) = 2 \theta_{11}(x_1) \theta_{11}(y_1) \theta_{11}(u_1) \theta_{11}(v_1) $$
where the notation is that of Mumford with $\theta_{11}$ the odd theta function and the argument $\tau$ of the theta functions has been suppressed.  $x_1, y_1, u_1, v_1$ are defined in terms of $x,y,u,v$ as in the question. 
First, note that if we specialize to $x=y=u=v=0$ we obtain Jacobi's abstruse identity
$$ (\theta_{00}(0))^4 - (\theta_{01}(0))^4 - (\theta_{10}(0))^4=0 $$
This equation can be interpreted in terms of characters of level one characters for affine $D_4$. There are four
irreducible representations, basic, vector and two spinor representations. The above equation is the statement that
the vector and spinor characters are equal,  $ \chi_{vec}(\tau) = \chi_{sp}(\tau) $ and this is a corollary of the existence of triality for the $D_4$ Dynkin diagram.  In open string theory $D_4=SO(8)$ appears
as the transverse part of the ten-dimensional Lorentz group, space-time bosons are associate to the vector representation and space-time fermions to one of the spinor representations and the above is the statement that there are equal numbers of fermions and bosons at each mass level as required by supersymmetry.
There are actually two free fermion constructions of level one affine $D_4$. One of these goes by the name
of  the RNS construction in string theory while the other is known as the GS construction, they two constructions
are related by triality in $D_4$. For details see section 5.2 of Vol.1 of "Superstring Theory" by Green, Schwarz and Witten. The Riemann relation can be interpreted as the statement
that the RNS and GS formalisms give the same answer for the quantity 
$\chi_{vec}(\tau,x,y,u,v)-\chi_{sp}(\tau,x,y,u,v)$ where
$$ \chi_r(\tau,x,y,u,v)=Tr_{V_r} q^{L_0-c/24} \exp(2 \pi i (x H_1 + y H_2 + u H_3 + v H_4)) $$
and the $H_i$ are the zero modes of currents in the Cartan sub algebra of $SO(8)$  The left-hand side of the Riemann relation above is the result of computing this difference of characters in the RNS formalism while the right side is the result of the computation in the GS formalism. Apologies for my physics accent. I'm sure there are others who can translate this into more mathematically precise language if desired. 
