I've been curious about what goes on under the hood in Linear Regression implementations, particularly those based on matrix factorization. The following is a combination of me exploring matrix factorization and learning PyTorch at the same time. I have a small amount of experience with PyTorch applied to Deep Learning but I would like to understand it better so that I can understand PyTorch implementations of academic papers.
I learned from the second answer on this stackexchange post that there are three ways to implement linear regression with matrix factorization: Cholesky Decomposition, QR Decomposition, and SVD. I chose to explore QR first because it was the technique I was least familiar with.
I began with the following code from rosettacode. It finds Q and R using Householder reflections. Q is orthonormal and R is upper triangular.
import numpy as np
import torch
from sklearn import datasets
from sklearn.linear_model import LinearRegression
import pandas as pd
def qr(A):
m, n = A.shape
Q = np.eye(m)
for i in range(n - (m == n)):
H = np.eye(m)
H[i:, i:] = make_householder(A[i:, i])
Q = np.dot(Q, H)
A = np.dot(H, A)
return Q, A
def make_householder(a):
v = a / (a[0] + np.copysign(np.linalg.norm(a), a[0]))
v[0] = 1
H = np.eye(a.shape[0])
H -= (2 / np.dot(v, v)) * np.dot(v[:, None], v[None, :])
return H
Translating from Numpy from PyTorch is generally pretty easy! Many of the common functions Numpy users are familiar with have analogues in PyTorch.
.size() is the PyTorch equivalent of .shape in Numpy. Note that it is a function called on the Tensor object (oh yeah, arrays are called Tensors in Torch) rather than an attribute of it as in numpy.ndarray.shape.
The implementation for .dot is different between the two libraries. In numpy, .dot will generally do what you want it to: a matrix mutiplication on 2-D arrays and an inner product on 1-D arrays. However, .dot in PyTorch is always treated as a 1-D inner product. It is suggested in this issue thread that the "numpy developers are not happy with the way dot function works," so it may be best to get used to using numpy.matmul and numpy.dot at the appropriate times ala torch.mm and torch.dot.
def qr_torch(A):
m, n = A.size()
Q = torch.eye(m)
for i in range(n - (m == n)):
H = torch.eye(m)
H[i:, i:] = make_householder_torch(A[i:, i])
Q = torch.mm(Q, H)
A = torch.mm(H, A)
return Q, A
def make_householder_torch(a):
v = a / (a[0] + np.copysign(torch.norm(a), a[0]))
v[0] = 1
H = torch.eye(a.size()[0])
H -= (2 / torch.dot(v, v)) * torch.mm(v[:, None], v[None, :])
return H
In order to compare the functions to known good solutions like sklearn, we need to actually produce coefficient estimates, not just Q and R matrices. To do this, I wrote a quick and dirty back-substitution algorithm based on the math from this clear, concise explanation of QR Decomposition applied to linear regession. Note that the z input for the backsub function should be QTy.
def backsub(r, z, n):
betas = np.empty(n)
for i in range(n-1, -1, -1):
betas[i] = (z[i] - sum([betas[j]*r[i,j] for j in range(n-1, i, -1)])) / r[i,i]
return betas
def backsub_torch(r, z, n):
betas = torch.FloatTensor(n, 1)
for i in range(n-1, -1, -1):
betas[i] = (z[i] - sum([betas[j]*r[i,j] for j in range(n-1, i, -1)])) / r[i,i]
return betas
I wanted at least one small dataset and one large dataset, with the idea being that once the PyTorch implementation is CUDA-fied, it should outperform the Numpy implementation on the large dataset because the time it spends loading the data onto the GPU is outweighed by the savings in computation time. For the small dataset, I chose Iris, a commonly used set available in most statistical packages like R and sklearn.
Iris is a dataset of flower measurements. It is composed of 150 observations of 4 features (Sepal Length, Sepal Width, Petal Length and Petal Width) which makes it a good candidate for comparing algorithm performance on a small dataset.
iris = datasets.load_iris()
iris_x = iris.data[:,1:].astype(np.float32)
iris_y = iris.data[:,0].astype(np.float32)
lm = LinearRegression()
lm.fit(iris_x, iris_y)
lm.intercept_, lm.coef_
sklearn takes care of this automatically, but for our qr function, we need to append a column of 1's at the front of our data matrix to represent the intercept.
iris_x_np = np.concatenate((np.ones((iris_x.shape[0],1)), iris_x), axis=1).astype(np.float32)
q, r = qr(iris_x_np)
betas = backsub(r, np.dot(q.T, iris_y), r.shape[1])
betas
x_torch_iris = torch.from_numpy(iris_x_np)
y_torch_iris = torch.from_numpy(iris_y)[:,None]
q_torch, r_torch = qr_torch(x_torch_iris)
betas = backsub_torch(r_torch, torch.mm(q_torch.t(), y_torch_iris), r_torch.size()[1])
betas
Great! Our code is working on Iris. For the large dataset I wanted, I tried the UCI Human Activity Recognition dataset. I chose this because all the data is numeric and I didn't want to spend time one-hot encoding categorical variables.
UCI is a dataset composed of 10299 observations of 561 features, all of which are measurements taken from the accelerometer and gyroscope inside a smartphone. The purpose of the data is to try to classify which movement a person is making based on the instrument readings.
UCIpath = '/home/vwrideout/PyTorchQR/UCIData/UCI HAR Dataset/train/'
dfx = pd.read_csv(UCIpath + 'X_train.txt', header = None, delim_whitespace=True)
dfy = pd.read_csv(UCIpath + 'y_train.txt', header = None)
har_x = dfx.as_matrix().astype(np.float32)
har_y = dfy.as_matrix().astype(np.float32)
har_x.shape
This dataset has 7352 rows and 561 columns. Pretty big!
lm = LinearRegression()
lm.fit(har_x, har_y)
lm.intercept_, lm.coef_
har_x_np = np.concatenate((np.ones((har_x.shape[0],1)), har_x), axis=1).astype(np.float32)
q, r = qr(har_x_np)
betas = backsub(r, np.dot(q.T, har_y), r.shape[1])
betas
Uh oh! What happened here? If you compare them, you can see that the coefficient estimates from our Numpy function start out looking like sklearn's, but by the end of the backsub function (which is the front of the arrays printed here... remember: back-substitution) the coefficients have exploded in size!
np.min(abs(np.diag(r)))
This looks like the culprit. The data is highly collinear, so we wind up with values close to zero on the diagonal of R. During back-substitution we divide by these values and all hell breaks loose. Let's try working with data that is not collinear instead. The Numpy factorization on UCI took hours to run, so I will simulate a dataset of comparable size to make sure the algorithm is sound but then scale it back to get reasonable running times.
I simulated a dataset even larger than UCI Human Activity Recognition in order to determine that our algorithm struggled with multicollinearity and not just the size of the data. It occured to me that floating point addition is imprecise and back-substitution does a boatload of those operations which could explain the behavior we saw above - correct coefficients at the beginning of back-substitution but very incorrect by the end.
sim_x = np.random.normal(size=(10000,600))
coefs = np.array(range(sim_x.shape[1])) * .1
sim_y = np.dot(sim_x, coefs) + np.random.normal(size=(10000,))
lm = LinearRegression()
lm.fit(sim_x, sim_y)
lm.intercept_, lm.coef_
sim_x_np = np.concatenate((np.ones((sim_x.shape[0],1)), sim_x), axis=1).astype(np.float32)
q, r = qr(sim_x_np)
betas = backsub(r, np.dot(q.T, sim_y), r.shape[1])
betas
Here we see that the coefficients produced by our algorithm match sklearn's. It looks like multicollinearity was the problem. Now let's find a simulated dataset size that takes minutes instead of hours to perform QR factorization on with Numpy. I'll try 3000 rows and 300 columns.
When you use the iPython magic %%timeit in a cell, it will run the cell several times and report the best result. We can use this to compare the efficiency of different algorithms.
sim_x = np.random.normal(size=(3000,300))
coefs = np.array(range(sim_x.shape[1])) * .1
sim_y = (np.dot(sim_x, coefs) + np.random.normal(size=(3000,))).astype(np.float32)
sim_x_np = np.concatenate((np.ones((sim_x.shape[0],1)), sim_x), axis=1).astype(np.float32)
%%timeit
q, r = qr(sim_x_np)
betas = backsub(r, np.dot(q.T, sim_y), r.shape[1])
betas
x_torch_sim = torch.from_numpy(sim_x_np)
y_torch_sim = torch.from_numpy(sim_y)[:,None]
%%timeit
q_torch, r_torch = qr_torch(x_torch_sim)
betas = backsub_torch(r_torch, torch.mm(q_torch.t(), y_torch_sim), r_torch.size()[1])
betas
Interesting... PyTorch is faster than Numpy for this operation without even using the GPU.
All it takes to move all the operations onto GPU is to add a .cuda() after each tensor declaration. That's the only difference between these cuda functions and the torch functions from earlier. Note below that I've included the setup code for x_cuda_sim and y_cuda_sim inside the %%timeit cell to account for the time it takes to load the data onto the GPU.
def qr_cuda(A):
m, n = A.size()
Q = torch.eye(m).cuda()
for i in range(n - (m == n)):
H = torch.eye(m).cuda()
H[i:, i:] = make_householder_cuda(A[i:, i])
Q = torch.mm(Q, H)
A = torch.mm(H, A)
return Q, A
def make_householder_cuda(a):
v = a / (a[0] + np.copysign(torch.norm(a), a[0]))
v[0] = 1
H = torch.eye(a.size()[0]).cuda()
H -= (2 / torch.dot(v, v)) * torch.mm(v[:, None], v[None, :])
return H
def backsub_cuda(r, z, n):
betas = torch.FloatTensor(n, 1).cuda()
for i in range(n-1, -1, -1):
betas[i] = (z[i] - sum([betas[j]*r[i,j] for j in range(n-1, i, -1)])) / r[i,i]
return betas
%%timeit
x_cuda_sim = x_torch_sim.cuda()
y_cuda_sim = y_torch_sim.cuda()
q_cuda, r_cuda = qr_cuda(x_cuda_sim)
betas = backsub_cuda(r_cuda, torch.mm(q_cuda.t(), y_cuda_sim), r_cuda.size()[1])
betas
That's an 18x speedup from the Numpy implementation! Now let's check my intuition that Numpy should perform better on a small dataset like Iris.
%%timeit
q, r = qr(iris_x_np)
betas = backsub(r, np.dot(q.T, iris_y), r.shape[1])
betas
%%timeit
q_torch, r_torch = qr_torch(x_torch_iris)
betas = backsub_torch(r_torch, torch.mm(q_torch.t(), y_torch_iris), r_torch.size()[1])
betas
%%timeit
x_cuda_iris = x_torch_iris.cuda()
y_cuda_iris = y_torch_iris.cuda()
q_cuda, r_cuda = qr_cuda(x_cuda_iris)
betas = backsub_cuda(r_cuda, torch.mm(q_cuda.t(), y_cuda_iris), r_cuda.size()[1])
betas
Here we see that the cuda-fied PyTorch algorithm took 1.56ms, which is almost 2.5 times slower than the 640 µs taken by the PyTorch algorithm without GPU acceleration.
It looks like my intuition was correct. GPU accelerated algorithms can more efficiently solve matrix factorization problems on large datasets, but with small datasets like Iris the overhead cost of moving the data onto the GPU outweighs the time savings gained during computation.
I was surprised to see that the algorithms using PyTorch tensors without GPU acceleration were slightly more efficient than the equivalent Numpy algorithms. I think it would be interesting to explore this more in the future. Do tensor declarations take the same amount of time as Numpy array declarations? Does it matter if we change the Numpy algorithm to use numpy.matmul where appropriate?