linear_regression_rainfall
In rTorch: R Bindings to 'PyTorch'

Source: https://medium.com/dsnet/linear-regression-with-pytorch-3dde91d60b50

Original title: Linear Regression and Gradient Descent from scratch in PyTorch

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(rTorch)

Select device

torch$manual_seed(0)

device = torch$device('cpu')

Training data

The training data can be represented using 2 matrices (inputs and targets), each with one row per observation, and one column per variable.

# Input (temp, rainfall, humidity)
inputs = np$array(list(list(73, 67, 43),
                   list(91, 88, 64),
                   list(87, 134, 58),
                   list(102, 43, 37),
                   list(69, 96, 70)), dtype='float32')

# Targets (apples, oranges)
targets = np$array(list(list(56, 70), 
                    list(81, 101),
                    list(119, 133),
                    list(22, 37), 
                    list(103, 119)), dtype='float32')

Convert to tensors

Before we build a model, we need to convert inputs and targets to PyTorch tensors.

# Convert inputs and targets to tensors
inputs = torch$from_numpy(inputs)
targets = torch$from_numpy(targets)

print(inputs)
print(targets)

The weights and biases can also be represented as matrices, initialized with random values. The first row of $w$ and the first element of $b$ are used to predict the first target variable, i.e. yield for apples, and, similarly, the second for oranges.

# random numbers for weights and biases. Then convert to double()
torch$set_default_dtype(torch$double)

w = torch$randn(2L, 3L, requires_grad=TRUE)  #$double()
b = torch$randn(2L, requires_grad=TRUE)      #$double()

print(w)
print(b)

Build the model

The model is simply a function that performs a matrix multiplication of the input $x$ and the weights $w$ (transposed), and adds the bias $b$ (replicated for each observation).

model <- function(x) {
  wt <- w$t()
  return(torch$add(torch$mm(x, wt), b))
}

Generate predictions

The matrix obtained by passing the input data to the model is a set of predictions for the target variables.

# Generate predictions
preds = model(inputs)
print(preds)

# Compare with targets
print(targets)

Because we've started with random weights and biases, the model does not a very good job of predicting the target variables.

Loss Function

We can compare the predictions with the actual targets, using the following method:

Calculate the difference between the two matrices (preds and targets).
Square all elements of the difference matrix to remove negative values.
Calculate the average of the elements in the resulting matrix.

The result is a single number, known as the mean squared error (MSE).

# MSE loss
mse = function(t1, t2) {
  diff <- torch$sub(t1, t2)
  mul <- torch$sum(torch$mul(diff, diff))
  return(torch$div(mul, diff$numel()))
}

# Compute loss
loss = mse(preds, targets)
print(loss)
# 46194
# 33060.8070

The resulting number is called the loss, because it indicates how bad the model is at predicting the target variables. Lower the loss, better the model.

Compute Gradients

With PyTorch, we can automatically compute the gradient or derivative of the loss w.r.t. to the weights and biases, because they have requires_grad set to True.

# Compute gradients
loss$backward()

The gradients are stored in the .grad property of the respective tensors.

# Gradients for weights
print(w)
print(w$grad)

# Gradients for bias
print(b)
print(b$grad)

A key insight from calculus is that the gradient indicates the rate of change of the loss, or the slope of the loss function w.r.t. the weights and biases.

If a gradient element is positive:
increasing the element's value slightly will increase the loss.
decreasing the element's value slightly will decrease the loss.
If a gradient element is negative,
increasing the element's value slightly will decrease the loss.
decreasing the element's value slightly will increase the loss.

The increase or decrease is proportional to the value of the gradient.

Finally, we'll reset the gradients to zero before moving forward, because PyTorch accumulates gradients.

# Reset the gradients
w$grad$zero_()
b$grad$zero_()

print(w$grad)
print(b$grad)

Adjust weights and biases using gradient descent

We'll reduce the loss and improve our model using the gradient descent algorithm, which has the following steps:

Generate predictions
Calculate the loss
Compute gradients w.r.t the weights and biases
Adjust the weights by subtracting a small quantity proportional to the gradient
Reset the gradients to zero

# Generate predictions
preds = model(inputs)
print(preds)

# Calculate the loss
loss = mse(preds, targets)
print(loss)

# Compute gradients
loss$backward()

print(w$grad)
print(b$grad)

# Adjust weights and reset gradients
with(torch$no_grad(), {
  print(w); print(b)    # requires_grad attribute remains
  w$data <- torch$sub(w$data, torch$mul(w$grad$data, torch$scalar_tensor(1e-5)))
  b$data <- torch$sub(b$data, torch$mul(b$grad$data, torch$scalar_tensor(1e-5)))

  print(w$grad$data$zero_())
  print(b$grad$data$zero_())
})

print(w)
print(b)

With the new weights and biases, the model should have a lower loss.

# Calculate loss
preds = model(inputs)
loss = mse(preds, targets)
print(loss)

Train for multiple epochs

To reduce the loss further, we repeat the process of adjusting the weights and biases using the gradients multiple times. Each iteration is called an epoch.

# Running all together
# Adjust weights and reset gradients
for (i in 1:100) {
  preds = model(inputs)
  loss = mse(preds, targets)
  loss$backward()
  with(torch$no_grad(), {
    w$data <- torch$sub(w$data, torch$mul(w$grad, torch$scalar_tensor(1e-5)))
    b$data <- torch$sub(b$data, torch$mul(b$grad, torch$scalar_tensor(1e-5)))

    w$grad$zero_()
    b$grad$zero_()
  })
}

# Calculate loss
preds = model(inputs)
loss = mse(preds, targets)
print(loss)

# predictions
preds

# Targets
targets