Matrix Magic: Building Machine Learning Models from the Ground Up#

Scalars, Vectors, Matrices, and Tensors#

• Scalars are single numbers, denoted typically by lowercase letters (e.g., x).
• Vectors are one-dimensional arrays of scalars, typically denoted by bold lowercase letters (e.g., x).
• Matrices are two-dimensional arrays of scalars, denoted by bold uppercase letters (e.g., X).
• Tensors generalize the concept of arrays to more than two dimensions.

Think of a 3D array of data (height × width × depth) as a third-order tensor. For most machine learning tasks, we primarily operate with vectors and matrices, occasionally moving into the world of higher-order tensors (such as in deep learning).

Matrix Operations#

Matrix operations are the building blocks of machine learning algorithms. Some key operations are:

Matrix Inversion and Rank#

• Matrix Inversion is crucial in deriving closed-form solutions (like the normal equation in linear regression). However, not all matrices are invertible.
• Rank of a matrix determines the maximum number of linearly independent row or column vectors. A square matrix is invertible if its rank is equal to the dimension of the matrix (i.e., full rank).

Eigenvalues and Eigenvectors#

For a square matrix A, an eigenvector v and corresponding eigenvalue λ satisfy:

A v = λ v

Eigen-decomposition and singular value decomposition (SVD) become critical in PCA (Principal Component Analysis) and other dimensionality-reduction techniques.

Problem Setup#

Suppose you have a dataset of m training examples, each with n features. These can be arranged into a matrix X (an m×n matrix). Each row of X corresponds to one sample, and you also have a vector y (m×1) containing the target values. The goal is to find a parameter vector θ (n×1) that predicts y given X using the linear model:

ŷ = Xθ

Batch Gradient Descent#

For linear regression, the cost function (often called the mean squared error) is:

J(θ) = (1 / (2m)) �?(h_θ(x^(i)) - y^(i))²

Where h_θ(x^(i)) is the prediction for the i-th example.

We seek to find θ that minimizes J(θ). Using gradient descent, we repeatedly update:

θ := θ - α (1 / m) Xᵀ(Xθ - y)

Here, α is the learning rate.

Vectorization Approach#

Vectorization is the practice of using matrix operations instead of explicit loops. This approach allows for more efficient computation and leverages optimized libraries (often backed by GPU acceleration) to handle large datasets swiftly. For linear regression, it’s typical to:

1. Compute predictions in one shot: Xθ.
2. Compute error vector: Xθ - y.
3. Multiply by Xᵀ and update θ.

Code Example: Implementing Linear Regression in Python#

Below is a simplified implementation of linear regression using NumPy.

1
import numpy as np
2

3
# Generate a simple dataset
4
np.random.seed(42)
5
m = 100  # number of samples
6
X = 2 * np.random.rand(m, 1)
7
y = 4 + 3 * X + np.random.randn(m, 1)  # y = 4 + 3x + noise
8

9
# Add a bias term by appending x0=1 to each instance
10
X_b = np.c_[np.ones((m, 1)), X]  # shape (100, 2)
11

12
# Hyperparameters
13
learning_rate = 0.1
14
n_iterations = 1000
15

16
# Initialization of theta
17
theta = np.random.randn(2, 1)
18

19
for iteration in range(n_iterations):
20
    gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
21
    theta = theta - learning_rate * gradients
22

23
print("Learned theta:")
24
print(theta)

Sigmoid Function and Probabilities#

Classification problems often require predicting categorical labels. Logistic regression, despite regression in its name, is used for binary classification (0 or 1). The key is the sigmoid function:

σ(z) = 1 / (1 + e⁻ᶻ)

Instead of directly predicting ŷ = Xθ, we model:

p = σ(Xθ)

where p is the probability of belonging to the positive class (label = 1).

Loss Function: Cross-Entropy#

For logistic regression, we use the cross-entropy loss (or log loss):

J(θ) = - (1/m) �?[y^(i) log(h_θ(x^(i))) + (1 - y^(i)) log(1 - h_θ(x^(i)))]

Gradient Descent for Logistic Regression#

Update rule using gradient descent:

θ := θ - α (1 / m) Xᵀ(σ(Xθ) - y)

Once again, vectorization is key to efficient computation.

Code Example: Logistic Regression Implementation#

1
import numpy as np
2

3
def sigmoid(z):
4
    return 1 / (1 + np.exp(-z))
5

6
# Synthetic data for binary classification
7
np.random.seed(42)
8
m, n = 100, 2
9
X = np.random.randn(m, n)
10
true_theta = np.array([[2], [-1.5]])
11
prob = sigmoid(X.dot(true_theta))
12
y = (prob > 0.5).astype(int)  # Convert probabilities to 0 or 1
13

14
# Add bias term
15
X_b = np.c_[np.ones((m, 1)), X]
16

17
# Hyperparameters
18
learning_rate = 0.1
19
n_iterations = 1000
20
theta = np.zeros((n+1, 1))
21

22
for i in range(n_iterations):
23
    z = X_b.dot(theta)
24
    predictions = sigmoid(z)
25
    gradients = 1/m * X_b.T.dot(predictions - y)
26
    theta -= learning_rate * gradients
27

28
print("Trained parameters:")
29
print(theta)

Ridge (L2) and Lasso (L1)#

Overfitting happens when a model fits the training data too closely and fails to generalize. Regularization addresses this by adding a penalty term to the cost function:

• Ridge (L2) adds (λ/2m)‖θ‖�?to the cost.
• Lasso (L1) adds (λ/m)‖θ‖₁ to the cost.

Here, ‖θ‖�?is the L2-norm squared, and ‖θ‖₁ is the sum of absolute values of θ’s components. Regularization encourages the parameters toward lower magnitudes (Ridge) or sparsity (Lasso).

Matrices and Overfitting#

Regularization can be seen through the lens of matrix operations. By shrinking the parameter vector, we effectively reduce the variance of predictions. In matrix form, a regularized gradient descent update includes an additional term λθ (for L2) or sign-based term for L1.

Forward Propagation#

Neural networks use multiple layers of transformations. Each layer typically carries out:

Z = W · A + b

followed by an activation function. Z is the weighted input, W is a weight matrix, A is the activation from the previous layer, and b is the bias vector. Matrix multiplication is the driving force here, making them extremely well-suited for GPU acceleration.

Backpropagation#

During backpropagation, partial derivatives of the loss function with respect to each weight and bias are computed. This process heavily relies on matrix calculus:

1. Compute forward pass.
2. Calculate loss.
3. Efficiently backpropagate through each layer to compute gradients.
4. Update weights.

From Vectors to Tensors#

As models become deeper, especially in convolutional neural networks (CNNs), the parameters and intermediate data are often 3D or 4D tensors (e.g., height × width × channels × batch size). The same principles of matrix multiplication still apply; we simply do generalized forms of these operations.

Code Example: Simple Neural Network in NumPy#

Below is a minimal demonstration of a single hidden-layer neural network to classify 2D points.

1
import numpy as np
2

3
# Generate synthetic data
4
np.random.seed(42)
5
m = 200
6
X = np.random.randn(m, 2)
7
y = (X[:, 0] * X[:, 1] > 0).astype(int).reshape(m, 1)
8

9
# Helper functions
10
def sigmoid(z):
11
    return 1 / (1 + np.exp(-z))
12

13
def sigmoid_derivative(a):
14
    return a * (1 - a)
15

16
# Network architecture
17
input_dim = 2
18
hidden_dim = 3
19
output_dim = 1
20

21
# Initialize weights
22
W1 = np.random.randn(input_dim, hidden_dim)
23
b1 = np.zeros((1, hidden_dim))
24
W2 = np.random.randn(hidden_dim, output_dim)
25
b2 = np.zeros((1, output_dim))
26

27
learning_rate = 0.1
28
n_iterations = 10000
29

30
for i in range(n_iterations):
31
    # Forward pass
32
    Z1 = X.dot(W1) + b1
33
    A1 = sigmoid(Z1)
34
    Z2 = A1.dot(W2) + b2
35
    A2 = sigmoid(Z2)
36

37
    # Compute loss (binary cross-entropy)
38
    loss = -np.mean(y * np.log(A2) + (1 - y) * np.log(1 - A2))
39

40
    # Backpropagation
41
    dA2 = A2 - y
42
    dZ2 = dA2 * sigmoid_derivative(A2)
43
    dW2 = A1.T.dot(dZ2) / m
44
    db2 = np.sum(dZ2, axis=0, keepdims=True) / m
45

46
    dA1 = dZ2.dot(W2.T)
47
    dZ1 = dA1 * sigmoid_derivative(A1)
48
    dW1 = X.T.dot(dZ1) / m
49
    db1 = np.sum(dZ1, axis=0, keepdims=True) / m
50

51
    # Parameter updates
52
    W2 -= learning_rate * dW2
53
    b2 -= learning_rate * db2
54
    W1 -= learning_rate * dW1
55
    b1 -= learning_rate * db1
56

57
    # Optional: print loss occasionally
58
    if i % 2000 == 0:
59
        print(f"Iteration {i}, loss: {loss:.4f}")
60

61
# Output final loss
62
print(f"Final loss: {loss:.4f}")

Convolutional Neural Networks (CNNs)#

CNNs are specialized for image data (or other spatial data). They use convolutional filters to detect features in an image:

Z = W * A + b

where * denotes convolution rather than matrix multiplication. However, under the hood, convolution can still be expressed using matrix operations via a process called “im2col,” which reorganizes patches of images into columns.

Recurrent Neural Networks (RNNs)#

RNNs process sequential data (like time series or text) by maintaining a hidden state vector that is updated at each time step. The transformations are matrix-based, and the “state” flows over multiple time steps. Variants include LSTM and GRU units, which manage the vanishing or exploding gradient problems through gating mechanisms.

Matrix Factorization and Recommendation Systems#

Recommendation systems often rely on matrix factorization to decompose user-item ratings. If R is a matrix where R(i, j) is the rating of user i for item j, we seek to factorize it into:

R �?P****Qᵀ

Here, P and Q are lower-dimensional matrices capturing latent user and item features. Gradient descent or alternating least squares can solve for P and Q, enabling predictions for missing ratings.

Data Preprocessing#

• Normalization/Standardization: Scale features to have zero mean and unit variance (or to a [0, 1] range).
• One-Hot Encoding for categorical features.
• Dimensionality Reduction (e.g., PCA) to reduce noise and speed up training.

Model Evaluation Metrics#

Beyond accuracy, you may need more nuanced metrics:

• Precision and Recall for binary classification.
• F1-score: harmonic mean of precision and recall.
• AUC-ROC: area under the ROC curve for binary classifiers.
• Mean Absolute Error (MAE) or Mean Squared Error (MSE) for regression tasks.

Scaling to Real-World Datasets#

• Stochastic Gradient Descent (SGD): Update parameters on each sample (or mini-batch) instead of the entire dataset.
• Vectorized Libraries: NumPy, TensorFlow, PyTorch, and others.
• Distributed Computing: Use specialized libraries for large datasets that don’t fit into memory.