Math for Machine Learning

Linear Algebra

Kousuke Ariga (kouariga@gmail.com, @kouariga)

[v1] January 17, 2018
[v2] January 18, 2019
[v3] April 10, 2019

The latest version is available here.

This notebook is meant to be a quick refresher of linear algebra and a brief introduction of NumPy (Python package for scientific computing), and it is by no means a through review. I assume that you have a prior experience of learning linear algebra such as taking an introductory course a while ago. The goal is to go over some of the important properties of matrices and showcase some of the NumPy methods through practical examples. We consider linear regression and three different solutions: algebraic, analytic, and geometric. I heavily cite and highly recommend Kolter's review notes on linear algebra [2].

Contents

  1. Introduction
    1.1 Boston house prices dataset
    1.2 Linear regression model
  2. Matrices and vectors
    2.1 Matrix
    2.2 Vector
  3. Basic operations on matrices and vectors
    3.1 Matrix addition
    3.2 Scalar matrix multiplication
    3.3 Matrix vector multiplication
    3.4 Matrix matrix multiplication
  4. Properties of matrices
    4.1 Linear independence
    4.2 Rank
    4.3 Inverse
    4.4 Transpose
    4.5 Span and range
    4.6 Norm
    4.7 Gradient
    4.8 Matrix calculus
    4.9 Dot products
    4.10 Projections
    4.11 Orthogonality
  5. Discussion and conclusions

Appendix
A. Linear algebra visualized
B. Transpose and 1-dimensional arrays in NumPy
C. Outer products in NumPy

Reference

1. Introduction

Why is linear algebra important in machine learning? Machine learning methods often involves a large amount of data, and linear algebra provides a clever way to analyze and manipulate them. To make the argument concrete, let's take a look at a sample dataset.

1.1 Boston house prices dataset

In [1]:
import time

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from IPython.display import Image, IFrame

np.set_printoptions(suppress=True, linewidth=120, precision=2)

We can load boston dataset from sklearn package, which is a very popular and easy to use machine learning package of Python. It implements many kinds of machine learning algorithms and utility functions. The loaded dataset has the following attributes.

In [2]:
boston = load_boston()
print(boston.__dir__())
dict_keys(['data', 'target', 'feature_names', 'DESCR', 'filename'])
In [3]:
print(boston.DESCR)
.. _boston_dataset:

Boston house prices dataset
---------------------------

**Data Set Characteristics:**  

    :Number of Instances: 506 

    :Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

    :Attribute Information (in order):
        - CRIM     per capita crime rate by town
        - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
        - INDUS    proportion of non-retail business acres per town
        - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
        - NOX      nitric oxides concentration (parts per 10 million)
        - RM       average number of rooms per dwelling
        - AGE      proportion of owner-occupied units built prior to 1940
        - DIS      weighted distances to five Boston employment centres
        - RAD      index of accessibility to radial highways
        - TAX      full-value property-tax rate per $10,000
        - PTRATIO  pupil-teacher ratio by town
        - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
        - LSTAT    % lower status of the population
        - MEDV     Median value of owner-occupied homes in $1000's

    :Missing Attribute Values: None

    :Creator: Harrison, D. and Rubinfeld, D.L.

This is a copy of UCI ML housing dataset.
https://archive.ics.uci.edu/ml/machine-learning-databases/housing/


This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.

The Boston house-price data of Harrison, D. and Rubinfeld, D.L. 'Hedonic
prices and the demand for clean air', J. Environ. Economics & Management,
vol.5, 81-102, 1978.   Used in Belsley, Kuh & Welsch, 'Regression diagnostics
...', Wiley, 1980.   N.B. Various transformations are used in the table on
pages 244-261 of the latter.

The Boston house-price data has been used in many machine learning papers that address regression
problems.   
     
.. topic:: References

   - Belsley, Kuh & Welsch, 'Regression diagnostics: Identifying Influential Data and Sources of Collinearity', Wiley, 1980. 244-261.
   - Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.

The data and target values are stored in arrays of type numpy.ndarray. In the data array, each row corresponds to a sample, a Boston suburb or town in this example, and each column corresponds to a feature that is described above. Note that numpy.ndarray is not just multi-dimensional array (or list in Python). It implements many useful numeric methods and indexing feature. Refer to the ndarray document and indexing document for the details. Here, I show the first 10 samples, each of which consists of 13 feature values, and some of their statistics by slicing the data array.

In [4]:
print(boston.feature_names)
print(boston.data[:10])

print('\nmean')
print(boston.data[:10].mean(axis=0))

print('\nvariance')
print(boston.data[:10].var(axis=0))
['CRIM' 'ZN' 'INDUS' 'CHAS' 'NOX' 'RM' 'AGE' 'DIS' 'RAD' 'TAX' 'PTRATIO' 'B' 'LSTAT']
[[  0.01  18.     2.31   0.     0.54   6.58  65.2    4.09   1.   296.    15.3  396.9    4.98]
 [  0.03   0.     7.07   0.     0.47   6.42  78.9    4.97   2.   242.    17.8  396.9    9.14]
 [  0.03   0.     7.07   0.     0.47   7.18  61.1    4.97   2.   242.    17.8  392.83   4.03]
 [  0.03   0.     2.18   0.     0.46   7.    45.8    6.06   3.   222.    18.7  394.63   2.94]
 [  0.07   0.     2.18   0.     0.46   7.15  54.2    6.06   3.   222.    18.7  396.9    5.33]
 [  0.03   0.     2.18   0.     0.46   6.43  58.7    6.06   3.   222.    18.7  394.12   5.21]
 [  0.09  12.5    7.87   0.     0.52   6.01  66.6    5.56   5.   311.    15.2  395.6   12.43]
 [  0.14  12.5    7.87   0.     0.52   6.17  96.1    5.95   5.   311.    15.2  396.9   19.15]
 [  0.21  12.5    7.87   0.     0.52   5.63 100.     6.08   5.   311.    15.2  386.63  29.93]
 [  0.17  12.5    7.87   0.     0.52   6.    85.9    6.59   5.   311.    15.2  386.71  17.1 ]]

mean
[  0.08   6.8    5.45   0.     0.49   6.46  71.25   5.64   3.4  269.    16.78 393.81  11.02]

variance
[   0.     48.66    7.06    0.      0.      0.25  297.44    0.5     2.04 1587.      2.53   14.51   68.29]

The target values are the following. Our task here is to predict the target value, or the "median value of owner-occupied homes in $1000's" in a Boston town, given its feature values such as "per capita crime rate by town" and "average number of rooms per dwelling."

In [5]:
print('MEDV')
print(boston.target[:10])
MEDV
[24.  21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9]

1.2 Linear regression

Linear regression is one of the simplest statistical models. It assumes that the target variable y is explained by a weighted sum of feature values x1,x2,,xn. In an equation,

yMEDV=wCRIMxCRIM+wZNxZN++wLSTATxLSTAT+b

where b is a bias term. Intuitively, wx terms define the relative up/down from the standard target value. This standard value is what the bias term accounts for. You may wonder if the relationship is really that simple.

"Essentially, all models are wrong, but some are useful."

George Box, 1987

Assuming that the linear regression model is valid and we know all the weights and bias, we can estimate a median house price of a Boston town from its feature values. The bad news is that we don't know the weights... The good news is that we have training samples (a set of features and target pair)! We want to find a set of weights such that the equation holds for the training samples. To this end, we can solve systems of equations.

{yMEDV(1)=wCRIMxCRIM(1)+wZNxZN(1)++wLSTATxLSTAT(1)+byMEDV(2)=wCRIMxCRIM(2)+wZNxZN(2)++wLSTATxLSTAT(2)+byMEDV(n)=wCRIMxCRIM(n)+wZNxZN(n)++wLSTATxLSTAT(n)+b

Great, we can solve it, ...or can we (more on this later)? Let's rewrite the equations with a better notation.

[yMEDV(1)yMEDV(2)yMEDV(n)]=[xCRIM(1)xZN(1)xLSTAT(1)1xCRIM(2)xZN(2)xLSTAT(2)1xCRIM(n)xZN(n)xLSTAT(n)1][wCRIMwZNwLSTATb]

More simply,

[y(1)y(2)y(n)]=[x(1)x(2)x(n)]w

or even...

y=Xw

Yes, this is beautiful. This notation is used in linear algebra, and it is a very powerful tool given to us to tackle machine learing problems. The objective here is to find a set of weights w that solves this equation. We call this process to learn from data.

2. Matrices and vectors

2.1 Matrix

A=[012345]AR2×3

A matrix is a rectangular array of numbers. The dimension of matrix is number of rows by number of columns. Aij is the i,j entry of A, which is in the i th row and the j th column.

In [6]:
A = np.array(np.arange(0, 6)).reshape((2, 3))

print(A)
print(A.shape)

for i in range(A.shape[0]):
    for j in range(A.shape[1]):
        print("{},{} entry: {}".format(i, j, A[i, j]))
[[0 1 2]
 [3 4 5]]
(2, 3)
0,0 entry: 0
0,1 entry: 1
0,2 entry: 2
1,0 entry: 3
1,1 entry: 4
1,2 entry: 5

2.2 Vector

y=[0246]yR4

A vector is a n×1 matrix. Here y is said to be a 4-dimensional vector because it has 4 elements in it. yi denotes the i th element of y.

In [7]:
y = np.array(2*np.arange(0, 4))
print(y)
print(y.shape)

for i in range(y.shape[0]):
    print("{} element: {}".format(i, y[i]))
[0 2 4 6]
(4,)
0 element: 0
1 element: 2
2 element: 4
3 element: 6

3. Basic operations on matrices and vectors

3. 1 Matrix addition

[102531]+[40.52501]=[50.541032]

The shapes have to be the same.

In [8]:
A = np.array([[1, 0],
              [2, 5],
              [3, 1]])
B = np.array([[4, 0.5],
              [2, 5],
              [0, 1]])
assert A.shape == B.shape
print(A + B)
[[ 5.   0.5]
 [ 4.  10. ]
 [ 3.   2. ]]

3.2 Scalar matrix multiplication

3[102531]=[3061593]
In [9]:
A = np.array([[1, 0],
              [2, 5],
              [3, 1]])
print(3*A)
[[ 3  0]
 [ 6 15]
 [ 9  3]]

3.3 Matrix vector multiplication

Ax=y

A:m×n matrix (m rows, n columns)

x:n×1 matrix (n-dimensional vector)

y:m×1 matrix (m-dimensional vector)

To get yi, multiply A's ith row with vector x element-wise, then add them up.

[121503041200][1321]=?

Hint: R3×4×R4×1=R3×1

In [10]:
A = np.array([[1, 2, 1, 5],
              [0, 3, 0, 4],
              [-1, -2, 0, 0]])
x = np.array([[1],
              [3],
              [2],
              [1]])

assert x.shape[0] == A.shape[1]
y = np.dot(A, x)
y = A.dot(x)  # another way to calculate the dot product

print(y)
[[14]
 [13]
 [-7]]

3.4 Matrix matrix multiplication

AB=C

A:=l×m matrix (l rows, m columns matrix)

B:=m×n matrix (m rows, n columns matrix)

C:=l×n matrix (l rows, n columns matrix)

[145326][18745623]=?

Hint: R3×2×R2×4=R3×4

Note AB and BA are not the same, i.e. matrix multiplication is NOT commutative. Actually, the latter is not necessarily defined. Check dimension.

In [11]:
A = np.array([[1, 4],
              [5, 3],
              [2, 6]])
B = np.array([[1, 8, 7, 4],
              [5, 6, 2, 3]])
print(A)
print(B)
[[1 4]
 [5 3]
 [2 6]]
[[1 8 7 4]
 [5 6 2 3]]
In [12]:
print(A.dot(B))
[[21 32 15 16]
 [20 58 41 29]
 [32 52 26 26]]
In [13]:
try:
    print(B.dot(A))
except ValueError as e:
    print(e)
shapes (2,4) and (3,2) not aligned: 4 (dim 1) != 3 (dim 0)

4. Properties of matrices

Is linear algebra all about saving papers? Definitely NO! Do you remember terminologies such as linear independence, rank, span, etc that you learned in the linear algebra course? Did you get the ideas? Being able to calculate them is important, but understanding the concept is more (at least as) important for the purpose of this course (we can use computers for calculation after all). Let's review the properties of matrices while solving the linear regression of Boston house price. The goal is to solve the following equation.

y=Xw

where XRm×n and m>n. m is greater than n because there are more samples than the number of features (remember rows are samples). In other words, X is a vertically long matrix.

4.1 Linear independence

Here, let's assume that all the features (columns of X) are linearly independent.

A set of vectors {x1,x2,,xn}Rm is said to be (linearly) independent if no vector can be represented as a linear combination of the remaining vectors. [2]

Otherwise, it is linearly dependent. For example, if we have temperature in Fahrenheit and in Celsius as two different features, the latter is represented in terms of the first as

F=95C+32.

Such features are linearly dependent. For another example, if we have categorical features like gender, we could have two columns one for male and the other for female. For male samples we have ones in the male column and zeros in the female column, and the opposite for female samples. Did you notice that we have a linear dependence here because these features can be represented in the form

F=M+1.

4.2 Rank

For a matrix ARm×n where m>n, if its columns are linearly independent, it is said to be full rank. Formally,

The column rank of a matrix ARm×n is the size of the largest subset of columns of A that constitute a linearly independent set. With some abuse of terminology, this is often referred to simply as the number of linearly independent columns of A. In the same way, the row rank is the largest number of rows of A that constitute a linearly independent set.

For any matrix ARm×n, it turns out that the column rank of A is equal to the row rank of A (though we will not prove this), and so both quantities are referred to collectively as the rank of A, denoted as rank(A).

For ARm×n, rank(A)min(m,n). If rank(A)=min(m,n), then A is said to be full rank. [2]

Therefore, the first statement holds.

4.3 Inverse

The inverse of a square matrix ARn×n is denoted A1, and is the unique matrix such that

A1A=I=AA1.

Note that not all matrices have inverses. Non-square matrices, for example, do not have inverses by definition. However, for some square matrices A, it may still be the case that A1 may not exist. In particular, we say that A is invertible or non-singular if A1 exists and non-invertible or singular otherwise. In order for a square matrix A to have an inverse A1, then A must be full rank. [2]

So, again, let's assume that the columns of X are linearly independent, i.e. X is full rank. Here's our first attempt to solve the equation for w.

The first approach (wrong)

y=Xww=X1y

Can we do this?

Remember that our X is a vertically long matrix i.e. non-square, and cannot be inverted by definition.

Therefore, we can't do w=X1y.

4.4 Transpose

By convention, an n-dimensional vector is often thought of as a matrix with n rows and 1 column, known as a column vector. If we want to explicitly represent a row vector — a matrix with 1 row and n columns — we typically write xT. (2)

y=[0123]yT=[0123]

The transpose can be generalized to matrices.

The transpose of a matrix results from "flipping" the rows and columns. Given a matrix ARm×n, its transpose, written ATRn×m, is the n×m matrix whose entries are given by

(AT)ij=Aji.

(2)

A=[123456]AT=[135246]

What is the dimension of ATA?

In [14]:
A = np.array([[1, 2],
              [3, 4],
              [5, 6]])
print(A)
print(np.dot(A.T, A))
[[1 2]
 [3 4]
 [5 6]]
[[35 44]
 [44 56]]

ATA is always a square matrix (Rn×mRm×n=Rn×n). If A is full rank, ATA is also invertible.

The second approach (algebraic)

y=XwXTy=XTXww=(XTX)1XTy

Note that the second line multiplies both sides by the transpose of X from the left. Then, we can solve for w because XTX is invertible.

This is a valid algebraic approach.

Why don't we get the intuition behind the solution. Consider the linear system

y=Xw

where

y=[122]X=[111213]w=[w1w2].

This equation is saying that y is a linear combination of column vectors of X.

[122]=w1[111]+w2[123]

However, there are no such weights.

4.5 Span

With linear algebra's terminology, y doesn't lie in the column space of X, or the space that column vectors of X spans. Formally,

The span of a set of vectors {x1,x2,,xn} is the set of all vectors that can be expressed as a linear combination of {x1,x2,,xn}. That is,

span({x1,,xn})={v:v=i=1nαixi,αiR}.

(2)

Especially, when x's are the columns of a matrix X, their span is said to be the range or the column space of X and denoted R(X).

Back to the equation, although the target vector y is 3-dimensional, there are only two column vectors that span the space, i.e. the range of X is just a 2-dimensional plane. Therefore, there certainly exists 3-dimensional vectors that don't lie on this plane, like the y above. Visually, it looks something like below.

In [15]:
Image('../images/linear-algebra/4.5.png', height=300, width=400)
Out[15]:

But we want to represent y in terms of xi's. The best we can do is to find a vector that lies in the range of X, but is also as close to y as possible.

4.6 Norm

This objective can be formulated by using norm by saying to find w that minimizes ||yXw||2.

A norm of a vector x is informally a measure of the “length” of the vector. For example, we have the commonly-used Euclidean or 2 norm,

x2=i=1nxi2.

Note that x22=xTx. (2)

If you take the norm of difference of vectors, it is a measure of distance between them. There are several types of norms, but another popular one is 1 norm. Given a vector xRn,

x1=i=0n|xi|

Let's use 2 norm as a measure of distance for now. For convinience, we can minimize the square of 2 norm without loss of generality. To find weights that minimizes yXw22, we can take its derivative with respect to w and set to zero. Easy, right?

4.7 Gradient

To this end, the notion of gradient, which is a natural extension of partial derivatives to a vector setting, comes in handy.

Suppose that f:Rm×nR is a function that takes as input a matrix A of size m×n and returns a real value. Then the gradient of f (with respect to ARm×n) is the matrix of partial derivatives, defined as:

Af(A)Rm×n=[f(A)A11f(A)A12f(A)A1nf(A)A21f(A)A22f(A)A2nf(A)Am1f(A)Am2f(A)Amn]

i.e., an m×n matrix with