Matrix representation of the OLS of an AR(1) process, Is there any precise way to express the OLS estimator of the centred error terms $\{u_t\}
_{t=1}^{n}$ that follows an AR(1) process? In other words, for
\begin{equation}
u_t=\rho u_{t-1}+\varepsilon_t,\quad \varepsilon_t\sim N(0,\sigma^2)
\end{equation}
is there a matrix representation for
\begin{equation}
\hat{\rho}=\frac{(1/n)\sum\limits_{t=1}^{n}u_tu_{t-1}}{(1/n)\sum\limits_{t=1}^{n}u_{t-1}^2}
\end{equation}
? I suspect there should be. However, I seem to fail to find it in Hamilton or other sources or derive an elegant expression myself.
Much appreciated in advance

Yes, there is an easy way. Expand the equations for every $n$.
\begin{equation}
u_2 = \rho u_1 + \epsilon_2 \\
u_3 = \rho u_2 + \epsilon_3 \\
\vdots \\
u_n = \rho u_{n-1} + \epsilon_n
\end{equation}
Let $\mathbf{y} = [u_2, \ldots, u_n]^T$ and $\mathbf{x} = [u_1, \ldots, u_{n-1}]^T$. Then,
\begin{align}
\mathbf{y} = \mathbf{x} \rho + \boldsymbol\epsilon
\end{align}
where $\boldsymbol\epsilon \sim \mathcal{N}(0, \sigma^2 \mathbf{I})$. Now, apply least squares.
\begin{align}
\hat{\rho} &= (\mathbf{x}^T \mathbf{x})^{-1} \mathbf{x}^T \mathbf{y} \\
&= \frac{\sum_{i=1}^{n-1} u_i u_{i+1}}{\sum_{i=1}^{n-1} u_i^2}
\end{align}
