Closed-form maximum likelihood weights for linear regression

You observe $N$ pairs $(x_n,y_n)$, with $x_n\in\mathbb{R}^D$, and model each response as $y_n=x_n^\top\theta+\varepsilon_n$ with i.i.d. Gaussian noise $\varepsilon_n\sim N(0,\sigma^2)$, so responses are conditionally independent given their inputs. Stack the inputs as rows of a design matrix $X$ and the responses into a vector $y$. Rather than using gradient descent, derive a closed-form expression for the maximum likelihood estimate

Hint: work with the log-likelihood.

$\theta_{\text{MLE}}=(X^\top X)^{-1}X^\top y$. The log-likelihood is a constant minus $\tfrac{1}{2\sigma^2}\|y-X\theta\|^2$, so maximising it is least squares. Setting the gradient $X^\top(y-X\theta)$ to zero gives the normal equations $X^\top X\theta=X^\top y$, solvable when $X^\top X$ is invertible (full column rank).