Is a sample covariance matrix always symmetric and positive definite? When computing the covariance matrix of a sample, is one then guaranteed to get a symmetric and positive-definite matrix?
Currently my problem has a sample of 4600 observation vectors and 24 dimensions.
 A: For a sample of vectors $x_i=(x_{i1},\dots,x_{ik})^\top$, with $i=1,\dots,n$, the sample mean vector is 
$$
  \bar{x}=\frac{1}{n} \sum_{i=1}^n x_i \, ,
$$ and the sample covariance matrix is
$$
  Q = \frac{1}{n} \sum_{i=1}^n (x_i-\bar{x})(x_i-\bar{x})^\top \, .
$$
For a nonzero vector $y\in\mathbb{R}^k$, we have
$$
  y^\top Qy = y^\top\left(\frac{1}{n} \sum_{i=1}^n (x_i-\bar{x})(x_i-\bar{x})^\top\right) y
$$
$$
 = \frac{1}{n} \sum_{i=1}^n y^\top (x_i-\bar{x})(x_i-\bar{x})^\top y
$$
$$
  = \frac{1}{n} \sum_{i=1}^n \left( (x_i-\bar{x})^\top y \right)^2 \geq 0 \, . \quad (*)
$$
Therefore, $Q$ is always positive semi-definite.
The additional condition for $Q$ to be positive definite was given in whuber's comment bellow. It goes as follows.
Define $z_i=(x_i-\bar{x})$, for $i=1,\dots,n$. For any nonzero $y\in\mathbb{R}^k$, $(*)$ is zero if and only if $z_i^\top y=0$, for each $i=1,\dots,n$. Suppose the set $\{z_1,\dots,z_n\}$ spans $\mathbb{R}^k$. Then, there are real numbers $\alpha_1,\dots,\alpha_n$ such that $y=\alpha_1 z_1 +\dots+\alpha_n z_n$. But then we have $y^\top y=\alpha_1 z_1^\top y + \dots +\alpha_n z_n^\top y=0$, yielding that $y=0$, a contradiction. Hence, if the $z_i$'s span $\mathbb{R}^k$, then $Q$ is positive definite. This condition is equivalent to $\mathrm{rank} [z_1 \dots z_n] = k$.
A: @Zen's answer plus @whuber's comment to @Konstantin's answer provide a complete proof. Nevertheless, I'll rephrase the proof by trying to place more statistical emphasis.
Indeed, one can say that the sample covariance matrix $S$ is always positive and semi-definite because it can be seen as the variance of a suitable univariate variable, which is always non-negative.
In detail, let $x_1,\ldots,x_n$ be the observed sample, with $x_i = (x_{i1},\ldots,x_{ik})^\top$, $i=1,\ldots,n$. The sample covariance matrix is
$$
Q = n^{-1}\sum_{i=1}^n(x_i-\bar x)(x_i-\bar x)^\top,
$$
where $\bar x=n^{-1}\sum_{i}x_i$ is the sample average.
Consider now any vector $a = (a_1,\ldots,a_k)^\top$ and take the $y_i$, linear combination of $x_i$ with coefficients $a_i$, i.e.
$$
y_i = a^\top x_i = a_1x_{11}+\cdots+a_{k}x_{ik},\quad\text{for all } i.
$$
Let $\bar y$ be the sample average of $y_i$'s and note that $\bar y = a^\top \bar x$. The variance of $y_i$ is
\begin{align*}
0\leq s_{y}^2 &= n^{-1}\sum_i(y_i-\bar y)^2 = n^{-1}\sum_{i}(y_i-\bar y)(y_i-\bar y)^\top\\
& = n^{-1}\sum_{i} (a^\top x_i - a^\top \bar x)(a^\top x_i - a^\top \bar x)\\
& = a^\top\left(n^{-1}\sum_{i} (x_i - a^\top \bar x)(x_i -\bar x)\right)a\\
& = a^\top S a.
\end{align*}
Since $a$ was arbitrary, this completes the proof.
A: Let
$$
X=
\begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1k} \\
x_{21} & x_{22} & \cdots & x_{2k} \\
\vdots & \vdots & \ddots & \vdots\\
x_{n1} & x_{n2} & \cdots & x_{nk}
\end{pmatrix}
$$
denote the data matrix whose $\left(i,j\right)$-th entry is the $i$-th measurement of the $j$-th variable (with $i \in \{1,\ldots, n\}, j \in \{1,\ldots,k \}$).

The sample covariance matrix $\mathcal S$ can be written as
$\mathcal S=n^{-1}X^\top  C_n X,$
where $C_n=I_n-n^{-1}\mathbb{1}_n\mathbb{1}_n^\top$ is the centering matrix.
Since $C_n$ is symmetric and idempotent, we also have $\mathcal S=n^{-1}X^\top C_n^\top  C_n X$.1 But with $Y\mathrel{:=}C_n X$ this becomes $\mathcal S=n^{-1}Y^\top Y$, which is generally positive semi-definite, and positive definite only if the columns of $Y$ are linearly independent.
This means that $\mathcal S$ is positive definite iff the centered measurement vectors of the $k$ variables, i.e. the vectors $\left(x_{1j}-\bar{x}_{.j},\ldots,x_{nj}-\bar{x}_{.j}\right)^\top$ indexed by $j$, are linearly independent.

1Another way to see that $\mathcal S$ can be written as $n^{-1}X^\top C_n^\top  C_n X$ is to interpret $X^\top C_n^\top  C_n X = \left(C_n X\right)^\top  \left(C_n X\right)$ as sum of outer products of the rows of the column-wise centered $X$ with itself.
A: A correct covariance matrix is always symmetric and positive *semi*definite.
The covariance between two variables is defied as $\sigma(x,y) = E [(x-E(x))(y-E(y))]$.
This equation doesn't change if you switch the positions of $x$ and $y$. Hence the matrix has to be symmetric.
It also has to be positive *semi-*definite because:
You can always find a transformation of your variables in a way that the covariance-matrix becomes diagonal. On the diagonal, you find the variances of your transformed variables which are either zero or positive, it is easy to see that this makes the transformed matrix positive semidefinite. However, since the definition of definity is transformation-invariant, it follows that the covariance-matrix is positive semidefinite in any chosen coordinate system.
When you estimate your covariance matrix (that is, when you calculate your sample covariance) with the formula you stated above, it will obv. still be symmetric. 
It also has to be positive semidefinite (I think), because for each sample, the pdf that gives each sample point equal probability has the sample covariance as its covariance (somebody please verify this), so everything stated above still applies.
A: I would add to the nice argument of Zen the following which explains why we often say that the covariance matrix is positive definite if $n-1\geq k$.
If $x_1,x_2,...,x_n$ are a random sample of a continuous probability distribution then $x_1,x_2,...,x_n$ are almost surely (in the probability theory sense) linearly independent. 
Now, $z_1,z_2,...,z_n$ are not linearly independent because $\sum_{i=1}^n z_i = 0$, but because of $x_1,x_2,...,x_n$ being a.s. linearly independent, $z_1,z_2,...,z_n$ a.s. span $\mathbb{R}^{n-1}$. If $n-1\geq k$, they also span $\mathbb{R}^k$. 
To conclude, if $x_1,x_2,...,x_n$ are a random sample of a continuous probability distribution and $n-1\geq k$, the covariance matrix is positive definite.
