---

Course: Math 535 - Mathematical Methods in Data Science (MMiDS)
Chapter: 2-Least squares: geometric, algebraic, and numerical aspects
Author: Sebastien Roch, Department of Mathematics, University of Wisconsin-Madison
Updated: July 15, 2024
Copyright: © 2024 Sebastien Roch

---

$\newcommand{\bmu}{\boldsymbol{\mu}}$ $\newcommand{\bSigma}{\boldsymbol{\Sigma}}$ $\newcommand{\bfbeta}{\boldsymbol{\beta}}$ $\newcommand{\bflambda}{\boldsymbol{\lambda}}$ $\newcommand{\bgamma}{\boldsymbol{\gamma}}$ $\newcommand{\bsigma}{{\boldsymbol{\sigma}}}$ $\newcommand{\bpi}{\boldsymbol{\pi}}$ $\newcommand{\btheta}{{\boldsymbol{\theta}}}$ $\newcommand{\bphi}{\boldsymbol{\phi}}$ $\newcommand{\balpha}{\boldsymbol{\alpha}}$ $\newcommand{\blambda}{\boldsymbol{\lambda}}$ $\renewcommand{\P}{\mathbb{P}}$ $\newcommand{\E}{\mathbb{E}}$ $\newcommand{\indep}{\perp\!\!\!\perp} \newcommand{\bx}{\mathbf{x}}$ $\newcommand{\bp}{\mathbf{p}}$ $\renewcommand{\bx}{\mathbf{x}}$ $\newcommand{\bX}{\mathbf{X}}$ $\newcommand{\by}{\mathbf{y}}$ $\newcommand{\bY}{\mathbf{Y}}$ $\newcommand{\bz}{\mathbf{z}}$ $\newcommand{\bZ}{\mathbf{Z}}$ $\newcommand{\bw}{\mathbf{w}}$ $\newcommand{\bW}{\mathbf{W}}$ $\newcommand{\bv}{\mathbf{v}}$ $\newcommand{\bV}{\mathbf{V}}$ $\newcommand{\bfg}{\mathbf{g}}$ $\newcommand{\bfh}{\mathbf{h}}$ $\newcommand{\horz}{\rule[.5ex]{2.5ex}{0.5pt}}$ $\renewcommand{\S}{\mathcal{S}}$ $\newcommand{\X}{\mathcal{X}}$ $\newcommand{\var}{\mathrm{Var}}$ $\newcommand{\pa}{\mathrm{pa}}$ $\newcommand{\Z}{\mathcal{Z}}$ $\newcommand{\bh}{\mathbf{h}}$ $\newcommand{\bb}{\mathbf{b}}$ $\newcommand{\bc}{\mathbf{c}}$ $\newcommand{\cE}{\mathcal{E}}$ $\newcommand{\cP}{\mathcal{P}}$ $\newcommand{\bbeta}{\boldsymbol{\beta}}$ $\newcommand{\bLambda}{\boldsymbol{\Lambda}}$ $\newcommand{\cov}{\mathrm{Cov}}$ $\newcommand{\bfk}{\mathbf{k}}$ $\newcommand{\idx}[1]{}$ $\newcommand{\xdi}{}$

Motivating example: predicting sales

The following dataset is from the excellent textbook [ISLP]. Quoting [ISLP, Section 2.1]:

Suppose that we are statistical consultants hired by a client to provide advice on how to improve sales of a particular product. The advertising data set consists of the sales of that product in 200 different markets, along with advertising budgets for the product in each of those markets for three different media: TV, radio, and newspaper. [...] It is not possible for our client to directly increase sales of the product. On the other hand, they can control the advertising expenditure in each of the three media. Therefore, if we determine that there is an association between advertising and sales, then we can instruct our client to adjust advertising budgets, thereby indirectly increasing sales. In other words, our goal is to develop an accurate model that can be used to predict sales on the basis of the three media budgets.

This a regression problem. That is, we want to estimate the relationship between an outcome variable and one or more predictors (or features). We load the data.

data = pd.read_csv('advertising.csv')
data.head()

	TV	radio	newspaper	sales
0	230.1	37.8	69.2	22.1
1	44.5	39.3	45.1	10.4
2	17.2	45.9	69.3	9.3
3	151.5	41.3	58.5	18.5
4	180.8	10.8	58.4	12.9

We will focus for now on the TV budget.

TV = data['TV'].to_numpy()
sales = data['sales'].to_numpy()

We make a scatterplot showing the relation between those two quantities.

plt.scatter(TV, sales, s=5, c='k')
plt.xlabel('TV'), plt.ylabel('sales')
plt.show()

There does seem to be a relationship between the two. Roughly a higher TV budget is linked to higher sales, although the correspondence is not perfect. To express the relationship more quantitatively, we seek a function $f$ such that

$$ y \approx f(\mathbf{x}) $$

where $\mathbf{x}$ denotes a sample TV budget and $y$ is the corresponding observed sales. We might posit for instance that there exists a true $f$ and that each observation is disrupted by some noise $\varepsilon$

$$ y = f(\mathbf{x}) + \varepsilon. $$

A natural way to estimate such an $f$ from data is $k$-nearest-neighbors ($k$-NN) regression$\idx{k-NN regression}\xdi$. Let the data be of the form $\{(\mathbf{x}_i, y_i)\}_{i=1}^n$, where in our case $\mathbf{x}_i$ is the TV budget of the $i$-th sample and $y_i$ is the corresponding sales. For each $\mathbf{x}$ (not necessarily in the data), we do the following:

• find the $k$ nearest $\mathbf{x}_i$'s to $\mathbf{x}$

• take an average of the corresponding $y_i$'s.

We implement this method in Python. We use the function numpy.argsort to sort an array and the function numpy.absolute to compute the absolute deviation. Our quick implementation here assumes that the $\mathbf{x}_i$'s are scalars.

def knnregression(x,y,k,xnew):
    n = len(x)
    closest = np.argsort([np.absolute(x[i]-xnew) for i in range(n)])
    return np.mean(y[closest[0:k]])

For $k=3$ and a grid of $1000$ points, we get the following approximation $\hat{f}$. Here the function numpy.linspace creates an array of equally spaced points.

k = 3
xgrid = np.linspace(TV.min(), TV.max(), num=1000)
yhat = [knnregression(TV,sales,k,xnew) for xnew in xgrid]

plt.scatter(TV, sales, s=5, c='b', alpha=0.5)
plt.plot(xgrid, yhat, 'r')
plt.xlabel('TV'), plt.ylabel('sales')
plt.show()

One downside of $k$-NN regression is that it does not give an easily interpretable relationship: if I increase my TV budget by $\Delta$ dollars, how is it expected to affect the sales? Another issue arises in high dimension where the counter-intuitive phenomena we discussed in a previous section can have a significant impact. Recall in particular the High-dimensional Cube Theorem. If we have $d$ predictors -- where $d$ is large -- and our data is distributed uniformly in a bounded region, then any given $\mathbf{x}$ will be far from any of our data points. In that case, the $y$-values of the closest $\mathbf{x}_i$'s may not be predictive. This is referred to as the Curse of Dimensionality$\idx{curse of dimensionality}\xdi$.

One way out is to make stronger assumptions on the function $f$. For instance, we can assume that the true relationship is (approximately) affine, that is, $y \approx \beta_0 + \beta_1 x$, or if we have $d$ predictors

$$ y \approx \beta_0 + \sum_{j=1}^d \beta_j x_j. $$

How do we estimate appropriate intercept and coefficients? The standard approach is to minimize the sum of the squared errors

$$ \sum_{i=1}^n \left(y_i - \left\{\beta_0 + \sum_{j=1}^d \beta_j (\mathbf{x}_{i})_j\right\}\right)^2, $$

where $(\mathbf{x}_{i})_j$ is the $j$-th entry of input vector $\mathbf{x}_i$ and $y_i$ is the corresponding $y$-value. This is called multiple linear regression.

It is a least-squares problem. We re-write it in a more convenient matrix form and combine $\beta_0$ with the other $\beta_i$'s by adding a dummy predictor to each sample. Let

$$ \mathbf{y} = \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}, \quad\quad A = \begin{pmatrix} 1 & \mathbf{x}_1^T \\ 1 & \mathbf{x}_2^T \\ \vdots & \vdots \\ 1 & \mathbf{x}_n^T \end{pmatrix} \quad\text{and}\quad \boldsymbol{\beta} = \begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_d \end{pmatrix}. $$

Then observe that

\begin{align*} \|\mathbf{y} - A \boldsymbol{\beta}\|^2 &= \sum_{i=1}^n \left(y_i - (A \boldsymbol{\beta})_i\right)^2\\ &= \sum_{i=1}^n \left(y_i - \left\{\beta_0 + \sum_{j=1}^d \beta_j (\mathbf{x}_{i})_j\right\}\right)^2. \end{align*}

The linear least-squares problem is then formulated as

$$ \min_{\boldsymbol{\beta}} \|\mathbf{y} - A \boldsymbol{\beta}\|^2. $$

In words, we are looking for a linear combination of the columns of $A$ that is closest to $\mathbf{y}$ in Euclidean distance. Indeed, minimizing the squared Euclidean distance is equivalent to minimizing its square root, as the latter in an increasing function.

One could solve this optimization problem through calculus (and we will come back to this approach later in the course), but understanding the geometric and algebraic structure of the problem turns out to provide powerful insights into its solution -- and that of many of problems in data science. It will also be an opportunity to review some basic linear-algebraic concepts along the way.

We will come back to the advertising dataset later in the chapter.

Background: review of vector spaces and matrix inverses

NUMERICAL CORNER: The plane $P$ made of all points $(x,y,z) \in \mathbb{R}^3$ that satisfy $z = x+y$ is a linear subspace. Indeed, $0 = 0 + 0$ so $(0,0,0) \in P$. And, for any $\mathbf{u}_1 = (x_1, y_1, z_1)$ and $\mathbf{u}_2 = (x_2, y_2, z_2)$ such that $z_1 = x_1 + y_1$ and $z_2 = x_2 + y_2$ and for any $\alpha \in \mathbb{R}$, we have

$$ \alpha z_1 + z_2 = \alpha (x_1 + y_1) + (x_2 + y_2) = (\alpha x_1 + x_2) + (\alpha y_1 + y_2). $$

That is, $\alpha \mathbf{u}_1 + \mathbf{u}_2$ satisfies the condition defining $P$ and therefore is itself in $P$. Note also that $P$ passes through the origin.

In this example, the linear subspace $P$ can be described alternatively as the collection of every vector of the form $(x, y, x+y)$.

We use plot_surface to plot it over a grid of points created using numpy.meshgrid.

x = np.linspace(0,1,num=101)
y = np.linspace(0,1,num=101)
X, Y = np.meshgrid(x, y)
print(X)

[[0.   0.01 0.02 ... 0.98 0.99 1.  ]
 [0.   0.01 0.02 ... 0.98 0.99 1.  ]
 [0.   0.01 0.02 ... 0.98 0.99 1.  ]
 ...
 [0.   0.01 0.02 ... 0.98 0.99 1.  ]
 [0.   0.01 0.02 ... 0.98 0.99 1.  ]
 [0.   0.01 0.02 ... 0.98 0.99 1.  ]]

print(Y)

[[0.   0.   0.   ... 0.   0.   0.  ]
 [0.01 0.01 0.01 ... 0.01 0.01 0.01]
 [0.02 0.02 0.02 ... 0.02 0.02 0.02]
 ...
 [0.98 0.98 0.98 ... 0.98 0.98 0.98]
 [0.99 0.99 0.99 ... 0.99 0.99 0.99]
 [1.   1.   1.   ... 1.   1.   1.  ]]

Z = X + Y
print(Z)

[[0.   0.01 0.02 ... 0.98 0.99 1.  ]
 [0.01 0.02 0.03 ... 0.99 1.   1.01]
 [0.02 0.03 0.04 ... 1.   1.01 1.02]
 ...
 [0.98 0.99 1.   ... 1.96 1.97 1.98]
 [0.99 1.   1.01 ... 1.97 1.98 1.99]
 [1.   1.01 1.02 ... 1.98 1.99 2.  ]]

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z, cmap='viridis')
plt.show()

Geometry of least squares: the orthogonal projection

NUMERICAL CORNER: To solve a linear system in Numpy, use numpy.linalg.solve. As an example, we consider the overdetermined system with

$$ A = \begin{pmatrix} 1 & 0\\ 0 & 1\\ 1 & 1 \end{pmatrix} \quad \text{and} \quad \mathbf{b} = \begin{pmatrix} 0\\ 0\\ 2 \end{pmatrix}. $$

We use numpy.ndarray.T for the transpose.

w1 = np.array([1., 0., 1.])
w2 = np.array([0., 1., 1.])
A = np.stack((w1, w2),axis=-1)
b = np.array([0., 0., 2.])
x = LA.solve(A.T @ A, A.T @ b)
print(x)

[0.66666667 0.66666667]

We can also use numpy.linalg.lstsq directly on the overdetermined system to compute the least-square solution.

x = LA.lstsq(A, b, rcond=None)[0]
print(x)

[0.66666667 0.66666667]

QR decomposition and Householder transformations

We implement the Gram-Schmidt algorithm in Python. For reasons that will become clear in the next subsection, we output both the $\mathbf{q}_j$'s and $r_{ij}$'s, each in matrix form. Here we use numpy.dot to compute inner products.

def gramschmidt(A):
    (n,m) = A.shape
    Q = np.zeros((n,m))
    R = np.zeros((m,m))
    for j in range(m):
        v = np.copy(A[:,j])
        for i in range(j):
            R[i,j] = np.dot(Q[:,i], A[:,j])
            v -= R[i,j]*Q[:,i]
        R[j,j] = LA.norm(v)
        Q[:,j] = v/R[j,j]
    return Q, R

NUMERICAL CORNER: Let's try a simple example.

w1 = np.array([1., 0., 1.])
w2 = np.array([0., 1., 1.])
A = np.stack((w1, w2),axis=-1)
print(A)

[[1. 0.]
 [0. 1.]
 [1. 1.]]

Q, R = gramschmidt(A)
print(Q)

[[ 0.70710678 -0.40824829]
 [ 0.          0.81649658]
 [ 0.70710678  0.40824829]]

print(R)

[[1.41421356 0.70710678]
 [0.         1.22474487]]

We implement back substitution in Python. In our naive implementation, we assume that the diagonal entries are not zero, which will suffice for our purposes.

def backsubs(R,b):
    m = b.shape[0]
    x = np.zeros(m)
    for i in reversed(range(m)):
        x[i] = (b[i] - np.dot(R[i,i+1:m],x[i+1:m]))/R[i,i]
    return x

Forward substitution is implemented similarly.

def forwardsubs(L,b):
    m = b.shape[0]
    x = np.zeros(m)
    for i in range(m):
        x[i] = (b[i] - np.dot(L[i,0:i],x[0:i]))/L[i,i]
    return x

We implement the QR approach to least squares.

def ls_by_qr(A, b):
    Q, R = gramschmidt(A)
    return backsubs(R, Q.T @ b)

NUMERICAL CORNER: We return to our simple overdetermined system example.

w1 = np.array([1., 0., 1.])
w2 = np.array([0., 1., 1.])
A = np.stack((w1, w2),axis=-1)
b = np.array([0., 0., 2.])

x = ls_by_qr(A, b)
print(x)

[0.66666667 0.66666667]

We implement the procedure above in Python. We will need the following function. For $\alpha \in \mathbb{R}$, let the sign of $\alpha$ be

$$ \mathrm{sign}(\alpha) = \begin{cases} 1 & \text{if $\alpha > 0$}\\ 0 & \text{if $\alpha = 0$}\\ -1 & \text{if $\alpha < 0$} \end{cases} $$

In Python, this is done using the function numpy.sign.

The following function constructs the upper triangular matrix $R$ by iteratively modifying the relevant block of $A$. On the other hand, computing the matrix $Q$ actually requires extra computational work that is often not needed. We saw that, in the context of the least-squares problem, we really only need to compute $Q^T \mathbf{b}$ for some input vector $\mathbf{b}$. This can be done at the same time that $R$ is constructed, as follows. The key point to note is that $Q^T \mathbf{b} = H_m \cdots H_1 \mathbf{b}$.

See here for an explanation of numpy.copy.

def householder(A, b):
    n, m = A.shape
    R = np.copy(A)
    Qtb = np.copy(b)
    for k in range(m):
    
        y = R[k:n,k]
        e1 = np.zeros(n-k)
        e1[0] = 1
        z = np.sign(y[0]) * LA.norm(y) * e1 + y
        z = z / LA.norm(z)
        
        R[k:n,k:m] = R[k:n,k:m] - 2 * np.outer(z, z) @ R[k:n,k:m]
        Qtb[k:n] = Qtb[k:n] - 2 * np.outer(z, z) @ Qtb[k:n]
    
    return R[0:m,0:m], Qtb[0:m]

In householder, we use both reflections defined above. We will not prove this here, but the particular choice made has good numerical properties. Quoting [TB, Lecture 10]:

Mathematically, either choice of sign is satisfactory. However, this is a case where numerical stability -- insensitivity to rounding errors -- dictates that one choice should be taken rather than the other. For numerical stability, it is desirable to reflect $\mathbf{x}$ to the vector $z \|\mathbf{x}\| \mathbf{e}_1$ that is not too close to $\mathbf{x}$ itself. [...] Suppose that [in the figure above] the angle between $H^+$ and the $\mathbf{e}_1$ axis is very small. Then the vector $\mathbf{v} = \|\mathbf{x}\| \mathbf{e}_1 - \mathbf{x}$ is much smaller than $\mathbf{x}$ or $\|\mathbf{x}\| \mathbf{e}_1$. Thus the calculation of $\mathbf{v}$ represents a subtraction of nearby quantities and will tend to suffer from cancellation errors.

NUMERICAL CORNER: We return to our overdetermined system example.

w1 = np.array([1., 0., 1.])
w2 = np.array([0., 1., 1.])
A = np.stack((w1, w2),axis=-1)
b = np.array([0., 0., 2.])
R, Qtb = householder(A, b)
x = backsubs(R, Qtb)
print(x)

[0.66666667 0.66666667]

One advantage of the Householder approach is that it produces a matrix $Q$ with very good orthogonality, i.e., $Q^T Q \approx I$. We give a quick example below comparing Gram-Schmidt and Householder. (The choice of matrix $A$ will become clearer when we discuss the singular value decomposition later in the chapter.)

seed = 535
rng = np.random.default_rng(seed)

n = 50
U, W = LA.qr(rng.normal(0,1,(n,n)))
V, W = LA.qr(rng.normal(0,1,(n,n)))
S = np.diag((1/2) ** np.arange(1,n+1))
A = U @ S @ V.T

Qgs, Rgs = gramschmidt(A)
print(LA.norm(A - Qgs @ Rgs))
print(LA.norm(Qgs.T @ Qgs - np.identity(n)))

1.4369568046009742e-16
19.745599060592102

As you can see above, the $Q$ and $R$ factors produced by the Gram-Schmidt algorithm do have the property that $QR \approx A$. However, $Q$ is far from orthogonal. (Recall that LA.norm computes the Frobenius norm introduced previously.)

On the other hand, Householder reflections perform much better in that respect as we show next. Here we use the implementation of Householder transformations in numpy.linalg.qr.

Qhh, Rhh = LA.qr(A)
print(LA.norm(A - Qhh @ Rhh))
print(LA.norm(Qhh.T @ Qhh - np.identity(n)))

4.739138228891714e-16
5.33506987519293e-15

Application: regression analysis

NUMERICAL CORNER: We test our least-squares method on simulated data. This has the advantage that we know the truth.

Suppose the truth is a linear function of one variable.

n, b0, b1 = 100, -1, 1
x = np.linspace(0,10,num=n)
y = b0 + b1*x

plt.scatter(x, y, s=3, c='k')
plt.show()

A perfect straight line is little too easy. So let's add some noise. That is, to each $y_i$ we add an independent random variable $\varepsilon_i$ with a standard Normal distribution (mean $0$, variance $1$).

seed = 535
rng = np.random.default_rng(seed)

y += rng.normal(0,1,n)

plt.scatter(x, y, s=5, c='k')
plt.show()

We form the matrix $A$ and use our least-squares code to solve for $\boldsymbol{\hat\beta}$. The function ls_by_qr, which we implemented previously, is in mmids.py, which is available on the GitHub of the book.

A = np.stack((np.ones(n),x),axis=-1)
coeff = mmids.ls_by_qr(A,y)
print(coeff)

[-1.03381171  1.01808039]

plt.scatter(x, y, s=5, c='b', alpha=0.5)
plt.plot(x, coeff[0]+coeff[1]*x, 'r')
plt.show()

NUMERICAL CORNER: Suppose the truth is in fact a degree-two polynomial of one variable with Gaussian noise.

n, b0, b1, b2 = 100, 0, 0, 1
x = np.linspace(0,10,num=n)
y = b0 + b1 * x + b2 * x**2 + 10*rng.normal(0,1,n)

plt.scatter(x, y, s=5, c='k')
plt.show()

We form the matrix $A$ and use our least-squares code to solve for $\boldsymbol{\hat\beta}$.

A = np.stack((np.ones(n), x, x**2), axis=-1)
coeff = mmids.ls_by_qr(A,y)
print(coeff)

[-2.76266982  1.01627798  0.93554204]

plt.scatter(x, y, s=5, c='b', alpha=0.5)
plt.plot(x, coeff[0] + coeff[1] * x + coeff[2] * x**2, 'r')
plt.show()

NUMERICAL CORNER: We return to the Advertising dataset from the [ISLP] textbook. We load the dataset again.

data = pd.read_csv('advertising.csv')

We will focus for now on the TV budget. We form the matrix $A$ and use our least-squares code to solve for $\boldsymbol{\beta}$.

TV = data['TV'].to_numpy()
sales = data['sales'].to_numpy()
n = np.size(TV)
A = np.stack((np.ones(n),TV),axis=-1)
coeff = mmids.ls_by_qr(A,sales)
print(coeff)

[7.03259355 0.04753664]

TVgrid = np.linspace(TV.min(), TV.max(), num=100)
plt.scatter(TV, sales, s=5, c='b', alpha=0.5)
plt.plot(TVgrid, coeff[0]+coeff[1]*TVgrid, 'r')
plt.show()

A degree-two polynomial might be a better fit.

A = np.stack((np.ones(n), TV, TV**2), axis=-1)
coeff = mmids.ls_by_qr(A,sales)
print(coeff)

[ 6.11412013e+00  6.72659270e-02 -6.84693373e-05]

plt.scatter(TV, sales, s=5, c='b', alpha=0.5)
plt.plot(TVgrid, coeff[0] + coeff[1] * TVgrid + coeff[2] * TVgrid**2, 'r')
plt.show()

The fit looks slightly better than the linear one. This is not entirely surprising though given that the linear model is a subset of the quadratic one. But, as we mentioned earlier, when adding more parameters we must now worry about overfitting the data. To illustrate, let's see what happens with a degree-$20$ polynomial fit.

deg = 20
A = np.stack([TV**i for i in range(deg+1)], axis=-1)
coeff = mmids.ls_by_qr(A,sales)
print(coeff)

[ 1.06538698e+00  6.72896471e-01 -1.53138969e-02 -2.74088516e-04
  1.83651714e-05 -3.40080020e-07  3.17915742e-09 -1.64042005e-11
  4.43633296e-14 -4.25654490e-17 -5.28727398e-20  1.11822932e-22
 -3.47096893e-27 -2.44665112e-30 -2.79435976e-33 -4.05263859e-36
 -6.83137511e-39 -1.27993830e-41 -2.59569760e-44 -5.59960687e-47
 -1.26949578e-49]

saleshat = np.sum([coeff[i] * TVgrid**i for i in range(deg+1)], axis=0)

plt.scatter(TV, sales, s=5, c='b', alpha=0.5)
plt.plot(TVgrid, saleshat, 'r')
plt.show()

The outcome now seems to vary wildly, seemingly driven by the randomness of the data.

CHAT & LEARN: Ask your favorite AI chatbot about using cross-validation to choose a suitable degree. Ask for code and apply it to this dataset. (Open In Colab) $\ddagger$