Proof of correctness of normal equation

Question

I was taking an online course and saw linear regression being by gradient descent The intuition behind why the method would work seemed plausible.

I tried understanding normal equation as to why setting individual partial derivates to 0 and solving the equations give the optimal values of theta, but that didn't ring a bell.

1. Why setting partial derivatives to zero, and solving the equations gives optimal value of theta.
   I also went through the following link. The part about
   The minimum is determined by calculating the partial derivatives of S(β1,β2) with respect to β1 and β2 and setting them to zero
   is still dicey to me.
   Why would it work. Why would it give the values of β1 and β2 for minimized cost function.
2. If we have multiple local optima, would the normal equation method give the global optima? If so, why?

Answer 1

This is really just a matter of using standard results in multivariate calculus. In the ordinary least-squares (OLS) problem, the objective function being minimised is:

$$F(\boldsymbol{\beta}) = ||\mathbf{y} - \mathbf{x} \boldsymbol{\beta} ||^2.$$

The first and second derivatives (gradiant and Hessian) are:

$$\begin{equation} \begin{aligned} \nabla F(\boldsymbol{\beta}) &= 2 [ (\mathbf{x}^\text{T} \mathbf{x}) \boldsymbol{\beta} - (\mathbf{x}^\text{T} \mathbf{y}) ], \\[6pt] \nabla^2 F(\boldsymbol{\beta}) &= 2 (\mathbf{x}^\text{T} \mathbf{x}). \\[6pt] \end{aligned} \end{equation}$$

The Gramian matrix $\mathbf{x}^\text{T} \mathbf{x}$ is non-negative definite. So long as the columns of the design matrix $\mathbb{x}$ are linearly independent (which is generally the case in regression problems) it is positive definite. In the latter case the objective function $F$ is strictly convex, which means that it has a global minimising value at its unique critical point, which is found by solving:

$$\nabla F(\hat{\boldsymbol{\beta}}) = \mathbf{0} \quad \quad \implies \quad \quad \hat{\boldsymbol{\beta}} = (\mathbf{x}^\text{T} \mathbf{x})^{-1} (\mathbf{x}^\text{T} \mathbf{y}).$$

In the least-squares problem, linear independence of the columns of the design matrix $\mathbb{x}$ ensures that we do not have multiple optima. In the case where we have linearly dependence in the design matrix, this means that there are excess variables that give "ridges" of minimising values, and the parameters are not identified. This is usually dealt with by removing columns of the design matrix until linear independence is obtained.

Answer 2

Why setting partial derivatives to zero, and solving the equations gives optimal value of theta. I also went through the following link. The part about The minimum is determined by calculating the partial derivatives of S(β1,β2) with respect to β1 and β2 and setting them to zero is still dicey to me. Why would it work. Why would it give the values of β1 and β2 for minimized cost function.

From calculus, if we set the gradient to zero we arrive at an extremum (see, e.g, here). Hence, setting the partial derivatives to zero w.r.t. the $\beta_i$'s is generally a good approach to optimise over some objective function.

If we have multiple local optima, would the normal equation method give the global optima? If so, why?

We do not have multiple local optima. The loss function for the linear regression problem is quadratic and (strictly) convex, so there exists only one, unique, global optimum. Hence, by the properties of the objective function the normal equations will yield that global optimum.

Answer 3

I just want to complement and integrate the very good answers provided above, just to provide a different point of view and provide a different intuitive understanding of the mathematical nature of the solution.

Given an OLS problem like the typical $y=X \beta + \epsilon$ where y is a $Tx1$ vector, $X$ is a $TxN$ matrix where the number of samples is T and the number of regressors is N, $\beta$ and $\epsilon$ are $Nx1$ vectors, you want to minimize the sum of squared errors. The problem can be seen as the resolution of an overdetermined set of T equations (each corresponding to a row of the precious representation): each equation is $y_{t}=x_{t} \beta$ for each row t (notice that $x_{t} $ is now a vector representing the t-th row of the matrix $X$, $y_{t}$ is a scalar). Since the system is overdetermined for most real-world applications where $y$ is not a linear combination of the columns of $X$, then there exists no exact solutions to this set of equations. However the best proxy for a solution for an overdetermined system is to choose the minimum norm solution, which is obtained by setting $\beta$ equal to the product of the pseudoinverse of the matrix $X$ and the vector of dependent variable $y$, where the pseudoinverse of $X$ is indeed $(X^{T}X)^{-1}X^{T}$. This is indeed the solution to the system that minimizes the squared norm of the vector $\epsilon = y-X\beta$. This is why you have that formula for the OLS estimator.

This is a just different way of interpreting the OLS estimator compared to the usual minimization of squared error or maximization of the corresponding likelihood. Just a different way of looking at the same mathematical problem.

I also add this concise source