How does lme4 optimize the variance parameters of random effects (theta vector)? In mixed-effects models within the lmer R function of lme4 calculating the vector $\theta$ seems to be a key step, necessary to obtain the matrix $\Sigma(\theta)$, which is later decomposed as $\Sigma = TSST'$.
The parameter vector $\theta$ is defined as "random-effects parameters estimates; these are parameterized as the relative Cholesky factors of each random effect term".
They are generated through iteration, and get incorporated into the $\Sigma(\theta)$ matrix as in the following example:
In the sleepstudy dataset in lme4 the 'Reaction' time in miliseconds depends on the number of days of sleep deprivation, and the random effect of the actual individual being tested or "Subject."
Fitting a probably nonsensical mixed-effects model as fit <- lmer(Reaction ~ 1 + (1 | Subject), sleepstudy) would result in the following T, S and $\Sigma$ matrices in relation to $\theta$:
First off, let's call $\theta$ with getME(fit,"theta") # 0.807831. From it we can calculate manually,
str(sleepstudy$Subject) # Factor w/ 18 levels
Hence,
T <- diag(1,18).
S <- T * 0.807831, and finally Sigma <- Sigma <- S %*% t(S).
It is clear, then, that $\theta$ is key in the lmer function. It also sounds as though there is no algebraic formula for it. So the question is, How do you get $\theta$ in a step-by-step, easy to follow sequence of operations?
 A: Parameterizing variance-covariance matrices is actually rather challenging; in principle an $n \times n$ symmetric matrix can be described by $n(n+1)/2$ parameters, but the need for the matrix to be positive (semi)definite, and the desire to use a parameterization that is convenient for nonlinear optimization (e.g. parameters are not highly correlated with each other, possibly unconstrained, or at worst constrained in a simple way -- see e.g. Pinheiro and Bates 1996).
In the parameterization that lme4 uses, the $\theta$ vector represents the (columnwise unpacking of) the lower triangle of the Cholesky factor of the variance-covariance matrix $\Sigma$: for example, in the 2 $\times$ 2 random-slopes case we have
$$
\Sigma = \left(
\begin{array}{cc}
\theta_1 & 0 \\
\theta_2 & \theta_3 
\end{array}
\right)
\left(
\begin{array}{cc}
\theta_1 & \theta_2 \\
0 & \theta_3 
\end{array}
\right) =
\left(
\begin{array}{cc}
\theta_1^2 & \theta_1 \theta_2 \\
\theta_1 \theta_2 & \theta_2^2 + \theta_3^2 
\end{array}
\right) = 
\left(
\begin{array}{cc}
\sigma_1^2 & \sigma_{12} \\
\sigma_{12} & \sigma_2^2 
\end{array}
\right)
$$
So in this case, for example, we can equate $\theta_1$ with $\sigma_1$ (the standard deviation of the intercept), $\theta_2$ with $\sigma_{12}/\sigma_1$ (the covariance scaled by the SD) ... $\theta_3$ is a little more complicated.  If there is more than one RE term, the $\theta$ vector is the concatenation of the Cholesky factors for each term; in particular, in a model with only random-intercept terms, the $\theta$ vector is just the vector of RE standard deviations.
You're absolutely right that there is no formula -- not even a simple iterative process -- for determining $\theta$.  The way lme4 works is that, given a value of $\theta$, we can determine $\Sigma$ and hence the likelihood
$$
\cal L(\theta) = \int L(y|\boldsymbol b,\boldsymbol \beta) L(\boldsymbol b|\Sigma(\theta)) \, d \boldsymbol b
$$
As described in the lme4 paper, we can use some fairly clever algebra to reduce the integral above to a calculation that can be done in a few steps of linear algebra; hence, we can fairly efficiently compute the (negative log-)likelihood for a specified value of $\theta$.
After that we have to use brute force (i.e. nonlinear optimization); we pick a starting value for $\theta$ (typically 1 for the diagonal elements and 0 for the off-diagonal elements, which corresponds to starting the relative variance-covariance matrix [i.e. the variance-covariance matrix of the random effects scaled by the residual variance] as an identity matrix). Then we use some derivative-free algorithm, typically the Nelder-Mead simplex or Powell's BOBYQA algorithm, to adjust the value of $\theta$ until we think we've found a maximum-likelihood solution.
