If you’ve ever wondered how R or python gives you regression coefficients, the answer is linear algebra! Linear algebra operations are essential to almost all modern methods for analyzing or modeling data. In this post I will first go through through some linear algebra basics, then build up to finding the ordinary least squares estimate for multiple linear regression. To try and make this self-contained, I will start from the basics. For the linear regression part, I will just assume knowledge of basic derivative rules, e.g., finding the derivative of \((a - bx)^2\) with respect to \(x\).

This post is about how to get linear regression coefficients. If you want to know what to do with them or how to interpret them afterwards, check out Kevin’s great tutorial!

Building up to matrices

We begin with the building block of scalars, which here will refer to any real number. A variable referring to a scalar will be a lower-case letter such as \(x\). For example, \(x\) can refer to the recording of a neuroimaging sensor at a particular time, or in general the value of a sample’s measurement.

We can arrange a list of scalars into a vector. We can enumerate the entries of a vector through a subscript, e.g., \(x_1\) refers to the first entry of a vector. Vectors will be denoted with a bold-faced lower-case letter, \(\mathbf{x}\), and can be visualized in the following way:

Where we have a list of \(n\) scalars. The number of scalars in a vector \(n\) is referred to as the dimensionality of the vector. A more mathematically formal way of stating this is that \(\mathbf{x} \in \mathbb{R}^n\), which states that the vector \(\mathbf{x}\) belongs to a set of all (real) vectors with dimension \(n\).

Importantly, vectors can be thought of as an arrow in \(n\)-dimensional space. The arrow goes from the origin of the space to the coordinates listed in the vector. A vector’s \(L_2\) norm or length is denoted by:

Which is the Euclidean distance from the origin to the vector tip.

(An aside: As human beings, we’re mostly limited to reasoning about 2-dimensional and 3-dimensional spaces. However, the data we collect is often higher-dimensional. The wonder of linear algebra is that it allows us us to reason and have intuitions about much higher-dimensional vectors and spaces.)

Basic vector operations

Some basic things we can do with a vector are:

• multiply the vector with a scalar
• add the vector to another vector
• find a scalar value called the dot product or inner product of the vector with another vector
• (there’s also an outer product, but we’ll come back to that in a bit)

Multiplying a vector with a scalar produces another vector, and can be thought of as scaling the vector in space (it points in the same direction but is larger or smaller):

We can add vectors, provided they have the same dimension:

We can define a dot product between two same-sized vectors, which is a scalar that indicates something about the angle between the vectors (\(\theta\)) in space.

Thus, if two vectors are oriented at \(90^\circ\) relative to each other - that is, they are orthogonal - their dot product is \(0\). This last formulation allows us to define the concept of angles in higher than 3-dimensional space.

The dot product of a vector with itself is its length squared:

Matrices

A list of vectors of the same dimension can form another mathematical object called a matrix, commonly denoted as upper-case bold letters, such as \(\mathbf{X}\):

In most data matrices, \(x_{n,p}\) is a scalar variable indicating the measurement of the \(p^{th}\) dimension of the \(n^{th}\) sample. As with vectors, we can formally state \(\textbf{X} \in \mathbb{R}^{n \times p}\), indicating that the matrix \(\textbf{X}\) has \(n\) rows and \(p\) columns. A more common way to state this is that the shape of \(\mathbf{X}\) is \((n, p)\).

We can think of matrices as a list of column vectors. For example, we could think of our data matrix as:

Here we have a list of \(p\) vectors \(\mathbf{x}_1, \dots, \mathbf{x}_p\), all of which are \(n\)-dimensional. In other words, each vector \(\mathbf{x}_i\) represents all of the samples’ data for whatever the \(i^{\textrm{th}}\) dimension was measuring.

We could also think of a matrix as a list of row vectors. A row of \(\mathbf{X}\) indicates all of the dimensions of a particular sample:

Here we have \(n\) row vectors (they are different than column vectors), which are denoted as \(\mathbf{x}_1^\top, \dots, \mathbf{x_n}^\top\). These \(n\) vectors represent one sample each, and are of dimension \(p\) which indicates that there were \(p\) dimensions measured.

If you’re wondering where the \(^\top\) symbol comes from, it’s actually the transpose of a matrix. It is denoted by a superscript \(\mathbf{X}^\top \in \mathbb{R}^{p \times n}\) and is an operation which turns each of the rows of \(\mathbf{X}\) into columns of \(\mathbf{X}^\top\). In other words, if \(\mathbf{X}\) is of shape \((n,p)\), then \(\mathbf{X}^\top\) will be of shape \((p,n)\):

Here we now have dimensions measured as columns and samples as rows. In essence, row vectors are transposes of column vectors.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
# working with vectors and matrices in numpy
import numpy as np 

# defining vectors (1d array)
x = np.array([0, 1, 2, 3])
y = np.array([1, 2, 3, 4])

# somewhat confusingly, vectors (1d arrays) in numpy are printed horizontally
print('x is:\t', x)
print('y is:\t', y)
print('x shape:', x.shape)
print('---------')

# operations
print('5x = \t', 5*x)
print('x + y =\t', x+y)
print('x * y =\t', x.dot(y))
print('---------')

# defining a matrix (2d array)
n, p = 5, 3 # 5 samples, 3 dims
X = np.random.randint(10, size=(n,p))   # random integer matrix
print('X looks like:\n', X)
print('X is of shape: ', X.shape)
print('---------')

# transpose
print('X transpose:\n', X.T)            # numpy syntax for transpose
print('X transpose shape: ', X.T.shape)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
x is:	 [0 1 2 3]
y is:	 [1 2 3 4]
x shape: (4,)
---------
5x = 	 [ 0  5 10 15]
x + y =	 [1 3 5 7]
x * y =	 20
---------
X looks like:
 [[4 9 6]
 [2 1 7]
 [0 1 3]
 [0 9 9]
 [0 5 6]]
X is of shape:  (5, 3)
---------
X transpose:
 [[4 2 0 0 0]
 [9 1 1 9 5]
 [6 7 3 9 6]]
X transpose shape:  (3, 5)

Matrix-vector multiplication

Multiplying a matrix with a vector results in a new vector. A way to intuit matrix-vector multiplication is to consider the matrix as a a “machine” which acts on the input vector and produces an output vector of a particular dimension. In this case, the matrix must have the dimensionality of the row vectors be equal to the dimensionality of the input vector.

There are two ways to frame matrix-vector multiplications, which result from thinking about the matrix as a either a collection of column vectors or row vectors.

• We can get the output vector by weighting the columns of the matrix with the entries of the input vector and summing these weighted vectors.
• Alternatively, we can get each element of the output vector by doing the dot product of the rows of the matrix with the input vector.

Both perspectives will be useful for the rest of this post. Either way, if we have a size \((n, p)\) matrix, that means we can do matrix-vector multiplication by giving it a \(p\)-dimensional vector, and our output will be an \(n\)-dimensional vector. To illustrate, let’s define a \((2 \times 3)\) matrix \(\mathbf{A}\) and a \(3\)-dimensional vector \(\mathbf{b}\):

The different colors (blue, orange, green) are different columns, and darker/lighter are different rows. Now, we are looking to find the output vector \(\mathbf{c = Ab}\):

How do we actually multiply these things together? Let’s see the first way of framing matrix-vector multiplication, which is to consider weighting the columns of \(\mathbf{A}\) with the elements of \(\mathbf{b}\):

For this case, imagine taking the elements of \(\mathbf{b}\) and putting them in front of the columns of \(\mathbf{A}\) as scalar multiples, then summing up those scaled columns. In some sense, the columns of matrix \(\mathbf{A}\) provide a set of basis patterns which are weighted and summed to produce \(\mathbf{c}\).

If we do a little rearranging of the above definition, we can get the other way of framing matrix-vector multiplication:

Where

That is, the elements of \(\mathbf{c}\) are the dot products of the row vectors of \(\mathbf{A}\) with \(\mathbf{b}\). For this case, imagine taking the vector \(\mathbf{b}\) and putting it sideways on top of the matrix \(\mathbf{A}\). Then, the first element of \(\mathbf{c}\) is the dot product of \(\mathbf{b}\) with the first row of \(\mathbf{A}\), the second element of \(\mathbf{c}\) is the dot product of \(\mathbf{b}\) with the second row of \(\mathbf{A}\), and so on.

To recap, and in general for multiplying an \((n,p)\) matrix and a \(p\) dimensional vector:

Either way, the output is an \(n\)-dimensional vector.

Vectors and vector operations as 1-dimensional matrices

Imagine we have a single row vector \(\mathbf{a}^\top\) and a column vector \(\mathbf{b}\), both of dimension \(n\). If we consider the vector \(\mathbf{a}^\top\) to be a size \((1, n)\) matrix and \(\mathbf{b}\) to be an \(n\) dimensional vector, we can see that the matrix-vector product \(\mathbf{a}^\top\mathbf{b}\) should be a size \((1 \times 1)\) vector, otherwise known as a scalar. As it turns out, we’ve already seen this before as the dot product (or inner product):

Let’s now consider the outer product, sometimes denoted \(\mathbf{a} \otimes \mathbf{b}\):

What might it mean to multiply an \((n,1)\) matrix with a \((1,n)\) matrix?

Matrix-matrix multiplication

We can now define how to multiply two matrices \(\mathbf{A}\) and \(\mathbf{B}\) together to get a matrix \(\mathbf{C = AB}\). In fact, matrix-matrix multiplication is equivalent to doing a series of matrix-vector multiplications. That is, imagine \(\mathbf{B}\) as a series of column vectors \(\mathbf{b}_1, \cdots, \mathbf{b}_k\). Then, each column vector \(\mathbf{c}_i\) will be equal to \(\mathbf{Ab}_i\):

Thus, the first column of \(\mathbf{C}\) is simply \(\mathbf{Ab}_1\). This reveals one rule of matrix-matrix multiplication: the so-called inner dimensions of the matrices must agree in order for the multiplication to be defined. That is, \(\mathbf{C = AB}\) is only defined when the matrix \(\mathbf{A}\) has shape \((n, k)\) and the matrix \(\mathbf{B}\) has shape \((k, p)\) for any \(n\) and \(p\). In other words, the rows of \(\mathbf{A}\) must be of the same dimension as the columns of \(\mathbf{B}\). Then, the output \(\mathbf{C}\) will be of shape \((n,p)\).

Importantly, this tells us that matrix multiplication is not commutative. That is, \(\mathbf{AB} = \mathbf{BA}\) is not generally true. In fact, if \(\mathbf{B}\) is size \((k,p)\) and \(\mathbf{A}\) is size \((n,k)\), then \(\mathbf{BA}\) isn’t even defined, as the inner dimensions do not agree. Thus, left multiplying a matrix is fundamentally different than right multiplying a matrix.

The identity matrix and matrix inverses

An identity matrix is a matrix with ones along the diagonal and zeros everywhere else. For example, size \(2\) identity matrix is defined as:

In general:

We can think of the identity matrix as the matrix version of a 1 for matrix multiplication. In other words, \(1*a = a\). Similarly, \(\mathbf{Ib} = \mathbf{b}\), and \(\mathbf{IA} = \mathbf{AI} = \mathbf{A}\).

The inverse of a matrix only exists for square matrices (those with size \((n,n)\)), and only certain square matrices at that. An inverse \(\mathbf{A}^{-1}\) of a matrix \(\mathbf{A}\) is defined such that \(\mathbf{A}^{-1} \mathbf{A} = \mathbf{A} \mathbf{A}^{-1} = \mathbf{I}\). If we have an equation like \(\mathbf{Ax} = \mathbf{b}\), and we wish to get \(\mathbf{x}\) alone, we can left multiply by \(\mathbf{A}^{-1}\):

1
2
3
4
5
6
7
8
9
10
11
# matrix-vector 
n, p = 2, 3                       # 2 rows, 3 cols
A = np.arange(n*p).reshape((n,p)) # fast way to get 2d array with ordered ints
b = np.arange(p)                  # vectors don't have a second dimension in numpy
c = A @ b                         # syntax for matrix multiplication in numpy

print('A:\n', A)
print('------')
print('b:\n', b)   # remember b is a column vector, despite how it's printed
print('------')
print('c:\n', c)   # same with c

1
2
3
4
5
6
7
8
9
A:
 [[0 1 2]
 [3 4 5]]
------
b:
 [0 1 2]
------
c:
 [ 5 14]

1
2
3
# let's try multiplying A with a 4-dim vector for fun
x = np.arange(4)
c = A.dot(x)  # another syntax for matrix multiplication in numpy

1
2
3
4
5
6
7
8
9
10
11
---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

<ipython-input-58-aa5665428299> in <module>()
      1 # let's try multiplying A with a 4-dim vector for fun
      2 x = np.arange(4)
----> 3 c = A.dot(x)  # another syntax for matrix multiplication in numpy


ValueError: shapes (5,3) and (4,) not aligned: 3 (dim 1) != 4 (dim 0)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# matrix-matrix
n, p, k = 5, 3, 2
A = np.arange(n*p).reshape((n,p))
B = np.arange(p*k).reshape((p,k))
C = A @ B

print('A shape:', A.shape, '\n', A)
print('------')
print('B shape:', B.shape, '\n', B)
print('------')
print('C shape:', C.shape, '\n', C)
print('------')

# let's get first col of C:
b_1 = B[:, 0]
c_1 = A.dot(b_1)
print('c_1:\n', c_1)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
A shape: (5, 3) 
 [[ 0  1  2]
 [ 3  4  5]
 [ 6  7  8]
 [ 9 10 11]
 [12 13 14]]
------
B shape: (3, 2) 
 [[0 1]
 [2 3]
 [4 5]]
------
C shape: (5, 2) 
 [[ 10  13]
 [ 28  40]
 [ 46  67]
 [ 64  94]
 [ 82 121]]
------
c_1:
 [10 28 46 64 82]

Estimating linear regression coefficients with matrix/vector operations

Let’s put all the above to work to help us do linear regression.

We are all familiar with the most simplified univariate model: \(y = mx + b\). To generalize this, we can consider more than one independent (input) variable by considering a model where we have a vector of \(n\) observations \(\mathbf{x}\) instead of a single observation \(x\) and a vector of \(n\) weights \(\mathbf{b}\) instead of the single weight \(m\). For instance, consider the model prediction for a single observation of \(n\) regressors \(x_1, x_2, \dots, x_p\):

There’s a lot of indices in that equation! We can simplify the notation using vectors, as this model is defined as the dot product of the regressors and the weights. We can also include the intercept \(b_0\) in the vector of weights \(\mathbf{b}\) by adding a single \(1\) to the front of \(\mathbf{x}_i\), so that both vectors are of size \((p+1)\):

Where

If we want to estimate the coefficients \(\mathbf{b}\) we can do a bit more linear algebra and some calculus to find the solution as an optimization problem. More specifically, we first have to define the objective, i.e., the thing we want to minimize. It stands to reason that we want to minimize the error of our model predictions over all our data. Then, we take the derivative of the error with respect to the parameters we want to find (the coefficients), and set that expression equal to 0. Finally, we find the values of \(\mathbf{b}\) which solve those derivative equations.

To begin, let’s define the error for sample \(i\) as:

Now we need to find values for our coefficients \(\mathbf{b}\) that minimize this loss over all samples. For a few reasons, we will minimize the sum of squared errors:

We could continue solving this optimization problem by taking the derivative with respect to \(b_0, \dots, b_p\), setting it to 0, and solving for the values of \(\mathbf{b}\). It would get a bit ugly with all the algebra, so the vector and matrix notation can help us out a bit here by stacking all of our \(n\) observations. Define the \((n, 1)\) vector \(\mathbf{y}\) as the vector of all the observed responses, the \((n, p+1)\) matrix \(\mathbf{X}\) as the design matrix of observed regressors, and the \((n, 1)\) vector \(\boldsymbol{\epsilon}\) as the error of all the samples:

Recall that the column of \(1\)’s is there to account for the intercept \(b_0\). Then,

That’s pretty much the univariate multiple linear regression model. We can now rewrite the optimization problem. Recall that we want to minimize the sum of squared errors. As it turns out, this is equivalent to the dot product of \(\boldsymbol{\epsilon}\) with itself:

Now our minimization now looks like this:

The setup is finally done. Now we substitute \(\boldsymbol{\epsilon} = \mathbf{y} - \mathbf{X} \mathbf{b}\), take the derivative, set to 0, and solve for our coefficients. If you’re worried about taking the derivative of an expression with vectors and matrices, it’s not actually all that bad. Vector derivative rules share many similarities to the univariate derivative rules we learned in basic calculus. First, let’s consider the derivative of \((y-xb)^2\) with respect to \(b\):

Now back to our objective function for multiple linear regression:

We’re almost there! Now we just set the derivative equal to 0 and solve for the values of \(\mathbf{b}\):

And that’s how we get the OLS estimator of the coefficients for multiple linear regression.

Multivariate multiple linear regression

We can further define multivariate multiple linear regression, where there are multiple measured dependent variables and independent variables. In matrix form, we get something like this:

The response variables are now in a matrix \(\mathbf{Y}\), shape \((n, k)\), to account for the \(k\)-dimensional response observed per sample. The coefficient weights are also a matrix \(\mathbf{B}\), shape \((p+1, k)\), to accomodate weights for each of the dimenions of the response variables. \(\mathbf{E}\) is now the error matrix, which gives us the error for all dimensions of all samples. Finally, the design matrix \(\mathbf{X}\) doesn’t change at all.

If we decide to fit the coefficient matrix \(\mathbf{B}\) by minimizing the squared loss, as we did for the univariate multiple regression, we get something that looks familiar:

Interpreting these coefficients might be hard, but estimating them can actually be pretty straightforward!

Finally, let’s see a toy example of this in action. I’ll do a simple linear regression model: \(\mathbf{y} = b_0 + b_1\mathbf{x}\).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
# simulating simple linear regression
import matplotlib.pyplot as plt   # for plotting

# pick random b_0 and b_1 - let's see how we recover these
b = np.array([3.4, 12.35])    # slope, intercept
b = b[:, np.newaxis]

# simulate n datapoints. the x is normally distributed
n = 30
data_mean = 0
data_std = 1
data = np.random.normal(data_mean, data_std, size=(n,1)) # make it 2d array for ease

# remember that since we have b_0 in our vector weights, we add a column of 1s to our data matrix
ones = np.ones((n,1))
x = np.hstack((ones, data)) # x is our design matrix

# add gaussian noise
noise_loc = 0
noise_scale = 5
e = np.random.normal(loc=noise_loc, scale=noise_scale, size=(n,1))

# simulate y values
y = x.dot(b) + e

# get estimates for b, predictions
b_hat = np.linalg.inv(x.T @ x) @ x.T @ y
y_hat = x.dot(b_hat)

# R^2
SS_res = e.T @ e
std = y - y.mean()
SS_tot = std.T @ std
r2 = 1 - (SS_res / SS_tot)

# plot data
fig, ax = plt.subplots()
ax.scatter(data, y) # plotting data
ax.set(xlabel='x', ylabel='y', title='simple linear regression')

# plot true line
b_0, b_1 = b
ax.plot(data, b_0 + (b_1 * data), color='k', label='ground truth')

# plot our fit
ax.plot(data, y_hat, color='r', label='our fit')
ax.legend()

print('true b:\n', b)
print('our estimate:\n', b_hat)
print('r^2:', r2)

1
2
3
4
5
6
7
true b:
 [[ 3.4 ]
 [12.35]]
our estimate:
 [[ 1.78464009]
 [14.5354385 ]]
r^2: [[0.80697382]]