Scikit-Learn & PyTorch: ML Algorithms & Models

Derivation of the Closed-Form Solution

Start with the SSE loss function:

Expand the expression:

Compute the gradient with respect to \( \theta \):

Set the gradient to zero to find the minimum:

Assuming \( X^T X \) is invertible, solve for \( \theta \):

Note: \( X^T X \) must be invertible (i.e., \( X \) must have full column rank). If not, regularization (e.g., Ridge Regression) is needed.

Derivation of the Gradient Descent Update Rule

The gradient of \( J(\theta) \) is:

In gradient descent, we update \( \theta \) in the direction opposite to the gradient (since we want to minimize \( J(\theta) \)):

This update is repeated until convergence.

When to Use Closed-Form Solution

• Small Datasets: The closed-form solution is computationally efficient for small datasets (e.g., \( n < 10,000 \) and \( d < 10,000 \)), as it computes the exact solution in one step.
• No Hyperparameters: Unlike gradient descent, the closed-form solution does not require tuning a learning rate or number of iterations.
• Full Rank \( X \): Use when \( X^T X \) is invertible (no multicollinearity).

When to Use Gradient Descent

• Large Datasets: Gradient descent (especially stochastic or mini-batch variants) is scalable to very large datasets where the closed-form solution is computationally infeasible.
• Online Learning: Gradient descent can update the model incrementally as new data arrives, making it suitable for streaming data.
• Non-Invertible \( X^T X \): Use when \( X \) is not full rank (e.g., in high-dimensional data where \( d > n \)).
• Non-Linear Models: Gradient descent is the foundation for training more complex models (e.g., neural networks) where closed-form solutions do not exist.

Implementation in PyTorch and Scikit-Learn

Scikit-Learn (Closed-Form):

from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X, y)
theta = model.coef_  # Parameters (excluding intercept)
intercept = model.intercept_  # Intercept term

Scikit-Learn (Gradient Descent):

from sklearn.linear_model import SGDRegressor
model = SGDRegressor(learning_rate='constant', eta0=0.01)
model.fit(X, y)
theta = model.coef_  # Parameters (excluding intercept)
intercept = model.intercept_  # Intercept term

PyTorch (Gradient Descent):

import torch
import torch.optim as optim

# Define model and loss
X_tensor = torch.tensor(X, dtype=torch.float32)
y_tensor = torch.tensor(y, dtype=torch.float32).reshape(-1, 1)
theta = torch.randn(X.shape[1], 1, requires_grad=True)
loss_fn = torch.nn.MSELoss()

# Gradient descent
optimizer = optim.SGD([theta], lr=0.01)
for epoch in range(1000):
    y_pred = X_tensor @ theta
    loss = loss_fn(y_pred, y_tensor)
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

Diagnosing Gradient Descent Issues

• Divergence: If the loss increases, reduce \( \alpha \) or check for feature scaling issues.
• Slow Convergence: If the loss decreases very slowly, increase \( \alpha \) (but not too much) or use adaptive methods like Adam.
• Oscillations: If the loss oscillates, reduce \( \alpha \) or use a learning rate schedule.
• Plateaus: If the loss plateaus, the model may be stuck in a saddle point. Try increasing \( \alpha \) or using momentum.

Derivation of Ridge Solution:

1. Start with the objective: \[ J(\beta) = \|y - X\beta\|_2^2 + \lambda \|\beta\|_2^2 \]
2. Expand the RSS: \[ J(\beta) = (y - X\beta)^T(y - X\beta) + \lambda \beta^T \beta \]
3. Take the gradient with respect to \( \beta \) and set to zero: \[ \nabla J(\beta) = -2X^T(y - X\beta) + 2\lambda \beta = 0 \]
4. Rearrange: \[ X^T y = X^T X \beta + \lambda \beta = (X^T X + \lambda I)\beta \]
5. Solve for \( \beta \): \[ \beta = (X^T X + \lambda I)^{-1} X^T y \]

Coordinate Descent for Lasso:

For each coefficient \( \beta_j \), the update rule (with other coefficients fixed) is:

• \( \tilde{y}_i^{(j)} = \sum_{k \neq j} x_{ik} \beta_k \) is the partial residual,
• \( S(z, \lambda) = \text{sign}(z)(|z| - \lambda)_+ \) is the soft-thresholding operator.

The soft-thresholding operator drives small coefficients to zero, enabling feature selection.

When to Use Ridge vs. Lasso:

• Ridge Regression:
  • When you have many features with small/medium effects.
  • When features are highly correlated (e.g., genomics, finance).
  • When you want to retain all features but reduce their impact.
• Lasso Regression:
  • When you suspect only a few features are important (sparse models).
  • For feature selection in high-dimensional data (e.g., text mining, bioinformatics).
  • When interpretability is crucial (smaller models).
• Elastic Net:
  • When you have many correlated features (e.g., gene expression data).
  • When \( p \gg n \) and you want to select groups of correlated features.

Scikit-Learn:

from sklearn.linear_model import Ridge, Lasso, ElasticNet
from sklearn.preprocessing import StandardScaler

# Standardize features (critical for regularization)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Ridge Regression
ridge = Ridge(alpha=1.0)
ridge.fit(X_train_scaled, y_train)
y_pred_ridge = ridge.predict(X_test_scaled)

# Lasso Regression
lasso = Lasso(alpha=0.1)
lasso.fit(X_train_scaled, y_train)
y_pred_lasso = lasso.predict(X_test_scaled)

# Elastic Net
elastic = ElasticNet(alpha=0.1, l1_ratio=0.5)
elastic.fit(X_train_scaled, y_train)
y_pred_elastic = elastic.predict(X_test_scaled)

PyTorch (Custom Implementation):

Below is a PyTorch implementation of Ridge Regression using gradient descent:

import torch
import torch.nn as nn
import torch.optim as optim

class RidgeRegression(nn.Module):
    def __init__(self, input_dim):
        super(RidgeRegression, self).__init__()
        self.linear = nn.Linear(input_dim, 1)

    def forward(self, x):
        return self.linear(x)

# Data
X_train = torch.randn(100, 10)  # 100 samples, 10 features
y_train = torch.randn(100, 1)

# Model, loss, optimizer
model = RidgeRegression(input_dim=10)
criterion = nn.MSELoss()
optimizer = optim.SGD(model.parameters(), lr=0.01, weight_decay=1.0)  # weight_decay = lambda

# Training loop
for epoch in range(1000):
    optimizer.zero_grad()
    outputs = model(X_train)
    loss = criterion(outputs, y_train)
    loss.backward()
    optimizer.step()

Note: In PyTorch, weight_decay in optimizers (e.g., SGD, Adam) corresponds to L2 regularization (Ridge). For Lasso, you would need to implement a custom loss function with an L1 penalty.

Common Questions and Answers:

1. Q: What is the difference between Ridge and Lasso regression?

   A: Ridge uses an L2 penalty (\( \|\beta\|_2^2 \)) and shrinks coefficients toward zero without setting them exactly to zero. Lasso uses an L1 penalty (\( \|\beta\|_1 \)) and can set some coefficients to exactly zero, performing feature selection. Ridge is better for handling multicollinearity, while Lasso is better for sparse models.

2. Q: Why is feature scaling important for regularized regression?

   A: Regularization penalties are sensitive to the scale of features. Without scaling, features with larger magnitudes dominate the penalty term, leading to biased coefficient estimates. Standardization (mean=0, variance=1) ensures all features contribute equally to the penalty.

3. Q: How do you choose the regularization parameter \( \lambda \)?

   A: Use cross-validation (e.g., k-fold CV) to evaluate model performance for different values of \( \lambda \). Select the \( \lambda \) that minimizes the validation error. In scikit-learn, this can be done using RidgeCV or LassoCV.

4. Q: What is the Bayesian interpretation of Ridge and Lasso?

   A: Ridge corresponds to placing a Gaussian prior on the coefficients, while Lasso corresponds to placing a Laplace prior. The regularization parameter \( \lambda \) is related to the variance of the prior distribution. The MAP estimate under these priors yields the Ridge and Lasso solutions.

5. Q: When would you use Elastic Net over Ridge or Lasso?

   A: Use Elastic Net when you have many correlated features and want to select groups of features together. It combines the benefits of Ridge (handling multicollinearity) and Lasso (feature selection) and is particularly useful in high-dimensional settings where \( p \gg n \).

6. Q: Why doesn't Lasso have a closed-form solution?

   A: The L1 penalty (\( \|\beta\|_1 \)) is not differentiable at \( \beta_j = 0 \), which prevents the derivation of a closed-form solution. Instead, iterative methods like coordinate descent or proximal gradient descent are used to solve the optimization problem.

Assume we have \( N \) independent observations \( \{(\mathbf{x}_i, y_i)\}_{i=1}^N \), where \( y_i \in \{0, 1\} \). The likelihood of observing the data given the parameters \( \mathbf{w} \) and \( b \) is: \[ L(\mathbf{w}, b) = \prod_{i=1}^N P(y_i \mid \mathbf{x}_i; \mathbf{w}, b) = \prod_{i=1}^N \sigma(\mathbf{w}^T \mathbf{x}_i + b)^{y_i} \cdot \left(1 - \sigma(\mathbf{w}^T \mathbf{x}_i + b)\right)^{1 - y_i} \]

Taking the natural logarithm of the likelihood simplifies the product into a sum: \[ \ell(\mathbf{w}, b) = \log L(\mathbf{w}, b) = \sum_{i=1}^N \left[ y_i \log(\sigma(\mathbf{w}^T \mathbf{x}_i + b)) + (1 - y_i) \log(1 - \sigma(\mathbf{w}^T \mathbf{x}_i + b)) \right] \]

To find the parameters \( \mathbf{w} \) and \( b \) that maximize the log-likelihood, we take the gradient of \( \ell(\mathbf{w}, b) \) with respect to \( \mathbf{w} \) and \( b \) and set it to zero. However, the resulting equations are nonlinear and do not have a closed-form solution. Instead, we use optimization techniques like gradient descent.

The gradient of the log-likelihood with respect to \( \mathbf{w} \) is: \[ \nabla_{\mathbf{w}} \ell(\mathbf{w}, b) = \sum_{i=1}^N \left( y_i - \sigma(\mathbf{w}^T \mathbf{x}_i + b) \right) \mathbf{x}_i \] Similarly, the gradient with respect to \( b \) is: \[ \nabla_{b} \ell(\mathbf{w}, b) = \sum_{i=1}^N \left( y_i - \sigma(\mathbf{w}^T \mathbf{x}_i + b) \right) \]

Using gradient descent, the parameters are updated iteratively: \[ \mathbf{w} := \mathbf{w} + \alpha \sum_{i=1}^N \left( y_i - \sigma(\mathbf{w}^T \mathbf{x}_i + b) \right) \mathbf{x}_i \] \[ b := b + \alpha \sum_{i=1}^N \left( y_i - \sigma(\mathbf{w}^T \mathbf{x}_i + b) \right) \] where \( \alpha \) is the learning rate.

import torch
import torch.nn as nn
import torch.optim as optim

# Define the model
class LogisticRegression(nn.Module):
    def __init__(self, input_dim):
        super(LogisticRegression, self).__init__()
        self.linear = nn.Linear(input_dim, 1)

    def forward(self, x):
        return torch.sigmoid(self.linear(x))

# Initialize model, loss, and optimizer
model = LogisticRegression(input_dim=2)
criterion = nn.BCELoss()  # Binary Cross-Entropy Loss
optimizer = optim.SGD(model.parameters(), lr=0.01)

# Example data (2 features, binary labels)
X = torch.tensor([[1.0, 2.0], [2.0, 3.0], [3.0, 4.0]], dtype=torch.float32)
y = torch.tensor([[0.0], [1.0], [1.0]], dtype=torch.float32)

# Training loop
for epoch in range(1000):
    # Forward pass
    outputs = model(X)
    loss = criterion(outputs, y)

    # Backward pass and optimization
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

    if (epoch + 1) % 100 == 0:
        print(f'Epoch [{epoch+1}/1000], Loss: {loss.item():.4f}')
    

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score

# Example data
X = [[1, 2], [2, 3], [3, 4], [4, 5]]
y = [0, 1, 1, 0]

# Split data into training and test sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.25, random_state=42)

# Initialize and train the model
model = LogisticRegression(solver='liblinear')
model.fit(X_train, y_train)

# Predict and evaluate
y_pred = model.predict(X_test)
print(f'Accuracy: {accuracy_score(y_test, y_pred):.2f}')

# Coefficients and intercept
print(f'Coefficients: {model.coef_}')
print(f'Intercept: {model.intercept_}')
    

Derivation of the Bias-Variance Decomposition

Let \( y = f(x) + \epsilon \), where \( \epsilon \) is the irreducible error with \( \mathbb{E}[\epsilon] = 0 \) and \( \text{Var}(\epsilon) = \sigma^2 \). The expected prediction error is:

Substitute \( y = f(x) + \epsilon \):

Expand the square:

Since \( \mathbb{E}[\epsilon] = 0 \), the cross term vanishes:

Note that \( \mathbb{E}[\epsilon^2] = \text{Var}(\epsilon) = \sigma^2 \). Now, focus on the first term:

Let \( \text{Bias}(\hat{f}(x)) = \mathbb{E}[\hat{f}(x)] - f(x) \). Then:

Expand the square:

The last term is zero because \( \mathbb{E}[\hat{f}(x) - \mathbb{E}[\hat{f}(x)]] = 0 \). Thus:

Putting it all together:

Practical Example: Polynomial Regression

Consider fitting a polynomial of degree \( d \) to data generated from a true function \( f(x) \).

• For \( d = 1 \) (linear model): High bias, low variance. The model is too simple and underfits the data.
• For \( d = 3 \): Bias and variance are balanced, leading to good generalization.
• For \( d = 10 \): Low bias, high variance. The model fits the training data very well but overfits, leading to poor performance on unseen data.

The optimal degree \( d \) can be selected using cross-validation to minimize the expected prediction error.

Visualizing the Bias-Variance Tradeoff

Consider the following plot of error vs. model complexity:

• The training error decreases monotonically as model complexity increases.
• The test error follows a U-shaped curve: it decreases initially (as bias decreases) but then increases (as variance dominates).
• The optimal model complexity minimizes the test error.

Example: Consider a binary classification problem with a node containing 4 samples of class A and 6 samples of class B.

Proportions: \( p_A = 0.4 \), \( p_B = 0.6 \)

Example: Using the same dataset as above (4 samples of class A and 6 samples of class B).

Example: Suppose we have a dataset \( D \) with 10 samples (4 class A, 6 class B) and a feature \( A \) with two possible values: \( v_1 \) and \( v_2 \). After splitting:

• Subset \( D_{v_1} \): 3 samples (2 class A, 1 class B).
• Subset \( D_{v_2} \): 7 samples (2 class A, 5 class B).

First, calculate the parent entropy (from earlier): \( Entropy(D) \approx 0.9710 \).

Next, calculate the weighted entropy of the children:

Step 1: Define the Parent Impurity

For a dataset \( D \) with \( C \) classes, the parent impurity (using Entropy) is:

Step 2: Split the Dataset on Feature \( A \)

After splitting on feature \( A \), the dataset is divided into subsets \( D_v \) for each value \( v \) of \( A \). The impurity of each subset is:

where \( p_{i,v} \) is the proportion of class \( i \) in subset \( D_v \).

Step 3: Calculate Weighted Child Impurity

The weighted average of the child impurities is:

Step 4: Compute Information Gain

The Information Gain is the difference between the parent impurity and the weighted child impurity:

• Classification Tasks: Decision trees are widely used for classification problems, such as spam detection, customer churn prediction, and medical diagnosis.
• Feature Selection: Information Gain can be used as a feature selection method to identify the most informative features in a dataset.
• Interpretability: Decision trees provide a white-box model, making them useful in domains where interpretability is crucial (e.g., healthcare, finance).
• Handling Non-Linear Relationships: Decision trees can capture non-linear relationships between features and the target variable without requiring feature scaling or transformation.

Derivation: Probability of a Sample Being Out-of-Bag

Consider a dataset with \( N \) samples. In a bootstrap sample, each sample is drawn independently with replacement, so the probability that a specific sample is not selected in one draw is \( 1 - \frac{1}{N} \).

For \( N \) draws, the probability that the sample is not selected at all is: \[ \left(1 - \frac{1}{N}\right)^N \] Taking the limit as \( N \to \infty \): \[ \lim_{N \to \infty} \left(1 - \frac{1}{N}\right)^N = e^{-1} \approx 0.368 \] Thus, roughly \( 36.8\% \) of the data is out-of-bag for each tree.

Derivation: Variance Reduction via Averaging

Let \( T_1(x), T_2(x), \dots, T_B(x) \) be \( B \) i.i.d. trees with variance \( \text{Var}(T_b(x)) = \sigma^2 \). The random forest prediction is the average of these trees: \[ \hat{f}(x) = \frac{1}{B} \sum_{b=1}^B T_b(x) \] The variance of \( \hat{f}(x) \) is: \[ \text{Var}(\hat{f}(x)) = \text{Var}\left(\frac{1}{B} \sum_{b=1}^B T_b(x)\right) = \frac{1}{B^2} \sum_{b=1}^B \text{Var}(T_b(x)) = \frac{\sigma^2}{B} \] This shows that the variance is reduced by a factor of \( B \).

1. Classification Tasks:

• Spam detection: Random forests can classify emails as spam or not spam based on features like word frequency, sender information, etc.
• Medical diagnosis: Predicting diseases (e.g., cancer) from patient data (e.g., age, biomarkers, imaging features).
• Fraud detection: Identifying fraudulent transactions in finance using features like transaction amount, location, and time.

2. Regression Tasks:

• House price prediction: Estimating the price of a house based on features like size, location, and amenities.
• Demand forecasting: Predicting product demand based on historical sales data and external factors (e.g., weather, holidays).
• Stock market analysis: Predicting stock prices or volatility using historical market data.

3. Feature Importance:

Random forests provide a measure of feature importance by calculating the total reduction in a criterion (e.g., Gini impurity or MSE) due to splits on a feature, averaged over all trees. This is useful for:

• Identifying the most relevant features in high-dimensional datasets.
• Feature selection for other models.
• Interpreting model predictions (e.g., in healthcare or finance).

4. Out-of-Bag (OOB) Error for Model Evaluation:

OOB error is a convenient way to estimate the generalization error of a random forest without needing a separate validation set. This is particularly useful when:

• The dataset is small, and splitting it into training and validation sets is not feasible.
• You want to avoid the computational cost of cross-validation.
• You need an unbiased estimate of the test error during model development.

Scikit-Learn Implementation:

from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
from sklearn.datasets import make_classification, make_regression
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, mean_squared_error

# Classification
X_clf, y_clf = make_classification(n_samples=1000, n_features=20, n_classes=2, random_state=42)
X_train_clf, X_test_clf, y_train_clf, y_test_clf = train_test_split(X_clf, y_clf, test_size=0.2, random_state=42)

clf = RandomForestClassifier(
    n_estimators=100,
    max_depth=10,
    min_samples_split=5,
    max_features='sqrt',
    random_state=42,
    oob_score=True
)
clf.fit(X_train_clf, y_train_clf)

y_pred_clf = clf.predict(X_test_clf)
print(f"Test Accuracy: {accuracy_score(y_test_clf, y_pred_clf):.4f}")
print(f"OOB Score: {clf.oob_score_:.4f}")

# Regression
X_reg, y_reg = make_regression(n_samples=1000, n_features=20, noise=0.1, random_state=42)
X_train_reg, X_test_reg, y_train_reg, y_test_reg = train_test_split(X_reg, y_reg, test_size=0.2, random_state=42)

reg = RandomForestRegressor(
    n_estimators=100,
    max_depth=10,
    min_samples_split=5,
    max_features='auto',
    random_state=42,
    oob_score=True
)
reg.fit(X_train_reg, y_train_reg)

y_pred_reg = reg.predict(X_test_reg)
print(f"Test MSE: {mean_squared_error(y_test_reg, y_pred_reg):.4f}")
print(f"OOB Score: {reg.oob_score_:.4f}")

PyTorch Implementation (Conceptual):

PyTorch does not have a built-in random forest implementation, but you can implement a decision tree and extend it to a random forest. Below is a conceptual outline:

import torch
import torch.nn as nn
import numpy as np
from sklearn.tree import DecisionTreeClassifier

class RandomForest(nn.Module):
    def __init__(self, n_estimators=100, max_depth=10, max_features='sqrt'):
        super(RandomForest, self).__init__()
        self.n_estimators = n_estimators
        self.max_depth = max_depth
        self.max_features = max_features
        self.trees = [DecisionTreeClassifier(max_depth=max_depth, max_features=max_features)
                      for _ in range(n_estimators)]

    def fit(self, X, y):
        for tree in self.trees:
            # Bootstrap sampling
            indices = np.random.choice(X.shape[0], X.shape[0], replace=True)
            X_boot = X[indices]
            y_boot = y[indices]
            tree.fit(X_boot, y_boot)

    def predict(self, X):
        predictions = np.array([tree.predict(X) for tree in self.trees])
        # Majority vote for classification
        return np.round(np.mean(predictions, axis=0)).astype(int)

# Example usage
X = np.random.rand(1000, 20)
y = np.random.randint(0, 2, 1000)
rf = RandomForest(n_estimators=100, max_depth=10)
rf.fit(X, y)
y_pred = rf.predict(X[:10])

For a full PyTorch implementation, you would need to implement the decision tree logic (e.g., Gini impurity, splitting criteria) from scratch, which is non-trivial.

Derivation of XGBoost's Tree Splitting Criterion:

XGBoost optimizes the following objective for a tree with \( T \) leaves:

• \( G_j = \sum_{i \in I_j} g_i \): Sum of gradients in leaf \( j \)
• \( H_j = \sum_{i \in I_j} h_i \): Sum of hessians in leaf \( j \)

Step-by-Step Derivation:

1. Start with the objective for a single tree: \[ \text{Obj} = \sum_{i=1}^n L(y_i, F_{m-1}(x_i) + h_m(x_i)) + \Omega(h_m) \]
2. Approximate the loss using a second-order Taylor expansion: \[ L(y_i, F_{m-1}(x_i) + h_m(x_i)) \approx L(y_i, F_{m-1}(x_i)) + g_i h_m(x_i) + \frac{1}{2} h_i h_m(x_i)^2 \]
3. Drop the constant term \( L(y_i, F_{m-1}(x_i)) \) and rewrite the objective: \[ \text{Obj} \approx \sum_{i=1}^n \left( g_i h_m(x_i) + \frac{1}{2} h_i h_m(x_i)^2 \right) + \Omega(h_m) \]
4. For a tree with \( T \) leaves, let \( w_j \) be the weight of leaf \( j \). The objective becomes: \[ \text{Obj} = \sum_{j=1}^T \left( G_j w_j + \frac{1}{2} (H_j + \lambda) w_j^2 \right) + \gamma T \]
5. Take the derivative with respect to \( w_j \) and set to zero to find the optimal weight: \[ \frac{\partial \text{Obj}}{\partial w_j} = G_j + (H_j + \lambda) w_j = 0 \implies w_j^* = -\frac{G_j}{H_j + \lambda} \]
6. Substitute \( w_j^* \) back into the objective to get the splitting criterion: \[ \text{Obj} = -\frac{1}{2} \sum_{j=1}^T \frac{G_j^2}{H_j + \lambda} + \gamma T \]

Derivation of AdaBoost's Weight Update:

AdaBoost minimizes the exponential loss \( L(y, F) = e^{-y F(x)} \). The weight update ensures that misclassified samples receive higher weights in the next iteration.

1. At iteration \( m \), the ensemble model is: \[ F_m(x) = F_{m-1}(x) + \alpha_m h_m(x) \]
2. The exponential loss is: \[ \text{Obj} = \sum_{i=1}^n e^{-y_i F_m(x_i)} = \sum_{i=1}^n e^{-y_i F_{m-1}(x_i)} e^{-y_i \alpha_m h_m(x_i)} \]
3. Let \( w_i^{(m)} = e^{-y_i F_{m-1}(x_i)} \). The objective becomes: \[ \text{Obj} = \sum_{i=1}^n w_i^{(m)} e^{-y_i \alpha_m h_m(x_i)} \]
4. Split the sum into correctly classified (\( y_i h_m(x_i) = 1 \)) and misclassified (\( y_i h_m(x_i) = -1 \)) samples: \[ \text{Obj} = e^{-\alpha_m} \sum_{i: y_i = h_m(x_i)} w_i^{(m)} + e^{\alpha_m} \sum_{i: y_i \neq h_m(x_i)} w_i^{(m)} \]
5. Let \( \epsilon_m = \sum_{i: y_i \neq h_m(x_i)} w_i^{(m)} \). The objective simplifies to: \[ \text{Obj} = e^{-\alpha_m} (1 - \epsilon_m) + e^{\alpha_m} \epsilon_m \]
6. Minimize the objective with respect to \( \alpha_m \): \[ \frac{\partial \text{Obj}}{\partial \alpha_m} = -e^{-\alpha_m} (1 - \epsilon_m) + e^{\alpha_m} \epsilon_m = 0 \] \[ \implies e^{2 \alpha_m} = \frac{1 - \epsilon_m}{\epsilon_m} \implies \alpha_m = \frac{1}{2} \ln \left( \frac{1 - \epsilon_m}{\epsilon_m} \right) \]
7. The weight update rule is derived from \( w_i^{(m+1)} = w_i^{(m)} e^{-y_i \alpha_m h_m(x_i)} \). For misclassified samples (\( y_i h_m(x_i) = -1 \)): \[ w_i^{(m+1)} = w_i^{(m)} e^{\alpha_m} \] For correctly classified samples (\( y_i h_m(x_i) = 1 \)): \[ w_i^{(m+1)} = w_i^{(m)} e^{-\alpha_m} \]

Use Cases for GBMs:

• Tabular Data: GBMs excel with structured/tabular data (e.g., CSV files with numerical and categorical features). Common in finance, healthcare, and marketing.
• Ranking: XGBoost and LightGBM are widely used in learning-to-rank tasks (e.g., search engines, recommendation systems).
• Anomaly Detection: GBMs can identify outliers by modeling the residual errors of normal data points.
• Feature Importance: GBMs provide interpretable feature importance scores, useful for understanding model decisions.
• Competitions: XGBoost and LightGBM are popular in Kaggle competitions due to their high performance and speed.

Example: Training XGBoost in Python (Scikit-Learn API):

from xgboost import XGBClassifier
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split

# Load data
data = load_breast_cancer()
X, y = data.data, data.target
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)

# Initialize and train XGBoost
model = XGBClassifier(
    n_estimators=100,
    learning_rate=0.1,
    max_depth=3,
    subsample=0.8,
    colsample_bytree=0.8,
    reg_alpha=0.1,
    reg_lambda=1.0,
    random_state=42
)
model.fit(X_train, y_train)

# Evaluate
accuracy = model.score(X_test, y_test)
print(f"Test Accuracy: {accuracy:.4f}")

Example: Training LightGBM with Categorical Features:

import lightgbm as lgb
import pandas as pd

# Load data with categorical features
data = pd.read_csv("data.csv")
X = data.drop("target", axis=1)
y = data["target"]

# Convert categorical columns to 'category' dtype
categorical_features = ["cat_feature1", "cat_feature2"]
for col in categorical_features:
    X[col] = X[col].astype("category")

# Create LightGBM dataset
train_data = lgb.Dataset(X, label=y, categorical_feature=categorical_features)

# Define parameters
params = {
    "objective": "binary",
    "metric": "binary_logloss",
    "boosting_type": "gbdt",
    "num_leaves": 31,
    "learning_rate": 0.05,
    "feature_fraction": 0.9,
    "bagging_fraction": 0.8,
    "bagging_freq": 5,
    "verbose": -1
}

# Train
model = lgb.train(params, train_data, num_boost_round=100)

Example: Training Gaussian Naive Bayes

Given a dataset with features \(\mathbf{X} = [x_1, x_2]\) and class labels \(y \in \{0, 1\}\), the steps to train the model are:

1. Compute the prior probabilities \(P(y=0)\) and \(P(y=1)\).
2. For each feature \(x_i\) and class \(y\), compute the mean \(\mu_{y,i}\) and variance \(\sigma_{y,i}^2\) of the feature values for that class.

For prediction, compute the posterior probability for each class using the Gaussian PDF and select the class with the highest probability.

Example: Text Classification with Multinomial Naive Bayes

Consider a dataset of documents labeled as "spam" or "not spam". Each document is represented as a bag-of-words vector \(\mathbf{x} = [x_1, x_2, \dots, x_n]\), where \(x_i\) is the count of word \(i\) in the document.

1. Compute the prior probabilities \(P(y=\text{spam})\) and \(P(y=\text{not spam})\).
2. For each word \(i\) and class \(y\), compute \(\theta_{y,i}\) (the probability of word \(i\) given class \(y\)) using Laplace smoothing.
3. For a new document, compute the posterior probability for each class and assign the class with the highest probability.

Example: Binary Text Classification with Bernoulli Naive Bayes

Consider a dataset of documents where each document is represented as a binary vector \(\mathbf{x} = [x_1, x_2, \dots, x_n]\), where \(x_i = 1\) if word \(i\) is present in the document and \(x_i = 0\) otherwise.

1. Compute the prior probabilities \(P(y=\text{spam})\) and \(P(y=\text{not spam})\).
2. For each word \(i\) and class \(y\), compute \(\theta_{y,i}\) (the probability of word \(i\) being present in class \(y\)) using Laplace smoothing.
3. For a new document, compute the posterior probability for each class and assign the class with the highest probability.

• Gaussian Naive Bayes: Medical diagnosis (e.g., classifying diseases based on continuous test results), anomaly detection, and real-valued sensor data.
• Multinomial Naive Bayes: Text classification (e.g., spam detection, sentiment analysis), document categorization, and recommendation systems.
• Bernoulli Naive Bayes: Binary text classification (e.g., presence/absence of keywords), author identification, and multi-label classification tasks.

Scikit-Learn Implementation:

from sklearn.naive_bayes import GaussianNB, MultinomialNB, BernoulliNB

# Gaussian Naive Bayes
gnb = GaussianNB()
gnb.fit(X_train, y_train)
y_pred = gnb.predict(X_test)

# Multinomial Naive Bayes
mnb = MultinomialNB(alpha=1.0)  # alpha is the smoothing parameter
mnb.fit(X_train, y_train)
y_pred = mnb.predict(X_test)

# Bernoulli Naive Bayes
bnb = BernoulliNB(alpha=1.0, binarize=0.5)  # binarize threshold for feature values
bnb.fit(X_train, y_train)
y_pred = bnb.predict(X_test)

Key Parameters in Scikit-Learn:

• alpha: Smoothing parameter (default=1.0).
• binarize (BernoulliNB): Threshold for binarizing features (default=0.0).
• fit_prior: Whether to learn class prior probabilities (default=True).

PyTorch Implementation (Custom Gaussian Naive Bayes):

While PyTorch does not have built-in Naive Bayes implementations, you can implement a custom Gaussian Naive Bayes model as follows:

import torch

class GaussianNaiveBayes:
    def __init__(self):
        self.classes_ = None
        self.mean_ = None
        self.var_ = None
        self.priors_ = None

    def fit(self, X, y):
        self.classes_ = torch.unique(y)
        n_classes = len(self.classes_)
        n_features = X.shape[1]

        self.mean_ = torch.zeros((n_classes, n_features))
        self.var_ = torch.zeros((n_classes, n_features))
        self.priors_ = torch.zeros(n_classes)

        for i, c in enumerate(self.classes_):
            X_c = X[y == c]
            self.mean_[i, :] = X_c.mean(dim=0)
            self.var_[i, :] = X_c.var(dim=0, unbiased=False)
            self.priors_[i] = X_c.shape[0] / X.shape[0]

    def predict(self, X):
        log_probs = []
        for i, c in enumerate(self.classes_):
            prior = torch.log(self.priors_[i])
            likelihood = -0.5 * torch.sum(torch.log(2. * torch.pi * self.var_[i, :]) +
                                         ((X - self.mean_[i, :]) ** 2) / self.var_[i, :], dim=1)
            log_prob = prior + likelihood
            log_probs.append(log_prob)

        log_probs = torch.stack(log_probs, dim=1)
        return self.classes_[torch.argmax(log_probs, dim=1)]

Note: This implementation computes log probabilities to avoid numerical underflow.

Step 1: Center the Data

Subtract the mean of each feature from the data to center it around the origin:

where \( \mu \in \mathbb{R}^{1 \times d} \) is the mean vector of the features.

Step 2: Compute the Covariance Matrix

Calculate the covariance matrix \( \Sigma \) as shown above. This matrix captures the relationships between features.

Step 3: Perform Eigenvalue Decomposition

Decompose \( \Sigma \) into its eigenvalues and eigenvectors:

where \( V \) is the matrix of eigenvectors (principal components) and \( \Lambda \) is the diagonal matrix of eigenvalues. The eigenvectors are sorted in descending order of their corresponding eigenvalues.

Step 4: Project the Data

Project the centered data \( X_{\text{centered}} \) onto the principal components to obtain the transformed data \( Z \):

where \( V_k \) contains the first \( k \) eigenvectors (columns of \( V \)).

Step 5: Compute Variance Explained

Calculate the variance explained by each principal component using the eigenvalues, as shown in the formulas above.

1. Dimensionality Reduction: PCA is widely used to reduce the number of features in a dataset while preserving as much variance as possible. This is useful for visualization (e.g., reducing to 2D or 3D) and speeding up downstream tasks like classification or regression.

2. Noise Reduction: By retaining only the principal components with the highest variance, PCA can filter out noise in the data, as noise typically contributes less to the variance.

3. Feature Extraction: PCA can transform the original features into a new set of uncorrelated features (principal components), which can improve the performance of machine learning models.

4. Anomaly Detection: Data points that lie far from the principal components (low variance directions) can be flagged as anomalies.

5. Data Compression: PCA can compress high-dimensional data (e.g., images) by storing only the principal components and their projections.

Below is a Python example using Scikit-Learn to perform PCA on the Iris dataset:

from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

# Load data
data = load_iris()
X = data.data
y = data.target

# Standardize the data
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Apply PCA
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X_scaled)

# Plot the results
plt.figure(figsize=(8, 6))
plt.scatter(X_pca[:, 0], X_pca[:, 1], c=y, cmap='viridis')
plt.xlabel('First Principal Component')
plt.ylabel('Second Principal Component')
plt.title('PCA of Iris Dataset')
plt.colorbar()
plt.show()

# Print variance explained
print("Variance explained by each component:", pca.explained_variance_ratio_)
print("Total variance explained:", sum(pca.explained_variance_ratio_))

Output:

• The plot shows the Iris dataset projected onto the first two principal components.
• The variance explained by each component is printed, e.g., [0.7296, 0.2285], meaning the first PC explains ~73% of the variance, and the second PC explains ~23%.

Below is a manual implementation of PCA using eigenvalue decomposition in NumPy:

import numpy as np

# Center the data
X_centered = X - np.mean(X, axis=0)

# Compute covariance matrix
cov_matrix = np.cov(X_centered, rowvar=False)

# Perform eigenvalue decomposition
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)

# Sort eigenvectors by eigenvalues (descending order)
sorted_indices = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[sorted_indices]
eigenvectors = eigenvectors[:, sorted_indices]

# Project data onto first 2 principal components
X_pca_manual = X_centered @ eigenvectors[:, :2]

# Print variance explained
variance_explained = eigenvalues / np.sum(eigenvalues)
print("Variance explained by each component:", variance_explained[:2])
print("Total variance explained:", sum(variance_explained[:2]))

Note: This manual implementation matches the Scikit-Learn output, demonstrating the underlying mathematics.

Derivation of Low-Rank Approximation:

1. Eckart-Young Theorem: The best rank-\( k \) approximation of \( A \) in the Frobenius norm is given by \( A_k \), where: \[ \| A - A_k \|_F = \min_{\text{rank}(B) \leq k} \| A - B \|_F = \sqrt{\sigma_{k+1}^2 + \dots + \sigma_{\min(m,n)}^2} \]
2. Truncated SVD: To compute \( A_k \), retain only the top \( k \) singular values and their corresponding singular vectors: \[ A_k = \sum_{i=1}^k \sigma_i u_i v_i^T \] where \( u_i \) and \( v_i \) are the \( i \)-th columns of \( U \) and \( V \), respectively.

Worked Example: Let \( A \) be a \( 4 \times 3 \) matrix with SVD:

Suppose \( U = I_4 \) (identity matrix) and \( V = I_3 \). The singular values are \( \sigma_1 = 3 \), \( \sigma_2 = 2 \), \( \sigma_3 = 1 \).

For \( k = 2 \), the low-rank approximation \( A_2 \) is:

The Frobenius norm error is:

1. Dimensionality Reduction (PCA):

   Principal Component Analysis (PCA) is a linear dimensionality reduction technique that uses SVD. Given a data matrix \( X \) (centered), its SVD is \( X = U \Sigma V^T \). The principal components are the columns of \( V \), and the projected data is \( U \Sigma \). Truncating to the top \( k \) singular values yields the best \( k \)-dimensional approximation of the data.

2. Image Compression:

   Images can be represented as matrices. By computing the SVD of an image matrix and retaining only the top \( k \) singular values, the image can be compressed with minimal loss of quality. The storage required is reduced from \( O(mn) \) to \( O(k(m + n)) \).

3. Latent Semantic Indexing (LSI):

   In natural language processing, LSI uses SVD to identify patterns in the relationships between terms and concepts in unstructured text. The term-document matrix is decomposed, and low-rank approximation is used to capture latent semantic structure.

4. Recommender Systems:

   SVD is used in collaborative filtering to factorize the user-item interaction matrix into latent factors. The low-rank approximation helps predict missing entries (e.g., user ratings) and make recommendations.

5. Noise Reduction:

   By truncating small singular values (which often correspond to noise), SVD can denoise data. This is useful in signal processing and image restoration.

6. Solving Linear Systems:

   For underdetermined or overdetermined systems \( Ax = b \), SVD can be used to compute the least-squares solution or the minimum-norm solution via the pseudoinverse.

PyTorch:

import torch

# Create a random matrix
A = torch.randn(4, 3)

# Compute SVD
U, S, V = torch.svd(A)

# Low-rank approximation (k=2)
k = 2
U_k = U[:, :k]
S_k = torch.diag(S[:k])
V_k = V[:, :k]
A_k = U_k @ S_k @ V_k.t()

print("Original matrix:\n", A)
print("Low-rank approximation:\n", A_k)

Scikit-Learn (for PCA):

from sklearn.decomposition import TruncatedSVD
import numpy as np

# Create a random matrix
X = np.random.rand(100, 10)  # 100 samples, 10 features

# Apply TruncatedSVD for dimensionality reduction (k=3)
svd = TruncatedSVD(n_components=3)
X_reduced = svd.fit_transform(X)

print("Original shape:", X.shape)
print("Reduced shape:", X_reduced.shape)
print("Explained variance ratio:", svd.explained_variance_ratio_)

Q1: What is the difference between SVD and PCA?

A: PCA is a dimensionality reduction technique that uses SVD as its computational backbone. Specifically, PCA involves centering the data matrix \( X \) and then computing its SVD: \( X = U \Sigma V^T \). The principal components are the columns of \( V \), and the projected data is \( U \Sigma \). SVD is a more general matrix factorization technique that can be applied to any matrix, while PCA is a specific application of SVD to data analysis.

Q2: How do you choose the rank \( k \) for low-rank approximation?

A: The choice of \( k \) depends on the application and the trade-off between approximation error and computational efficiency. Common methods include:

• Retaining singular values above a certain threshold (e.g., \( \sigma_i > \epsilon \)).
• Choosing \( k \) such that a certain fraction of the total variance is preserved (e.g., 95%). The total variance is the sum of squares of the singular values, and the preserved variance is the sum of squares of the top \( k \) singular values.
• Using the "elbow method" to identify a knee point in the singular value spectrum, where adding more components yields diminishing returns.

Q3: What is the relationship between SVD and the Moore-Penrose pseudoinverse?

A: The Moore-Penrose pseudoinverse \( A^+ \) of a matrix \( A \) can be computed using its SVD. If \( A = U \Sigma V^T \), then:

where \( \Sigma^+ \) is obtained by taking the reciprocal of each non-zero singular value in \( \Sigma \) and transposing the resulting matrix. The pseudoinverse is used to solve linear systems \( Ax = b \) in the least-squares sense when \( A \) is not square or is rank-deficient.

Q4: Why is SVD useful for recommender systems?

A: In recommender systems, the user-item interaction matrix (e.g., user ratings) is often sparse and incomplete. SVD can factorize this matrix into latent factors representing users and items. The low-rank approximation helps predict missing entries by capturing underlying patterns in the data. This is the basis for collaborative filtering techniques like FunkSVD.

Q5: How does SVD help in noise reduction?

A: In many applications, small singular values correspond to noise in the data, while larger singular values capture the signal. By truncating the SVD and retaining only the top \( k \) singular values, the reconstructed matrix \( A_k \) will have reduced noise. This is because the discarded singular values (and their corresponding singular vectors) contribute less to the overall structure of the data.

Example: EM Algorithm for GMM

Consider a 1D dataset with \( N = 1000 \) points generated from a mixture of two Gaussians. Initialize \( K = 2 \), \( \pi_1 = \pi_2 = 0.5 \), \( \mu_1 = 0 \), \( \mu_2 = 1 \), and \( \Sigma_1 = \Sigma_2 = 1 \).

E-Step:

For each data point \( x_n \), compute the responsibilities:

M-Step:

Update parameters:

Iterate until convergence (e.g., change in log-likelihood is below a threshold).

Example: Model Selection with AIC/BIC

Suppose we fit GMMs with \( K = 1, 2, 3, 4 \) to a dataset and obtain the following log-likelihoods and parameter counts:

The model with \( K = 3 \) has the lowest AIC, while \( K = 2 \) has the lowest BIC (due to the stronger penalty for complexity). The choice depends on whether you prioritize fit (AIC) or simplicity (BIC).

Ward's method aims to minimize the increase in total within-cluster variance when merging two clusters. The within-cluster variance for a cluster \( C \) is:

where \( \bar{x} \) is the centroid of \( C \). The increase in variance when merging \( C_i \) and \( C_j \) is:

where \( C_k = C_i \cup C_j \). Using the identity for the variance of merged clusters:

Thus, the increase in variance is:

Ward's distance is the square root of this increase, scaled by 2 for consistency with other linkage methods:

• Biology: Hierarchical clustering is widely used in genomics for gene expression analysis, where it helps identify groups of genes with similar expression patterns.
• Document Clustering: Used in natural language processing to group similar documents (e.g., news articles or research papers) based on their content.
• Image Segmentation: Hierarchical clustering can segment images into regions with similar pixel intensities or textures.
• Customer Segmentation: Businesses use hierarchical clustering to group customers based on purchasing behavior or demographic data.
• Phylogenetics: Used to construct phylogenetic trees that represent evolutionary relationships between species.

from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram, linkage
import matplotlib.pyplot as plt
import numpy as np

# Sample data
X = np.array([[1, 2], [1, 4], [1, 0], [4, 2], [4, 4], [4, 0]])

# Agglomerative Clustering
clustering = AgglomerativeClustering(n_clusters=2, affinity='euclidean', linkage='ward')
clustering.fit(X)
print("Cluster labels:", clustering.labels_)

# Dendrogram
Z = linkage(X, method='ward')
plt.figure(figsize=(10, 5))
dendrogram(Z)
plt.title('Dendrogram')
plt.show()
    

Key Parameters:

• n_clusters: Number of clusters to find (None for hierarchical representation).
• affinity: Distance metric (e.g., 'euclidean', 'manhattan').
• linkage: Linkage criterion ('ward', 'complete', 'average', 'single').

1. What is the difference between agglomerative and divisive hierarchical clustering?

   Answer: Agglomerative clustering is a bottom-up approach where each data point starts in its own cluster, and clusters are merged iteratively. Divisive clustering is a top-down approach where all data points start in one cluster, and the cluster is recursively split into smaller clusters.

2. How do you choose the number of clusters in hierarchical clustering?

   Answer: The number of clusters is typically chosen by inspecting the dendrogram and cutting it at a height that yields a desired number of clusters. Alternatively, domain knowledge or metrics like the elbow method (for within-cluster variance) can be used.

3. What are the advantages and disadvantages of single linkage vs. complete linkage?

   Answer:

   • Single Linkage:
     • Advantages: Can detect non-elliptical clusters and is less sensitive to outliers.
     • Disadvantages: Prone to chaining, which can lead to long, straggly clusters.
   • Complete Linkage:
     • Advantages: Tends to produce compact, spherical clusters.
     • Disadvantages: Sensitive to outliers and may not perform well with non-spherical clusters.

4. Explain Ward's method for hierarchical clustering.

   Answer: Ward's method minimizes the increase in total within-cluster variance when merging two clusters. It is equivalent to minimizing the sum of squared distances between points and their cluster centroids. The distance between two clusters is calculated as the square root of the increase in variance when they are merged.

5. How does the Lance-Williams formula generalize linkage criteria?

   Answer: The Lance-Williams formula provides a unified way to update distances between clusters after a merge. It expresses the distance between a new cluster \( C_k \) (formed by merging \( C_i \) and \( C_j \)) and another cluster \( C_l \) as a weighted combination of the distances \( D(C_i, C_l) \), \( D(C_j, C_l) \), and \( D(C_i, C_j) \). The weights depend on the linkage criterion (e.g., single, complete, average, Ward).

Example: DBSCAN Clustering Process

1. Initialize: Mark all points as unvisited.
2. Iterate: For each unvisited point \( p \):
   • Mark \( p \) as visited.
   • Find \( N_\epsilon(p) \).
   • If \( |N_\epsilon(p)| < \text{MinPts} \), mark \( p \) as noise (temporarily).
   • Else:
     • Create a new cluster \( C \) and add \( p \) to \( C \).
     • For each point \( q \) in \( N_\epsilon(p) \):
       • If \( q \) is unvisited, mark it as visited and find \( N_\epsilon(q) \). If \( |N_\epsilon(q)| \geq \text{MinPts} \), add \( N_\epsilon(q) \) to the seed set.
       • If \( q \) is not yet a member of any cluster, add \( q \) to \( C \).
3. Terminate: When all points are visited, the algorithm terminates.

Worked Example:

Consider the following 2D dataset with MinPts = 3 and ε = 2:

Points: A(1,1), B(1.5,1.5), C(5,5), D(5.5,5.5), E(6,6), F(10,10), G(3,4), H(4,3)
    

1. Start with point A:
   • \( N_2(A) = \{A, B\} \) (size = 2 < 3) → Mark A as noise (temporarily).
2. Visit point B:
   • \( N_2(B) = \{A, B, G\} \) (size = 3 ≥ 3) → Core point. Create cluster C₁ and add B.
   • Add A and G to C₁ (both are border points).
   • For G: \( N_2(G) = \{B, G, H\} \) (size = 3 ≥ 3) → Core point. Add H to C₁.
3. Visit point C:
   • \( N_2(C) = \{C, D, E\} \) (size = 3 ≥ 3) → Core point. Create cluster C₂ and add C.
   • Add D and E to C₂.
4. Visit point F:
   • \( N_2(F) = \{F\} \) (size = 1 < 3) → Mark F as noise.
5. Final Clusters: C₁ = {A, B, G, H}, C₂ = {C, D, E}, Noise: {F}.

PyTorch Implementation Note:

While DBSCAN is not natively implemented in PyTorch, you can use PyTorch for distance computations and then apply DBSCAN logic. Here’s a minimal example:

import torch
from sklearn.neighbors import NearestNeighbors

# Generate synthetic data
X = torch.randn(100, 2)

# Compute pairwise distances (PyTorch)
distances = torch.cdist(X, X)

# Use scikit-learn for DBSCAN (or implement custom logic)
from sklearn.cluster import DBSCAN
dbscan = DBSCAN(eps=0.5, min_samples=5, metric='precomputed')
clusters = dbscan.fit_predict(distances.numpy())
    

Example: Suppose we have true values \( \mathbf{y} = [3, -0.5, 2] \) and predicted values \( \mathbf{\hat{y}} = [2.5, 0.0, 2.1] \). The MSE is calculated as:

Example (Binary Cross-Entropy): For a single sample with true label \( y = 1 \) and predicted probability \( \hat{y} = 0.9 \):

For \( y = 0 \) and \( \hat{y} = 0.1 \):

Incorrect predictions (e.g., \( y = 1, \hat{y} = 0.1 \)) yield higher loss: \( -\log(0.1) \approx 2.3026 \).

Example: Consider two samples:

• Sample 1: \( y_1 = 1 \), \( \hat{y}_1 = 1.5 \) → \( \max(0, 1 - 1 \cdot 1.5) = \max(0, -0.5) = 0 \)
• Sample 2: \( y_2 = -1 \), \( \hat{y}_2 = 0.3 \) → \( \max(0, 1 - (-1) \cdot 0.3) = \max(0, 1.3) = 1.3 \)

The hinge loss for these samples is \( \frac{0 + 1.3}{2} = 0.65 \).

Example (Discrete): Let \( P = [0.6, 0.4] \) and \( Q = [0.5, 0.5] \). Then:

Note that \( D_{KL}(Q \parallel P) \approx 0.0204 \), illustrating the asymmetry.

SGD: Often used in large-scale machine learning problems where the dataset is too large to compute the full gradient. It is simple to implement and works well with a properly tuned learning rate.

Momentum: Helps accelerate SGD in the relevant direction and dampens oscillations. It is particularly useful in cases where the loss function has high curvature or noisy gradients.

Adam: Widely used in deep learning due to its adaptive learning rate properties. It is particularly effective for problems with sparse gradients and non-stationary objectives.

RMSprop: Effective for recurrent neural networks (RNNs) and problems with non-convex loss landscapes. It helps in handling the vanishing and exploding gradient problems.

Learning Rate Schedules: Useful for fine-tuning the learning process. Step decay is commonly used in training deep neural networks, while 1Cycle policy has been shown to achieve faster convergence and better performance in some cases.

Derivation of Output Size for Convolution:

1. Start with an input of size \(W \times W\).
2. Add padding \(P\) to each side, increasing the effective input size to \((W + 2P) \times (W + 2P)\).
3. The kernel of size \(K \times K\) slides over the padded input with stride \(S\).
4. The number of possible positions the kernel can take along one dimension is: \[ \frac{(W + 2P) - K}{S} + 1 \] The floor function \(\lfloor \cdot \rfloor\) is applied to ensure the result is an integer.

Effect of Stride and Padding:

• Stride = 1, Padding = 0 (Valid Convolution): \[ \text{Output Size} = \left\lfloor \frac{5 - 3 + 0}{1} \right\rfloor + 1 = 3 \] For a \(5 \times 5\) input and \(3 \times 3\) kernel.
• Stride = 2, Padding = 1 (Same Convolution): \[ \text{Output Size} = \left\lfloor \frac{5 - 3 + 2}{2} \right\rfloor + 1 = 3 \] The output size matches the input size (\(5 \times 5\)) due to padding.

Image Classification: CNNs are the backbone of modern image classification tasks (e.g., ResNet, VGG). Kernels learn to detect edges, textures, and patterns, while pooling reduces spatial dimensions to focus on high-level features.

Object Detection: CNNs (e.g., YOLO, Faster R-CNN) use convolutional layers to generate feature maps for detecting and localizing objects in images.

Semantic Segmentation: Architectures like U-Net use CNNs to classify each pixel in an image, enabling applications like medical image analysis and autonomous driving.

Natural Language Processing (NLP): 1D CNNs are used for text classification (e.g., sentiment analysis) by treating sequences as 1D grids.

Vanishing Gradients Example:

Consider a sequence of length \( T = 100 \) and \( W_{hh} \) with singular values \( \sigma \approx 0.9 \). The gradient term becomes:

This demonstrates how gradients vanish exponentially with sequence length.

Exploding Gradients Example:

If \( W_{hh} \) has a singular value \( \sigma \approx 1.1 \), the gradient term becomes:

This leads to numerical instability and overflow.

LSTM Gradient Flow:

The cell state \( C_t \) in LSTMs allows gradients to flow unchanged through time:

By setting \( f_t \approx 1 \), the gradient can propagate over long sequences without vanishing.

Consider a sequence of 3 tokens, each embedded into a 4-dimensional space. The input embeddings are:

We use learned weight matrices \(W_Q, W_K, W_V \in \mathbb{R}^{4 \times 4}\) to project \(X\) into queries, keys, and values:

Assume \(W_Q = W_K = W_V = I\) (identity matrix) for simplicity. Then:

Compute the attention scores:

Scale by \(\sqrt{d_k} = \sqrt{4} = 2\):

Apply softmax to each row:

Finally, compute the output:

The dot product \(QK^T\) grows with the dimension \(d_k\). For large \(d_k\), the dot products can become very large, pushing the softmax into regions with extremely small gradients. To counteract this, the dot product is scaled by \(\sqrt{d_k}\):

Assume \(q\) and \(k\) are random vectors with mean 0 and variance 1. The dot product \(q \cdot k\) has mean 0 and variance \(d_k\). Thus, the dot product grows as \(O(\sqrt{d_k})\), and scaling by \(\sqrt{d_k}\) ensures the variance remains 1, stabilizing training.

Multi-head attention allows the model to jointly attend to information from different representation subspaces. For example, one head might focus on syntactic relationships, while another captures semantic dependencies. This is analogous to having multiple "experts" specializing in different aspects of the data.

Derivation of the ELBO

The goal of variational inference is to maximize the log-likelihood of the observed data \(\log p_\theta(\mathbf{x})\). However, this is intractable for complex models. Instead, we maximize a lower bound on the log-likelihood, known as the Evidence Lower Bound (ELBO).

1. Start with the log-likelihood:

   \[ \log p_\theta(\mathbf{x}) = \log \int p_\theta(\mathbf{x}, \mathbf{z}) d\mathbf{z} \]

2. Introduce the approximate posterior \(q_\phi(\mathbf{z}|\mathbf{x})\):

   \[ \log p_\theta(\mathbf{x}) = \log \int q_\phi(\mathbf{z}|\mathbf{x}) \frac{p_\theta(\mathbf{x}, \mathbf{z})}{q_\phi(\mathbf{z}|\mathbf{x})} d\mathbf{z} \]

3. Apply Jensen's inequality (since \(\log\) is concave):

   \[ \log p_\theta(\mathbf{x}) \geq \int q_\phi(\mathbf{z}|\mathbf{x}) \log \frac{p_\theta(\mathbf{x}, \mathbf{z})}{q_\phi(\mathbf{z}|\mathbf{x})} d\mathbf{z} = \mathcal{L}(\theta, \phi; \mathbf{x}) \]

4. Rewrite the ELBO:

   \[ \mathcal{L}(\theta, \phi; \mathbf{x}) = \mathbb{E}_{q_\phi(\mathbf{z}|\mathbf{x})}[\log p_\theta(\mathbf{x}|\mathbf{z})] - \text{KL}(q_\phi(\mathbf{z}|\mathbf{x}) \| p(\mathbf{z})) \]

Reparameterization Trick Derivation

The reparameterization trick allows gradients to flow through the stochastic sampling step in the VAE. Here's how it works:

1. Assume \(q_\phi(\mathbf{z}|\mathbf{x}) = \mathcal{N}(\mu, \sigma^2)\). Sampling \(\mathbf{z}\) directly from this distribution is not differentiable with respect to \(\mu\) and \(\sigma\).

2. Instead, express \(\mathbf{z}\) as a deterministic function of \(\mu\), \(\sigma\), and a random noise variable \(\epsilon \sim \mathcal{N}(0, I)\):

   \[ \mathbf{z} = \mu + \sigma \odot \epsilon \]

3. Now, the gradient of \(\mathbf{z}\) with respect to \(\mu\) and \(\sigma\) can be computed, as \(\epsilon\) is independent of \(\mu\) and \(\sigma\).

4. This reparameterization allows the use of standard backpropagation to train the VAE.

1. Anomaly Detection

VAEs can learn a compressed representation of "normal" data. During inference, if a new data point has a high reconstruction error, it is likely an anomaly. This is useful in fraud detection, manufacturing defect detection, and medical diagnosis.

2. Data Generation

VAEs can generate new data samples by sampling from the latent space and passing the samples through the decoder. This is used in applications like image generation, text generation, and drug discovery.

3. Dimensionality Reduction

VAEs can be used for non-linear dimensionality reduction, similar to PCA but with the ability to capture more complex data structures. The latent space can be used for visualization or as features for downstream tasks.

4. Denoising

VAEs can be trained to reconstruct clean data from noisy inputs, making them useful for image denoising, speech enhancement, and other signal processing tasks.

5. Semi-Supervised Learning

VAEs can be extended to semi-supervised learning tasks, where the model leverages both labeled and unlabeled data to improve performance on tasks like classification.

VAE Model Architecture

import torch
import torch.nn as nn
import torch.nn.functional as F

class VAE(nn.Module):
    def __init__(self, input_dim, hidden_dim, latent_dim):
        super(VAE, self).__init__()
        # Encoder
        self.fc1 = nn.Linear(input_dim, hidden_dim)
        self.fc_mu = nn.Linear(hidden_dim, latent_dim)
        self.fc_logvar = nn.Linear(hidden_dim, latent_dim)

        # Decoder
        self.fc2 = nn.Linear(latent_dim, hidden_dim)
        self.fc3 = nn.Linear(hidden_dim, input_dim)

    def encode(self, x):
        h = F.relu(self.fc1(x))
        return self.fc_mu(h), self.fc_logvar(h)

    def reparameterize(self, mu, logvar):
        std = torch.exp(0.5 * logvar)
        eps = torch.randn_like(std)
        return mu + eps * std

    def decode(self, z):
        h = F.relu(self.fc2(z))
        return torch.sigmoid(self.fc3(h))

    def forward(self, x):
        mu, logvar = self.encode(x)
        z = self.reparameterize(mu, logvar)
        return self.decode(z), mu, logvar
    

Loss Function

The VAE loss is the negative ELBO, which consists of the reconstruction loss and the KL divergence.

def vae_loss(recon_x, x, mu, logvar):
    # Reconstruction loss (binary cross-entropy)
    BCE = F.binary_cross_entropy(recon_x, x, reduction='sum')

    # KL divergence
    KLD = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())

    return BCE + KLD
    

Training Loop

model = VAE(input_dim=784, hidden_dim=400, latent_dim=20)
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)

for epoch in range(epochs):
    for batch_idx, (data, _) in enumerate(train_loader):
        data = data.view(-1, 784)  # Flatten the data
        optimizer.zero_grad()
        recon_batch, mu, logvar = model(data)
        loss = vae_loss(recon_batch, data, mu, logvar)
        loss.backward()
        optimizer.step()
    

Example: Q-Learning in a Grid World

Consider a 2x2 grid world where the agent starts at the top-left corner and the goal is to reach the bottom-right corner. The agent can move up, down, left, or right. Each step incurs a reward of -1, except reaching the goal which gives a reward of +10.

Initialize \(Q(s,a)\) arbitrarily (e.g., to zero). For each episode:

1. Choose an action \(a_t\) in state \(s_t\) using an exploration strategy (e.g., \(\epsilon\)-greedy).
2. Observe the reward \(r_{t+1}\) and next state \(s_{t+1}\).
3. Update \(Q(s_t, a_t)\) using the Q-learning update rule.
4. Repeat until \(s_{t+1}\) is the terminal state.

Example: REINFORCE for CartPole

In the CartPole environment, the agent must balance a pole on a cart by moving left or right. The policy \(\pi_\theta(a|s)\) can be represented by a neural network with parameters \(\theta\). The steps are:

1. Initialize the policy parameters \(\theta\) randomly.
2. Generate an episode by following \(\pi_\theta(a|s)\).
3. For each step \(t\) in the episode, compute the return \(G_t\).
4. Update \(\theta\) using the REINFORCE update rule.
5. Repeat for multiple episodes.

Example: A2C for LunarLander

In the LunarLander environment, the agent must land a spacecraft on a landing pad. The actor-critic method can be implemented as follows:

1. Initialize the actor (\(\pi_\theta\)) and critic (\(V_w\)) networks.
2. For each episode:
1. Sample an action \(a_t \sim \pi_\theta(\cdot|s_t)\).
2. Observe the reward \(r_{t+1}\) and next state \(s_{t+1}\).
3. Compute the TD error \(\delta_t = r_{t+1} + \gamma V_w(s_{t+1}) - V_w(s_t)\).
4. Update the critic: \(w \leftarrow w + \beta \delta_t \nabla_w V_w(s_t)\).
5. Update the actor: \(\theta \leftarrow \theta + \alpha \nabla_\theta \log \pi_\theta(a_t|s_t) \delta_t\).
3. Repeat until convergence.

Applications of Reinforcement Learning:

• Robotics: Training robots to perform tasks like grasping objects, walking, or navigating environments (e.g., using DDPG or PPO).
• Game Playing: Achieving superhuman performance in games like Go (AlphaGo), Chess (AlphaZero), or video games (DQN for Atari).
• Autonomous Vehicles: Decision-making for self-driving cars, including lane-keeping, obstacle avoidance, and route planning.
• Finance: Algorithmic trading, portfolio management, and risk assessment.
• Healthcare: Personalized treatment planning, drug discovery, and resource allocation in hospitals.
• Recommendation Systems: Dynamic recommendation of content or products based on user interactions.

Q-Learning with PyTorch (DQN):

import torch
import torch.nn as nn
import torch.optim as optim
import numpy as np
from collections import deque
import random

class DQN(nn.Module):
    def __init__(self, state_dim, action_dim):
        super(DQN, self).__init__()
        self.fc1 = nn.Linear(state_dim, 64)
        self.fc2 = nn.Linear(64, 64)
        self.fc3 = nn.Linear(64, action_dim)

    def forward(self, x):
        x = torch.relu(self.fc1(x))
        x = torch.relu(self.fc2(x))
        return self.fc3(x)

class DQNAgent:
    def __init__(self, state_dim, action_dim):
        self.model = DQN(state_dim, action_dim)
        self.target_model = DQN(state_dim, action_dim)
        self.target_model.load_state_dict(self.model.state_dict())
        self.optimizer = optim.Adam(self.model.parameters(), lr=0.001)
        self.memory = deque(maxlen=10000)
        self.batch_size = 64
        self.gamma = 0.99
        self.epsilon = 1.0
        self.epsilon_min = 0.01
        self.epsilon_decay = 0.995

    def remember(self, state, action, reward, next_state, done):
        self.memory.append((state, action, reward, next_state, done))

    def act(self, state):
        if np.random.rand() <= self.epsilon:
            return random.randrange(action_dim)
        state = torch.FloatTensor(state).unsqueeze(0)
        q_values = self.model(state)
        return torch.argmax(q_values).item()

    def replay(self):
        if len(self.memory) < self.batch_size:
            return
        batch = random.sample(self.memory, self.batch_size)
        states, actions, rewards, next_states, dones = zip(*batch)

        states = torch.FloatTensor(np.array(states))
        actions = torch.LongTensor(actions).unsqueeze(1)
        rewards = torch.FloatTensor(rewards)
        next_states = torch.FloatTensor(np.array(next_states))
        dones = torch.FloatTensor(dones)

        current_q = self.model(states).gather(1, actions)
        next_q = self.target_model(next_states).max(1)[0].detach()
        target_q = rewards + (1 - dones) * self.gamma * next_q

        loss = nn.MSELoss()(current_q.squeeze(), target_q)
        self.optimizer.zero_grad()
        loss.backward()
        self.optimizer.step()

        if self.epsilon > self.epsilon_min:
            self.epsilon *= self.epsilon_decay

    def update_target_model(self):
        self.target_model.load_state_dict(self.model.state_dict())
    

Policy Gradients with PyTorch (REINFORCE):

import torch
import torch.nn as nn
import torch.optim as optim
import numpy as np

class PolicyNetwork(nn.Module):
    def __init__(self, state_dim, action_dim):
        super(PolicyNetwork, self).__init__()
        self.fc1 = nn.Linear(state_dim, 64)
        self.fc2 = nn.Linear(64, 64)
        self.fc3 = nn.Linear(64, action_dim)

    def forward(self, x):
        x = torch.relu(self.fc1(x))
        x = torch.relu(self.fc2(x))
        x = torch.softmax(self.fc3(x), dim=-1)
        return x

class REINFORCEAgent:
    def __init__(self, state_dim, action_dim):
        self.policy = PolicyNetwork(state_dim, action_dim)
        self.optimizer = optim.Adam(self.policy.parameters(), lr=0.001)
        self.gamma = 0.99
        self.saved_log_probs = []
        self.rewards = []

    def act(self, state):
        state = torch.FloatTensor(state).unsqueeze(0)
        probs = self.policy(state)
        m = torch.distributions.Categorical(probs)
        action = m.sample()
        self.saved_log_probs.append(m.log_prob(action))
        return action.item()

    def update(self):
        R = 0
        policy_loss = []
        returns = []
        for r in self.rewards[::-1]:
            R = r + self.gamma * R
            returns.insert(0, R)
        returns = torch.FloatTensor(returns)
        returns = (returns - returns.mean()) / (returns.std() + 1e-9)  # Normalize
        for log_prob, R in zip(self.saved_log_probs, returns):
            policy_loss.append(-log_prob * R)
        self.optimizer.zero_grad()
        policy_loss = torch.cat(policy_loss).sum()
        policy_loss.backward()
        self.optimizer.step()
        del self.rewards[:]
        del self.saved_log_probs[:]
    

Actor-Critic with PyTorch (A2C):

import torch
import torch.nn as nn
import torch.optim as optim
import numpy as np

class ActorCritic(nn.Module):
    def __init__(self, state_dim, action_dim):
        super(ActorCritic, self).__init__()
        self.fc1 = nn.Linear(state_dim, 64)
        self.fc2 = nn.Linear(64, 64)

        # Actor head
        self.actor = nn.Linear(64, action_dim)

        # Critic head
        self.critic = nn.Linear(64, 1)

    def forward(self, x):
        x = torch.relu(self.fc1(x))
        x = torch.relu(self.fc2(x))
        action_probs = torch.softmax(self.actor(x), dim=-1)
        state_value = self.critic(x)
        return action_probs, state_value

class A2CAgent:
    def __init__(self, state_dim, action_dim):
        self.model = ActorCritic(state_dim, action_dim)
        self.optimizer = optim.Adam(self.model.parameters(), lr=0.001)
        self.gamma = 0.99

    def act(self, state):
        state = torch.FloatTensor(state).unsqueeze(0)
        probs, state_value = self.model(state)
        m = torch.distributions.Categorical(probs)
        action = m.sample()
        log_prob = m.log_prob(action)
        return action.item(), log_prob, state_value

    def update(self, log_probs, state_values, rewards):
        R = 0
        policy_loss = []
        value_loss = []
        returns = []
        for r in rewards[::-1]:
            R = r + self.gamma * R
            returns.insert(0, R)
        returns = torch.FloatTensor(returns)
        returns = (returns - returns.mean()) / (returns.std() + 1e-9)  # Normalize

        for log_prob, value, R in zip(log_probs, state_values, returns):
            advantage = R - value.item()
            policy_loss.append(-log_prob * advantage)
            value_loss.append(nn.MSELoss()(value, torch.FloatTensor([R])))

        self.optimizer.zero_grad()
        loss = torch.stack(policy_loss).sum() + torch.stack(value_loss).sum()
        loss.backward()
        self.optimizer.step()
    

Derivation of the Bellman Equation for \(V^\pi(s)\)

Starting from the definition of the value function:

We can split the sum into the immediate reward and the future rewards:

Using the linearity of expectation and the Markov property:

Recognizing that the expectation inside is the value function at \(s'\):

Value Iteration Algorithm

1. Initialize \(V(s)\) arbitrarily (e.g., \(V(s) = 0\) for all \(s \in S\)).
2. Repeat until convergence:
   1. For each state \(s \in S\), update: \[ V(s) \leftarrow \max_a \sum_{s'} P(s'|s,a) \left[ R(s,a,s') + \gamma V(s') \right] \]
3. Derive the optimal policy: \[ \pi^*(s) = \arg\max_a \sum_{s'} P(s'|s,a) \left[ R(s,a,s') + \gamma V^*(s') \right] \]

Worked Example: Value Iteration on a Simple MDP

Consider an MDP with two states \(S = \{s_1, s_2\}\), one action \(A = \{a\}\), and the following transition and reward:

• \(P(s_1|s_1,a) = 0.5\), \(P(s_2|s_1,a) = 0.5\), \(R(s_1,a,s_1) = 0\), \(R(s_1,a,s_2) = 1\)
• \(P(s_1|s_2,a) = 0\), \(P(s_2|s_2,a) = 1\), \(R(s_2,a,s_2) = 2\)

Let \(\gamma = 0.9\). Initialize \(V(s_1) = V(s_2) = 0\).

Iteration 1:

• \(V(s_1) = \max_a [0.5(0 + 0.9 \cdot 0) + 0.5(1 + 0.9 \cdot 0)] = 0.5\)
• \(V(s_2) = \max_a [1.0(2 + 0.9 \cdot 0)] = 2\)

Iteration 2:

• \(V(s_1) = \max_a [0.5(0 + 0.9 \cdot 0.5) + 0.5(1 + 0.9 \cdot 2)] = 1.625\)
• \(V(s_2) = \max_a [1.0(2 + 0.9 \cdot 2)] = 3.8\)

This process continues until \(V(s)\) converges.

Example: ARIMA(1, 1, 1)

The model can be written as:

Expanding the left-hand side:

Rearranging terms:

Example: SARIMA(1, 1, 1)(1, 1, 1)[12]

The model can be written as:

Expanding the left-hand side:

This results in a complex model with both non-seasonal and seasonal terms.

Example: Local Level Model

A simple state space model where the state \( \alpha_t \) represents the level of the series:

Here, \( T_t = 1 \), \( R_t = 1 \), \( Z_t = 1 \), \( Q_t = \sigma_\eta^2 \), and \( H_t = \sigma_\epsilon^2 \).

1. ARIMA:

• Forecasting stock prices or sales data where trends and autocorrelations are present.
• Modeling temperature or other environmental data with clear temporal dependencies.

2. SARIMA:

• Forecasting retail sales with strong seasonal patterns (e.g., holiday sales).
• Modeling electricity demand, which exhibits daily, weekly, and yearly seasonality.

3. State Space Models:

• Tracking the position and velocity of an object (e.g., in robotics or aerospace).
• Econometric modeling, where latent factors (e.g., "business confidence") drive observed data.
• Signal processing, where the goal is to filter noise from a signal.

1. Object Tracking: Kalman filters are widely used in radar and computer vision for tracking the position and velocity of objects (e.g., aircraft, vehicles, or pedestrians). The state vector might include position, velocity, and acceleration, while measurements come from sensors like radar or cameras.

2. Navigation Systems: In GPS and inertial navigation systems, Kalman filters fuse noisy sensor data (e.g., accelerometers, gyroscopes, GPS) to estimate the position, velocity, and orientation of a vehicle or aircraft.

3. Economics and Finance: Kalman filters are used to estimate hidden states in economic models (e.g., the "true" value of a stock price obscured by market noise) or to track time-varying parameters in financial time series.

4. Robotics: In simultaneous localization and mapping (SLAM), Kalman filters estimate the robot's pose and the positions of landmarks in the environment using noisy sensor data.

Consider a car moving in a straight line with constant velocity. The state vector is \(\mathbf{x}_k = [x_k, \dot{x}_k]^T\), where \(x_k\) is the position and \(\dot{x}_k\) is the velocity. The state transition model is:

where \(\Delta t\) is the time step and \(\sigma_w^2\) is the process noise variance. The measurement model is:

where \(\sigma_v^2\) is the measurement noise variance.

Initialization:

Prediction Step:

Update Step: Suppose the measurement at \(k=1\) is \(z_1 = 2\) with \(\sigma_v^2 = 1\).

Derivation of the Forward Algorithm:

The forward variable \( \alpha_t(i) \) represents the probability of observing the partial sequence \( o_1, o_2, ..., o_t \) and being in state \( s_i \) at time \( t \).

1. Initialization:

   At \( t=1 \), the probability of being in state \( s_i \) and observing \( o_1 \) is:

   \[ \alpha_1(i) = P(o_1, q_1 = s_i | \lambda) = P(q_1 = s_i) P(o_1 | q_1 = s_i) = \pi_i b_i(o_1) \]
2. Recursion:

   For \( t > 1 \), the probability of being in state \( s_j \) at time \( t \) and observing \( o_t \) can be computed by summing over all possible previous states \( s_i \):

   \[ \alpha_t(j) = P(o_1, o_2, ..., o_t, q_t = s_j | \lambda) = \sum_{i=1}^N P(o_1, o_2, ..., o_t, q_{t-1} = s_i, q_t = s_j | \lambda) \]

   Using the Markov property and the definition of \( a_{ij} \) and \( b_j(o_t) \):

   \[ \alpha_t(j) = \sum_{i=1}^N \alpha_{t-1}(i) a_{ij} b_j(o_t) \]
3. Termination:

   The probability of the entire observation sequence is the sum of the forward variables at time \( T \):

   \[ P(O | \lambda) = \sum_{i=1}^N \alpha_T(i) \]

Derivation of the Backward Algorithm:

The backward variable \( \beta_t(i) \) represents the probability of observing the partial sequence \( o_{t+1}, o_{t+2}, ..., o_T \) given that the state at time \( t \) is \( s_i \).

1. Initialization:

   At \( t=T \), there are no more observations, so:

   \[ \beta_T(i) = 1 \]
2. Recursion:

   For \( t < T \), the probability can be computed by summing over all possible next states \( s_j \):

   \[ \beta_t(i) = P(o_{t+1}, o_{t+2}, ..., o_T | q_t = s_i, \lambda) = \sum_{j=1}^N P(o_{t+1}, o_{t+2}, ..., o_T, q_{t+1} = s_j | q_t = s_i, \lambda) \]

   Using the Markov property and the definition of \( a_{ij} \) and \( b_j(o_{t+1}) \):

   \[ \beta_t(i) = \sum_{j=1}^N a_{ij} b_j(o_{t+1}) \beta_{t+1}(j) \]

Derivation of the Viterbi Algorithm:

The Viterbi algorithm finds the most likely sequence of hidden states by keeping track of the maximum probability path to each state at each time step.

1. Initialization:

   At \( t=1 \), the probability of the most likely path ending in state \( s_i \) is:

   \[ \delta_1(i) = \pi_i b_i(o_1) \]

   The backpointer \( \psi_1(i) \) is initialized to 0 since there is no previous state.

2. Recursion:

   For \( t > 1 \), the probability of the most likely path ending in state \( s_j \) at time \( t \) is:

   \[ \delta_t(j) = \max_{1 \leq i \leq N} \left[ \delta_{t-1}(i) a_{ij} \right] b_j(o_t) \]

   The backpointer \( \psi_t(j) \) stores the state \( s_i \) that maximized the above probability:

   \[ \psi_t(j) = \arg\max_{1 \leq i \leq N} \left[ \delta_{t-1}(i) a_{ij} \right] \]
3. Termination:

   The probability of the most likely path is the maximum of the \( \delta_T(i) \) values:

   \[ P^* = \max_{1 \leq i \leq N} \delta_T(i) \]

   The final state in the most likely path is:

   \[ q_T^* = \arg\max_{1 \leq i \leq N} \delta_T(i) \]
4. Path Backtracking:

   The most likely path is obtained by backtracking from \( q_T^* \) using the backpointers \( \psi_t \):

   \[ q_t^* = \psi_{t+1}(q_{t+1}^*), \quad t = T-1, T-2, ..., 1 \]

1. Speech Recognition:

HMMs are widely used in speech recognition systems. The hidden states represent phonemes or words, and the observations are acoustic features extracted from the speech signal. The Viterbi algorithm is used to find the most likely sequence of words given the acoustic observations.

2. Part-of-Speech Tagging:

In natural language processing, HMMs can be used to assign part-of-speech tags to words in a sentence. The hidden states are the part-of-speech tags, and the observations are the words in the sentence. The Forward-Backward algorithm can be used to compute the probability of each tag for a given word, and the Viterbi algorithm can find the most likely sequence of tags.

3. Bioinformatics:

HMMs are used in bioinformatics for gene prediction and sequence alignment. For example, in gene prediction, the hidden states represent different regions of a DNA sequence (e.g., exons, introns, intergenic regions), and the observations are the nucleotide sequences. The Viterbi algorithm can be used to find the most likely path through the hidden states, effectively predicting the gene structure.

4. Financial Time Series Analysis:

HMMs can model financial time series data where the hidden states represent different market regimes (e.g., bull market, bear market), and the observations are the financial returns. The Forward-Backward algorithm can be used to compute the probability of being in each regime at any given time, and the Viterbi algorithm can identify the most likely sequence of regimes.

Derivation of the Chain Rule for BNs:

1. Start with the joint probability \( P(X_1, X_2, \dots, X_n) \).
2. Apply the chain rule of probability: \[ P(X_1, X_2, \dots, X_n) = P(X_1) P(X_2 \mid X_1) P(X_3 \mid X_1, X_2) \dots P(X_n \mid X_1, \dots, X_{n-1}) \]
3. By the Markov property of BNs, each variable \( X_i \) is conditionally independent of its non-descendants given its parents \( \text{Pa}(X_i) \). Thus, the conditional probabilities simplify to: \[ P(X_i \mid X_1, \dots, X_{i-1}) = P(X_i \mid \text{Pa}(X_i)) \]
4. Substitute back to get the BN chain rule: \[ P(X_1, X_2, \dots, X_n) = \prod_{i=1}^n P(X_i \mid \text{Pa}(X_i)) \]

Derivation of d-Separation (Example):

Consider the BN: \( A \rightarrow B \rightarrow C \) and \( A \rightarrow D \leftarrow C \). Show that \( A \perp\!\!\!\perp C \mid B \).

1. Identify paths between \( A \) and \( C \):
   • Path 1: \( A \rightarrow B \rightarrow C \) (chain).
   • Path 2: \( A \rightarrow D \leftarrow C \) (collider).
2. For \( A \perp\!\!\!\perp C \mid B \), all paths must be blocked by \( B \):
   • Path 1 is blocked because \( B \) is observed (chain rule).
   • Path 2 is blocked because \( D \) is a collider and neither \( D \) nor its descendants are observed.
3. Thus, \( A \perp\!\!\!\perp C \mid B \).

Medical Diagnosis: BNs are used to model relationships between diseases and symptoms. For example:

• Nodes: Diseases (e.g., "Flu"), symptoms (e.g., "Fever"), and test results.
• Edges: Conditional dependencies (e.g., "Flu" causes "Fever").
• Inference: Compute \( P(\text{Disease} \mid \text{Symptoms}) \) to assist diagnosis.

Spam Filtering: BNs can model the probability of an email being spam based on features like word frequencies or sender reputation. Inference is used to classify emails as spam or not spam.

Genetics: BNs model inheritance patterns and the probability of genetic disorders given family history. For example, computing \( P(\text{Disease} \mid \text{Parental Genotypes}) \).

Robotics: BNs are used for sensor fusion and decision-making under uncertainty. For example, a robot may use a BN to estimate its location given noisy sensor data.

Example: Suppose \( p(x) = \mathcal{N}(x; 0, 1) \) and \( q(x) = \mathcal{N}(x; 1, 1) \). We want to estimate \( \mathbb{E}_{p}[x^2] \).

1. Sample \( x_i \sim q(x) \), i.e., \( x_i \sim \mathcal{N}(1, 1) \).
2. Compute the importance weights \( w(x_i) = \frac{p(x_i)}{q(x_i)} = \exp\left( - \frac{1}{2} (x_i^2 - (x_i - 1)^2) \right) \).
3. Compute the estimator: \( \hat{\mathbb{E}}_{p}[x^2] = \frac{1}{N} \sum_{i=1}^{N} x_i^2 w(x_i) \).

Example: Sampling from a Gaussian distribution \( \pi(x) = \mathcal{N}(x; 0, 1) \) using a symmetric proposal distribution \( q(x' | x) = \mathcal{N}(x'; x, \sigma^2) \).

1. Initialize \( x_0 \).
2. For each iteration:
   1. Propose \( x' \sim \mathcal{N}(x_t, \sigma^2) \).
   2. Compute \( \alpha = \min\left(1, \frac{\pi(x')}{\pi(x_t)}\right) = \min\left(1, \exp\left( -\frac{1}{2} (x'^2 - x_t^2) \right)\right) \).
   3. Accept \( x' \) with probability \( \alpha \).

Bayesian Inference: MCMC is widely used in Bayesian statistics to sample from posterior distributions, especially when the posterior is not analytically tractable. For example, in hierarchical models or complex likelihoods.

Reinforcement Learning: Monte Carlo methods are used in reinforcement learning for policy evaluation, where the goal is to estimate the value function of a given policy by averaging sampled returns.

Computer Graphics: Monte Carlo integration is used in rendering algorithms to compute global illumination by simulating the transport of light.

Physics: Monte Carlo methods are used to simulate systems with a large number of coupled degrees of freedom, such as in statistical mechanics or quantum chromodynamics.

Finance: Importance sampling is used to price complex financial derivatives and to estimate risk measures like Value at Risk (VaR).

Derivation: Tail Dependence for the Clayton Copula

For the bivariate Clayton copula \(C_{\theta}(u, v) = (u^{-\theta} + v^{-\theta} - 1)^{-1/\theta}\), the lower tail dependence coefficient is derived as follows:

1. Compute \(C(u, u)\): \[ C(u, u) = (2u^{-\theta} - 1)^{-1/\theta}. \]
2. Compute the limit: \[ \lambda_L = \lim_{u \to 0^+} \frac{C(u, u)}{u} = \lim_{u \to 0^+} \frac{(2u^{-\theta} - 1)^{-1/\theta}}{u}. \]
3. Simplify the expression: \[ \lambda_L = \lim_{u \to 0^+} \left( 2 - u^{\theta} \right)^{-1/\theta} = 2^{-1/\theta}. \]

Practical Application: Risk Modeling in Finance

Copulas are widely used in finance to model dependencies between asset returns, especially in risk management and portfolio optimization. For example:

• Value-at-Risk (VaR): Copulas help model the joint distribution of asset returns to estimate the VaR of a portfolio, accounting for tail dependencies.
• Credit Risk: The Clayton copula is often used to model default dependencies due to its lower-tail dependence, capturing the likelihood of joint defaults during market downturns.
• Insurance: The Gumbel copula is used to model extreme events (e.g., natural disasters) due to its upper-tail dependence.

Example: Suppose we model the joint distribution of two stock returns using a Clayton copula with \(\theta = 2\). The lower tail dependence is:

This indicates a 70.7% probability that one stock will experience a large loss given that the other does, highlighting the importance of tail dependence in risk assessment.

Python Example: Fitting a Clayton Copula

import numpy as np
from copulae import ClaytonCopula

# Generate synthetic data with lower-tail dependence
np.random.seed(42)
n = 1000
theta_true = 2.0
cop = ClaytonCopula(theta=theta_true, dim=2)
data = cop.random(n)

# Fit the Clayton copula
clayton = ClaytonCopula(dim=2)
clayton.fit(data)
print(f"Estimated theta: {clayton.params[0]:.3f}")  # Should be close to 2.0

This example generates synthetic data from a Clayton copula with \(\theta = 2\) and fits the copula to the data, recovering the parameter.

Example: Consider the following survival data (time in months, event indicator: 1 = event, 0 = censored):

Compute the Kaplan-Meier estimate at each event time:

1. At \( t = 2 \): \( d_1 = 1 \), \( n_1 = 5 \). \[ \hat{S}(2) = 1 - \frac{1}{5} = 0.8 \]
2. At \( t = 5 \): \( d_2 = 1 \), \( n_2 = 3 \) (subject at \( t=3 \) is censored, so not at risk at \( t=5 \)). \[ \hat{S}(5) = 0.8 \times \left(1 - \frac{1}{3}\right) = 0.8 \times \frac{2}{3} \approx 0.533 \]
3. At \( t = 8 \): \( d_3 = 1 \), \( n_3 = 2 \) (subject at \( t=10 \) is still at risk). \[ \hat{S}(8) = 0.533 \times \left(1 - \frac{1}{2}\right) = 0.533 \times 0.5 = 0.267 \]

The final Kaplan-Meier curve is a step function with values 1 (at \( t=0 \)), 0.8 (at \( t=2 \)), 0.533 (at \( t=5 \)), and 0.267 (at \( t=8 \)).

Example: Suppose we have the following data for 3 subjects:

Compute the partial likelihood for \( \boldsymbol{\beta} = (\beta_1, \beta_2) \):

1. At \( t = 2 \): Risk set \( R(2) = \{1, 2, 3\} \). \[ \text{Numerator} = \exp(\beta_1 \cdot 1 + \beta_2 \cdot 0) = e^{\beta_1} \] \[ \text{Denominator} = e^{\beta_1} + e^{\beta_2} + e^{\beta_1 + \beta_2} \]
2. At \( t = 5 \): Risk set \( R(5) = \{2, 3\} \) (subject 1 has already experienced the event). \[ \text{Numerator} = \exp(\beta_1 \cdot 1 + \beta_2 \cdot 1) = e^{\beta_1 + \beta_2} \] \[ \text{Denominator} = e^{\beta_2} + e^{\beta_1 + \beta_2} \]
3. Partial likelihood: \[ L(\beta_1, \beta_2) = \frac{e^{\beta_1}}{e^{\beta_1} + e^{\beta_2} + e^{\beta_1 + \beta_2}} \times \frac{e^{\beta_1 + \beta_2}}{e^{\beta_2} + e^{\beta_1 + \beta_2}} \]

The log-partial likelihood is maximized numerically to estimate \( \beta_1 \) and \( \beta_2 \).

Applications of Survival Analysis:

• Medical Research: Analyzing time until death, relapse, or recovery (e.g., clinical trials for cancer treatments).
• Engineering: Modeling time until failure of mechanical components (reliability analysis).
• Economics: Studying duration of unemployment or time until loan default.
• Social Sciences: Analyzing time until marriage, divorce, or recidivism.
• Customer Analytics: Predicting churn (time until a customer stops using a service).

Kaplan-Meier Estimator with lifelines:

from lifelines import KaplanMeierFitter
import pandas as pd

# Example data
data = pd.DataFrame({
    'time': [2, 3, 5, 8, 10],
    'event': [1, 0, 1, 1, 0]
})

# Fit Kaplan-Meier estimator
kmf = KaplanMeierFitter()
kmf.fit(data['time'], event_observed=data['event'])

# Plot survival function
kmf.plot_survival_function()
plt.title('Kaplan-Meier Survival Curve')
plt.show()

Cox Proportional Hazards Model with lifelines:

from lifelines import CoxPHFitter

# Example data
data = pd.DataFrame({
    'time': [2, 3, 5, 8, 10],
    'event': [1, 0, 1, 1, 0],
    'age': [50, 60, 45, 55, 65],
    'treatment': [1, 0, 1, 0, 1]
})

# Fit Cox model
cph = CoxPHFitter()
cph.fit(data, duration_col='time', event_col='event', formula='age + treatment')

# Print summary
cph.print_summary()

# Plot coefficients
cph.plot()

Cox Model with scikit-survival:

from sksurv.linear_model import CoxPHSurvivalAnalysis
from sksurv.datasets import load_whas500
from sklearn.model_selection import train_test_split

# Load example data
X, y = load_whas500()

# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Fit Cox model
model = CoxPHSurvivalAnalysis()
model.fit(X_train, y_train)

# Evaluate
print("Concordance index:", model.score(X_test, y_test))

Example: Tuning the hyperparameters \(C\) (regularization) and \(\gamma\) (kernel coefficient) for an SVM with:

• \(C \in \{0.1, 1, 10\}\)
• \(\gamma \in \{0.01, 0.1, 1\}\)

The grid search evaluates all \(3 \times 3 = 9\) combinations.

Example: Using the same SVM hyperparameters as above, random search might sample the following combinations (assuming \(k = 5\)):

• (\(C = 0.1\), \(\gamma = 0.1\))
• (\(C = 10\), \(\gamma = 0.01\))
• (\(C = 1\), \(\gamma = 1\))
• (\(C = 0.1\), \(\gamma = 1\))
• (\(C = 10\), \(\gamma = 1\))

Example: Bayesian optimization for tuning the learning rate \(\eta\) and number of layers \(L\) of a neural network:

1. Initialize with a few random evaluations of \(f(\eta, L)\).
2. Fit a GP surrogate model to the observed data.
3. Use EI to select the next \((\eta, L)\) to evaluate.
4. Evaluate \(f(\eta, L)\) and update the GP model.
5. Repeat until convergence or a budget is exhausted.

1. Grid Search in scikit-learn:

from sklearn.model_selection import GridSearchCV
from sklearn.svm import SVC

param_grid = {
    'C': [0.1, 1, 10],
    'gamma': [0.01, 0.1, 1],
    'kernel': ['rbf', 'linear']
}

svm = SVC()
grid_search = GridSearchCV(svm, param_grid, cv=5)
grid_search.fit(X_train, y_train)

print("Best parameters:", grid_search.best_params_)
print("Best score:", grid_search.best_score_)

2. Random Search in scikit-learn:

from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import loguniform

param_dist = {
    'C': loguniform(1e-2, 1e2),
    'gamma': loguniform(1e-3, 1e1),
    'kernel': ['rbf', 'linear']
}

random_search = RandomizedSearchCV(svm, param_dist, n_iter=20, cv=5)
random_search.fit(X_train, y_train)

print("Best parameters:", random_search.best_params_)
print("Best score:", random_search.best_score_)

3. Bayesian Optimization with Optuna:

import optuna
from sklearn.svm import SVC
from sklearn.model_selection import cross_val_score

def objective(trial):
    C = trial.suggest_float('C', 1e-2, 1e2, log=True)
    gamma = trial.suggest_float('gamma', 1e-3, 1e1, log=True)
    kernel = trial.suggest_categorical('kernel', ['rbf', 'linear'])

    svm = SVC(C=C, gamma=gamma, kernel=kernel)
    score = cross_val_score(svm, X_train, y_train, cv=5).mean()
    return score

study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=50)

print("Best parameters:", study.best_params)
print("Best score:", study.best_value)

Solgi's PyPI Genetic Algorithm (geneticalgorithm) with scikit-learn:

# Install Ryan (Mohammad) Solgi's package:
# pip install geneticalgorithm

import numpy as np
from geneticalgorithm import geneticalgorithm as ga
from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestClassifier

# Objective must be minimized in this package, so we return negative CV accuracy.
def objective(x):
    n_estimators = int(x[0])               # [50, 500]
    max_depth = int(x[1])                  # [1, 30]
    min_samples_split = int(x[2])          # [2, 20]

    model = RandomForestClassifier(
        n_estimators=n_estimators,
        max_depth=max_depth,
        min_samples_split=min_samples_split,
        random_state=42
    )
    score = cross_val_score(model, X_train, y_train, cv=5, scoring="accuracy").mean()
    return -score

varbound = np.array([[50, 500], [1, 30], [2, 20]])

algorithm_param = {
    'max_num_iteration': 80,
    'population_size': 60,
    'mutation_probability': 0.1,
    'elit_ratio': 0.05,      # key elitist parameter
    'crossover_probability': 0.8,
    'parents_portion': 0.3,
    'crossover_type': 'uniform',
    'max_iteration_without_improv': 15
}

model = ga(
    function=objective,
    dimension=3,
    variable_type='int',
    variable_boundaries=varbound,
    algorithm_parameters=algorithm_param
)
model.run()

best_solution = model.output_dict['variable']
best_cv_acc = -model.output_dict['function']
print(best_solution, best_cv_acc)

Derivation: Why k-Fold CV Reduces Variance

The variance of the k-fold CV estimate can be derived as follows:

1. Assume each fold's score \(\text{Score}_i\) is an independent and identically distributed (i.i.d.) random variable with variance \(\sigma^2\).
2. The average CV score is \(\text{CV Score} = \frac{1}{k} \sum_{i=1}^{k} \text{Score}_i\).
3. The variance of the average is: \[ \text{Var}(\text{CV Score}) = \text{Var}\left(\frac{1}{k} \sum_{i=1}^{k} \text{Score}_i\right) = \frac{1}{k^2} \sum_{i=1}^{k} \text{Var}(\text{Score}_i) = \frac{1}{k^2} \cdot k \sigma^2 = \frac{\sigma^2}{k}. \]

Thus, increasing k reduces the variance of the CV estimate.

Step-by-Step: Stratified k-Fold in Practice

Given a dataset with classes \(C = \{c_1, c_2, ..., c_m\}\), stratified k-fold ensures:

1. Calculate the proportion of each class in the full dataset: \[ p_c = \frac{\text{Count}(y = c)}{N}, \quad \text{where } N \text{ is the total number of samples.} \]
2. For each fold \(i\), ensure the proportion of class \(c\) in the fold is \(p_c\).
3. Randomly sample (without replacement) from each class to construct the folds.

This preserves the class distribution in every fold, reducing bias in performance estimates for imbalanced datasets.

Step-by-Step: Time Series CV (Rolling Window)

For a time series dataset with \(T\) observations:

1. Define the forecast horizon \(h\) (e.g., predict the next 5 time steps).
2. Define the step size \(s\) (e.g., move the window forward by 5 time steps).
3. For each fold \(i\) (from 1 to \(k\)):
   • Training set: First \(T - h - (k - i) \cdot s\) observations.
   • Validation set: Next \(h\) observations after the training set.
4. Slide the window forward by \(s\) time steps for the next fold.

This ensures the temporal order is preserved, and the model is evaluated on "future" data.

When to Use k-Fold CV:

• Small to medium-sized datasets where maximizing data usage is critical.
• Datasets with no temporal dependencies or class imbalance.
• Hyperparameter tuning (e.g., using GridSearchCV in scikit-learn).

Example (scikit-learn):

from sklearn.model_selection import KFold, cross_val_score
from sklearn.ensemble import RandomForestClassifier

X, y = load_data()
kfold = KFold(n_splits=5, shuffle=True, random_state=42)
model = RandomForestClassifier()
scores = cross_val_score(model, X, y, cv=kfold, scoring='accuracy')
print(f"Mean CV Accuracy: {scores.mean():.3f}")

When to Use Stratified k-Fold CV:

• Imbalanced datasets (e.g., fraud detection, rare disease classification).
• Multi-class classification problems where class distribution matters.

Example (scikit-learn):

from sklearn.model_selection import StratifiedKFold

skfold = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
scores = cross_val_score(model, X, y, cv=skfold, scoring='f1_macro')
print(f"Mean CV F1-Score: {scores.mean():.3f}")

When to Use Time Series CV:

• Forecasting tasks (e.g., stock prices, weather prediction).
• Any dataset where observations are temporally ordered.

Example (scikit-learn):

from sklearn.model_selection import TimeSeriesSplit

tscv = TimeSeriesSplit(n_splits=5)
for train_index, val_index in tscv.split(X):
    X_train, X_val = X[train_index], X[val_index]
    y_train, y_val = y[train_index], y[val_index]
    model.fit(X_train, y_train)
    score = model.score(X_val, y_val)
    print(f"Fold Score: {score:.3f}")

Example: Lasso in scikit-learn

from sklearn.linear_model import Lasso
from sklearn.datasets import make_regression

# Generate synthetic data
X, y = make_regression(n_samples=100, n_features=20, n_informative=5, noise=0.5)

# Fit Lasso model
lasso = Lasso(alpha=0.1)
lasso.fit(X, y)

# Selected features (non-zero coefficients)
selected_features = [i for i, coef in enumerate(lasso.coef_) if coef != 0]
print("Selected features:", selected_features)

Example: Mutual Information in scikit-learn

from sklearn.feature_selection import mutual_info_classif, mutual_info_regression
from sklearn.datasets import load_iris

# Classification example
X, y = load_iris(return_X_y=True)
mi_scores = mutual_info_classif(X, y)
print("Mutual Information Scores (Classification):", mi_scores)

# Regression example
X, y = make_regression(n_samples=100, n_features=10, n_informative=3, noise=0.5)
mi_scores = mutual_info_regression(X, y)
print("Mutual Information Scores (Regression):", mi_scores)

Example: RFE in scikit-learn

from sklearn.feature_selection import RFE
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import load_breast_cancer

# Load data
X, y = load_breast_cancer(return_X_y=True)

# Create a model (Logistic Regression)
model = LogisticRegression(max_iter=1000)

# Create RFE object
rfe = RFE(estimator=model, n_features_to_select=5)

# Fit RFE
rfe.fit(X, y)

# Selected features
selected_features = [i for i, selected in enumerate(rfe.support_) if selected]
print("Selected features:", selected_features)

# Feature rankings (1 = selected, higher = eliminated earlier)
print("Feature rankings:", rfe.ranking_)

1. High-Dimensional Data (e.g., Genomics, Text Data):

• Lasso is widely used in genomics to identify a small subset of genes associated with a disease.
• Mutual information is used in text classification to select the most informative words.

2. Model Interpretability:

• Lasso and RFE produce sparse models, making them easier to interpret.
• Mutual information can identify non-linear relationships that are not captured by linear models.

3. Preprocessing for Other Models:

• Feature selection can improve the performance of models sensitive to irrelevant features (e.g., k-NN, SVM).
• Reducing the feature space can speed up training for computationally expensive models (e.g., deep learning).

Example: Consider a 2D minority class sample \( x_i = [1, 2] \) and its nearest neighbor \( x_{zi} = [3, 4] \). If \( \lambda = 0.5 \), the synthetic sample is:

Example: For a dataset with 1000 samples where 900 belong to the majority class and 100 to the minority class:

The minority class is given a weight of 5, making misclassifications of minority samples 5 times more costly.

Example: For a dataset with \( n = 100 \), the normalization factor \( c(100) \) is approximately 5.187. If the average path length \( E(h(x)) \) for a sample is 3, its anomaly score is:

Scores close to 1 indicate anomalies, while scores close to 0 indicate normal instances.

1. Fraud Detection: In credit card transactions, fraudulent transactions are rare (minority class). Techniques like SMOTE or class weighting can improve the detection of fraudulent transactions.

2. Medical Diagnosis: Diseases like cancer are rare in the general population. Imbalanced learning techniques can help in building models that accurately predict the presence of such diseases.

3. Manufacturing Defect Detection: Anomaly detection algorithms can identify defective products on an assembly line, where defects are rare but critical to detect.

Example: Building a Computational Graph in PyTorch

import torch

# Create tensors with gradient tracking
x = torch.tensor(2.0, requires_grad=True)
w = torch.tensor(3.0, requires_grad=True)
b = torch.tensor(1.0, requires_grad=True)

# Forward pass: y = w * x + b
y = w * x + b

# Backward pass: compute gradients
y.backward()

# Gradients are now available
print(f"dy/dx = {x.grad}")  # 3.0
print(f"dy/dw = {w.grad}")  # 2.0
print(f"dy/db = {b.grad}")  # 1.0
    

Explanation:

1. PyTorch builds a computational graph during the forward pass.
2. When y.backward() is called, PyTorch traverses the graph backward using the chain rule.
3. Gradients are accumulated in the .grad attribute of leaf tensors.

Example: Multi-Layer Perceptron (MLP) Gradient Flow

Consider a simple MLP with one hidden layer:

Where \(\sigma\) is the ReLU activation. The gradient of the loss \(L\) with respect to \(W_1\) is:

Expanding each term:

PyTorch's autograd handles this computation automatically.

Example: Jacobian Computation in PyTorch

import torch

def f(x):
    return torch.stack([x[0] ** 2, x[1] * x[0]])

x = torch.tensor([2.0, 3.0], requires_grad=True)
y = f(x)

# Compute Jacobian
jacobian = []
for i in range(y.shape[0]):
    grad_output = torch.zeros_like(y)
    grad_output[i] = 1.0
    gradients = torch.autograd.grad(y, x, grad_outputs=grad_output, retain_graph=True)
    jacobian.append(gradients[0])

jacobian = torch.stack(jacobian)
print("Jacobian:")
print(jacobian)
    

Output:

Jacobian:
tensor([[4., 0.],
        [3., 2.]])
    

This matches the analytical Jacobian:

Example: Custom Autograd Function

import torch

class Exp(torch.autograd.Function):
    @staticmethod
    def forward(ctx, x):
        ctx.save_for_backward(x)
        return x.exp()

    @staticmethod
    def backward(ctx, grad_output):
        x, = ctx.saved_tensors
        return grad_output * x.exp()

# Usage
x = torch.tensor(1.0, requires_grad=True)
y = Exp.apply(x)
y.backward()
print(f"dy/dx = {x.grad}")  # e^1 = 2.718...
    

Example: Custom Transformer for Log Scaling

from sklearn.base import BaseEstimator, TransformerMixin
import numpy as np

class LogScaler(BaseEstimator, TransformerMixin):
    def __init__(self, add_epsilon=True):
        self.add_epsilon = add_epsilon

    def fit(self, X, y=None):
        # No parameters to learn; return self
        return self

    def transform(self, X):
        if self.add_epsilon:
            X = X + 1e-6  # Avoid log(0)
        return np.log(X)
    

This transformer applies a log transformation to the input data, optionally adding a small epsilon to avoid numerical issues.

Example: Applying Different Transformers to Numerical and Categorical Features

from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.impute import SimpleImputer

# Define transformers for numerical and categorical features
numerical_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler', StandardScaler())
])

categorical_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='most_frequent')),
    ('onehot', OneHotEncoder(handle_unknown='ignore'))
])

# Apply transformers to specific columns
preprocessor = ColumnTransformer(
    transformers=[
        ('num', numerical_transformer, ['age', 'income']),
        ('cat', categorical_transformer, ['gender', 'country'])
    ],
    remainder='drop'  # Drop columns not specified
)
    

Application 1: End-to-End Machine Learning Workflow

Pipelines are essential for deploying machine learning models, as they encapsulate the entire preprocessing and modeling workflow. For example:

from sklearn.ensemble import RandomForestClassifier
from sklearn.pipeline import Pipeline

# Define the full pipeline
model = Pipeline(steps=[
    ('preprocessor', preprocessor),  # ColumnTransformer from earlier
    ('classifier', RandomForestClassifier())
])

# Fit and predict
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
    

This ensures that the same preprocessing steps are applied during training and prediction, avoiding data leakage.

Application 2: Hyperparameter Tuning with Pipelines

Pipelines can be used with GridSearchCV or RandomizedSearchCV to tune hyperparameters for both preprocessing steps and the estimator. For example:

from sklearn.model_selection import GridSearchCV

param_grid = {
    'preprocessor__num__imputer__strategy': ['mean', 'median'],
    'classifier__n_estimators': [50, 100, 200]
}

grid_search = GridSearchCV(model, param_grid, cv=5)
grid_search.fit(X_train, y_train)
    

Example: Consider a model with 3 features \(F = \{1, 2, 3\}\). The SHAP value for feature 1 is computed as:

This averages the marginal contribution of feature 1 across all possible feature subsets.

Example: For a linear interpretable model \(g(z) = w_0 + \sum_{i=1}^M w_i z_i\), LIME minimizes:

where \(Z\) is the perturbed dataset, and \(\pi_x(z) = \exp(-D(x, z)^2 / \sigma^2)\) is an exponential kernel with distance \(D\) and bandwidth \(\sigma\).

Example: For a dataset with 1000 samples and a model \(f\), the PDP for feature \(j\) at value \(x_j = 0.5\) is computed as:

This averages the model's predictions when feature \(j\) is set to 0.5 for all samples.

1. Debugging Models:

• SHAP/LIME can identify if a model is relying on spurious correlations (e.g., a hospital's zip code instead of patient symptoms).
• PDPs can reveal non-monotonic relationships (e.g., a drug's efficacy increasing then decreasing with dosage).

2. Regulatory Compliance:

• SHAP values can provide "right to explanation" under GDPR by quantifying feature contributions.
• LIME explanations can be presented to non-technical stakeholders (e.g., "The model denied your loan because of your credit score and debt-to-income ratio").

3. Feature Engineering:

• PDPs can guide feature transformations (e.g., log-transforming a feature with a non-linear relationship).
• SHAP can identify redundant features (e.g., two highly correlated features with low SHAP values).

4. Fairness Audits:

• SHAP can detect bias by comparing feature contributions across demographic groups.
• LIME can explain individual cases of discrimination (e.g., "The model gave a lower score to this applicant because of their gender").

1. SHAP with scikit-learn (TreeSHAP for Random Forest):

import shap
from sklearn.ensemble import RandomForestClassifier

# Train a model
model = RandomForestClassifier().fit(X_train, y_train)

# Explain the model's predictions
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# Visualize the first prediction's explanation
shap.force_plot(explainer.expected_value[0], shap_values[0][0,:], X_test.iloc[0,:])

# Summary plot of global feature importance
shap.summary_plot(shap_values, X_test)

2. LIME with scikit-learn:

import lime
import lime.lime_tabular
from sklearn.ensemble import RandomForestClassifier

# Train a model
model = RandomForestClassifier().fit(X_train, y_train)

# Explain an instance
explainer = lime.lime_tabular.LimeTabularExplainer(
    X_train.values,
    feature_names=X_train.columns,
    class_names=['class_0', 'class_1'],
    mode='classification'
)
exp = explainer.explain_instance(
    X_test.iloc[0].values,
    model.predict_proba,
    num_features=10
)

# Show explanation
exp.show_in_notebook()

3. Partial Dependence Plots with scikit-learn:

from sklearn.inspection import PartialDependenceDisplay
from sklearn.ensemble import GradientBoostingRegressor

# Train a model
model = GradientBoostingRegressor().fit(X_train, y_train)

# Plot PDP for feature 0
PartialDependenceDisplay.from_estimator(
    model,
    X_train,
    features=[0],
    feature_names=X_train.columns
)

# Plot 2D PDP for features 0 and 1
PartialDependenceDisplay.from_estimator(
    model,
    X_train,
    features=[(0, 1)],
    feature_names=X_train.columns
)

4. SHAP with PyTorch (DeepSHAP for Neural Networks):

import shap
import torch
import torch.nn as nn

# Define a simple neural network
class Net(nn.Module):
    def __init__(self):
        super(Net, self).__init__()
        self.fc1 = nn.Linear(10, 5)
        self.fc2 = nn.Linear(5, 1)

    def forward(self, x):
        x = torch.relu(self.fc1(x))
        x = self.fc2(x)
        return x

model = Net()
model.load_state_dict(torch.load('model.pth'))
model.eval()

# Create a SHAP explainer
background = torch.randn(100, 10)  # Background dataset
explainer = shap.DeepExplainer(model, background)
shap_values = explainer.shap_values(torch.tensor(X_test.values, dtype=torch.float32))

# Plot summary
shap.summary_plot(shap_values, X_test)

Interpretation: Delta tells you how the portfolio reacts to an infinitesimal move right now. A deep hedging policy decides what trade to make now while accounting for future costs, future risk, remaining time, current inventory, and the fact that you may rebalance again later.

Honest summary:

• In a Black-Scholes style frictionless continuous-time world, AI adds little conceptual value because the classical solution is already known.
• In realistic desk settings with costs, discrete rebalancing, multiple instruments, and constraints, AI can matter a lot because the problem becomes a difficult dynamic optimization problem.