I just want to elucidate a bit on proofs provided by @roger and @acharuva.
Two popular ways to prove that a function $f$ is convex are to prove that:
\begin{equation}
f(tx_1+(1-t)x_2)-tf(x_1)-(1-t)f(x_2) \leq 0
\end{equation}
for $t \in [0,1]$ and to prove that the second derivative is non-negative for the entire support (domain of $f$). The first is the definition of convexity and the second is a well-known theorem (Proof here http://www.princeton.edu/~aaa/Public/Teaching/ORF523/S16/ORF523_S16_Lec7_gh.pdf)
.
For multivariate functions this means proving that the Hessian is PSD (A matrix A is PSD if A is symmetric and $u^TAu\geq0$ irrespective of $u$).
In case of OLS or the linear regression cost function we have:
\begin{equation}
J(\theta) = \frac 1 2 {(X \theta -Y)}^2 = \frac12( \theta^TX^TX\theta -2Y^TX\theta + Y^TY)
\end{equation}
\begin{equation}
\frac {\partial J(\theta)} {\partial \theta} = X^T(X\theta-Y)
\end{equation}
[because $\nabla_Atr(ABA^TC)=C^TAB^T+CAB$; here A=$\theta^T$, B=$X^TX$, C=$I$]
\begin{equation}
\frac {\partial^2 J(\theta)} {\partial \theta^2} = X^TX
\end{equation}
$X^TX$ is a PSD matrix ($u^TX^TXu=\Vert Xu\Vert^2$) and hence the cost function is convex with respect to $\theta$. That's the first proof.
The 2nd proof goes like this. It requires us to prove that
\begin{equation}
J(t\theta_1+(1-t)\theta_2)-tJ(\theta_1)-(1-t)J(\theta_2)\leq0
\end{equation}
$\forall \theta_1, \theta_2$ and $\forall t \in [0,1]$.
We will call this $E_1$. Note that here
\begin{equation}
J(\theta) = \frac12( \theta^TX^TX\theta -2Y^TX\theta + Y^TY)
\end{equation}
The full expansion of $E_1$ will be slightly verbose but we can make our life easy by considering that for any linear function $f(x)=ax+b$,
\begin{equation}
f(tx_1+(1-t)x_2)=tf(x_1)+(1-t)f(x_2)
\end{equation}
You can convince yourself of this by expanding both sides.
This means the term $\frac12( -2Y^TX\theta + Y^TY)$ of $J(\theta)$ will be $0$ in $E_1$. Let us now look at $E_1$ with only the first term of $J(\theta)$ and drop the $\frac12$ for convenience.
\begin{equation}
LHS = {(t\theta_1+(1-t)\theta_2)}^TX^TX(t\theta_1+(1-t)\theta_2)-t\theta_1^TX^TX\theta_1-(1-t)\theta_2^TX^TX\theta_2
\end{equation}
\begin{equation}
LHS = t^2\theta_1^TX^TX\theta_1+(1-t)^2\theta_2^TX^TX\theta_2+2t(1-t)\theta_1^TX^TX\theta_2 -t\theta_1^TX^TX\theta_1-(1-t)\theta_2^TX^TX\theta_2
\end{equation}
[Note that $\theta_1^TX^TX\theta_2 $=$(X\theta_1)^T(X\theta_2)$=$(X\theta_2)^T(X\theta_1)$=$\theta_2^TX^TX\theta_1$]
\begin{equation}
LHS = (t^2-t)\theta_1^TX^TX\theta_1+((1-t)^2-(1-t))\theta_2^TX^TX\theta_2+2t(1-t)\theta_1^TX^TX\theta_2
\end{equation}
\begin{equation}
LHS = -t(1-t)[\theta_1^TX^TX\theta_1+\theta_2^TX^TX\theta_2-2\theta_1^TX^TX\theta_2]
\end{equation}
[Note that $(1-t)^2-(1-t) = t^2-t=-t(1-t)$]
\begin{equation}
LHS = -t(1-t)[(\theta_1-\theta_2)^TX^TX(\theta_1-\theta_2)]
\end{equation}
\begin{equation}
LHS = -t(1-t) \Vert X(\theta_1-\theta_2) \Vert^2
\end{equation}
which is always $\leq0$.
This completes our second proof.
If you are having trouble following the second proof, try to first prove that a scalar quadratic function $f(x)=ax^2+bx+c$ is convex for $a>0$. I promise that exercise will be helpful (since our function for $J(\theta)$ is also quadratic in $\theta$).