  We present a fast direct algorithm for computing symmetric factorizations,
i.e. $A = WW^T$, of symmetric positive-definite hierarchical matrices with
weak-admissibility conditions. The computational cost for the symmetric
factorization scales as $\mathcal{O}(n \log^2 n)$ for hierarchically
off-diagonal low-rank matrices. Once this factorization is obtained, the cost
for inversion, application, and determinant computation scales as
$\mathcal{O}(n \log n)$. In particular, this allows for the near optimal
generation of correlated random variates in the case where $A$ is a covariance
matrix. This symmetric factorization algorithm depends on two key ingredients.
First, we present a novel symmetric factorization formula for low-rank updates
to the identity of the form $I+UKU^T$. This factorization can be computed in
$\mathcal{O}(n)$ time if the rank of the perturbation is sufficiently small.
Second, combining this formula with a recursive divide-and-conquer strategy,
near linear complexity symmetric factorizations for hierarchically structured
matrices can be obtained. We present numerical results for matrices relevant to
problems in probability \& statistics (Gaussian processes), interpolation
(Radial basis functions), and Brownian dynamics calculations in fluid mechanics
(the Rotne-Prager-Yamakawa tensor).
