Special properties in the direct solving of sparse symmetric linear systems In the area of computational solving of large sparse linear systems, some solvers specialize only on symmetric sparse matrices, be it positive definite or indefinite as compared to general (non-symmetric) sparse systems solver. 
What mathematical properties does symmetric sparse matrix possess that makes it computationally more efficient to solve as compared to a nonsymmetric matrix?
 A: A symmetric matrix has a number of nice properties that solvers (both direct and iterative) take advantage of.  Since your question seems aimed at large, sparse systems, I presume you are also most interested in iterative solvers, namely, Krylov subspace methods.
For simplicity, suppose we have the system $A\textbf{x} = \textbf{b}$, and we wish to solve for $\textbf{x}$.  $A$ is large and sparse.
For a general $A$ and a general Krylov method, we generate an orthonormal basis $V_m$ spanning the space $\mathcal{K}_m(A,b) := \text{span} \{\textbf{b}, A\textbf{b}, A^2\textbf{b},... A^{m-1}\textbf{b}\}$.  The process of generating this basis is called the Arnoldi method, and it is essentially the modified Gram-Schmidt procedure applied to the basis $\{\textbf{b}, A\textbf{b},... A^{m-1}\textbf{b}\}$.  We obtain successive basis vectors $\textbf{v}_{m+1}$ as follows:
\begin{align}
\textbf{q}_{m+1} = A\textbf{v}_m - \sum_{k=1}^m (\textbf{v}^*_k A\textbf{v}_m) \textbf{v}_k
\\
\textbf{v}_{m+1} = \dfrac{\textbf{q}_{m+1}}{||\textbf{q}_{m+1}||}
\end{align}
This process can be summed up by the following relation, known as the Arnoldi relation:
$$
AV_m = V_m H_m + h_{m+1,m} \textbf{v}_{m+1},
$$
where $H_m$ is an upper Hessenberg matrix of dimension $m \times m$.
Notice that for a general $A$, we need access to all previous basis vectors $\textbf{v}_{k}$, $k \leq m$.  When $A$ is symmetric, however, the algebra reduces everything to a short-term recurrence-- meaning, we only need access to the previous two vectors in order to find the next orthonormal one.  This trick is summed up in the Lanczos procedure, and saves considerable storage.
But wait, there's more!  The point of Krylov methods is to reduce the large matrix $A$ to a low-rank approximation, namely the Hessenberg matrix $H_m$.  We then apply a direct method (via $QR$ or $LU$ factorizations) to solve the smaller approximate system generated by $H_m$.  (I'm leaving out details here, because they lead to different methods depending on choices made with the residual.)  For symmetric $A$, it can be shown that $H_m$ is actually tridiagonal.  The cost of directly solving a tridiagonal system is very low-- see, for example, the Thomas algorithm of complexity $O(n)$ and contrast with the typical $O(n^3)$ required for solvers based on $QR$ or $LU$ for general matrices.  Other methods may entail factoring $H_m$ instead of inverting it-- such factorizations (like Cholesky for symmetric positive definite matrices) require half the computational effort of their nonsymmetric counterparts ($LU$).
In summary, then, the algebraic relations arising from symmetric matrices reduce to short-term recursions and tridiagonal systems, and therefore less storage and less complexity.  A naive way to think of this is that since the matrix is symmetric, a computer only has to look at half of it-- thus cutting the workload down at least by half.
A: A great advantage of the sparse Cholesky factorization for symmetric and positive definite matrices is that you do not need to do pivoting for numerical stability but only focus on the symbolic diagonal pivoting to minimize fill-in. So you can completely separate the symbolic and numeric factorization and reuse the structure of the triangular factor when only the values of the original matrix change but not its structure. This is not the case for symmetric indefinite and nonsymmetric sparse direct solvers where you need to do pivoting not only to reduce fill-in but also for numerical stability reasons. Hence the symbolic and numeric factorization are not completely separable.
