xbe.at

Deep learning / Backpropagation and optimization

Overcoming Vanishing Gradients in Deep Neural Networks

10 min read, 14 sections

Contents
  1. Understanding Vanishing Gradients Mathematically
  2. Implementing Gradient Flow Visualization
  3. ReLU Implementation from Scratch
  4. Batch Normalization Implementation
  5. Residual Network Building Block
  6. Xavier/Glorot Weight Initialization
  7. Gradient Clipping Implementation
  8. LSTM Cell for Long-Term Dependencies
  9. Practical Example - Time Series Prediction with Gradient Monitoring
  10. Real-world Application - Stock Price Prediction with Gradient Monitoring
  11. Results Analysis for Stock Price Prediction
  12. GRU Implementation for Gradient Flow Comparison
  13. Comparative Analysis: LSTM vs GRU vs Vanilla RNN
  14. Additional Resources

Understanding Vanishing Gradients Mathematically

The vanishing gradient problem can be understood through the chain rule of calculus in neural networks. As gradients flow backward through many layers, repeated multiplication of small values causes exponential decrease, particularly with sigmoid and tanh activation functions.

import numpy as np

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def sigmoid_derivative(x):
    s = sigmoid(x)
    return s * (1 - s)

# Demonstrate gradient decay
x = np.linspace(-5, 5, 100)
derivatives = sigmoid_derivative(x)
print(f"Max derivative value: {np.max(derivatives):.4f}")  # Always < 0.25

Implementing Gradient Flow Visualization

A practical implementation to visualize how gradients diminish through layers using a simple feed-forward network. This helps understand why deeper networks struggle with learning in early layers.

import numpy as np
import matplotlib.pyplot as plt

class VanishingGradientDemo:
    def __init__(self, n_layers=10):
        self.weights = [np.random.randn(1, 1) * 0.1 for _ in range(n_layers)]
        
    def forward(self, x):
        activations = []
        for w in self.weights:
            x = np.tanh(np.dot(x, w))
            activations.append(x)
        return activations

    def backward(self, activations):
        gradients = []
        grad = 1.0
        for a in reversed(activations):
            grad *= (1 - a**2)  # tanh derivative
            gradients.append(abs(grad))
        return gradients[::-1]

# Demonstration
demo = VanishingGradientDemo(10)
x = np.array([[1.0]])
activations = demo.forward(x)
gradients = demo.backward(activations)

print("Gradient magnitudes through layers:")
for i, g in enumerate(gradients):
    print(f"Layer {i+1}: {g[0][0]:.10f}")

ReLU Implementation from Scratch

The Rectified Linear Unit (ReLU) activation function helps mitigate vanishing gradients by maintaining constant gradients for positive inputs, unlike sigmoid or tanh functions which compress gradients.

import numpy as np

class ReLU:
    def forward(self, x):
        self.input = x
        return np.maximum(0, x)
    
    def backward(self, grad_output):
        grad_input = grad_output.copy()
        grad_input[self.input <= 0] = 0
        return grad_input

# Demonstration
relu = ReLU()
x = np.array([-2, -1, 0, 1, 2])
forward_output = relu.forward(x)
backward_output = relu.backward(np.ones_like(x))

print(f"Input: {x}")
print(f"Forward: {forward_output}")
print(f"Gradient: {backward_output}")

Batch Normalization Implementation

class BatchNormalization:
    def __init__(self, epsilon=1e-5):
        self.epsilon = epsilon
        self.gamma = 1.0
        self.beta = 0.0
        
    def forward(self, x):
        self.input = x
        self.mean = np.mean(x, axis=0)
        self.var = np.var(x, axis=0)
        
        # Normalize
        self.x_norm = (x - self.mean) / np.sqrt(self.var + self.epsilon)
        
        # Scale and shift
        out = self.gamma * self.x_norm + self.beta
        return out
    
    def backward(self, grad_output):
        # Simplified backward pass
        N = float(self.input.shape[0])
        grad_input = self.gamma * grad_output / np.sqrt(self.var + self.epsilon)
        return grad_input

# Example usage
bn = BatchNormalization()
x = np.random.randn(32, 10)  # Batch of 32 samples, 10 features
normalized = bn.forward(x)
print(f"Input mean: {np.mean(x):.4f}, var: {np.var(x):.4f}")
print(f"Output mean: {np.mean(normalized):.4f}, var: {np.var(normalized):.4f}")

Residual Network Building Block

A residual connection provides a direct path for gradients to flow backward, preventing vanishing gradients in very deep networks. This implementation shows a basic ResNet block with skip connection.

import numpy as np

class ResNetBlock:
    def __init__(self, input_channels, output_channels):
        self.weights1 = np.random.randn(input_channels, output_channels) * 0.01
        self.weights2 = np.random.randn(output_channels, output_channels) * 0.01
        
    def forward(self, x):
        self.input = x
        # First convolution + ReLU
        self.layer1_out = np.maximum(0, np.dot(x, self.weights1))
        # Second convolution
        self.layer2_out = np.dot(self.layer1_out, self.weights2)
        # Skip connection
        return self.layer2_out + x
    
    def backward(self, grad_output):
        # Gradient flows through both paths
        grad_skip = grad_output
        grad_main = grad_output
        return grad_skip + grad_main

# Example usage
block = ResNetBlock(64, 64)
x = np.random.randn(32, 64)  # Batch of 32 samples, 64 features
output = block.forward(x)
print(f"Input shape: {x.shape}")
print(f"Output shape: {output.shape}")

Xavier/Glorot Weight Initialization

Proper weight initialization is crucial for preventing vanishing gradients. Xavier initialization scales weights based on the number of input and output connections to maintain consistent gradient flow.

import numpy as np

def xavier_init(shape, gain=1.0):
    fan_in = shape[0] if len(shape) > 1 else shape[0]
    fan_out = shape[1] if len(shape) > 1 else shape[0]
    
    # Calculate standard deviation based on Xavier formula
    std = gain * np.sqrt(2.0 / (fan_in + fan_out))
    
    # Generate weights from normal distribution
    weights = np.random.normal(0, std, shape)
    return weights

# Compare gradient flow with different initializations
def compare_gradients(input_size, hidden_size, n_layers):
    # Standard initialization
    std_weights = [np.random.randn(hidden_size, hidden_size) for _ in range(n_layers)]
    
    # Xavier initialization
    xavier_weights = [xavier_init((hidden_size, hidden_size)) for _ in range(n_layers)]
    
    x = np.random.randn(1, hidden_size)
    
    # Forward pass with both initializations
    def forward(weights):
        activations = []
        h = x
        for w in weights:
            h = np.tanh(np.dot(h, w))
            activations.append(h)
        return activations
    
    std_activations = forward(std_weights)
    xavier_activations = forward(xavier_weights)
    
    print("Gradient magnitudes comparison:")
    print("Layer  Standard    Xavier")
    print("-" * 30)
    for i in range(n_layers):
        print(f"{i+1:2d}    {np.std(std_activations[i]):8.6f}  {np.std(xavier_activations[i]):8.6f}")

compare_gradients(100, 100, 5)

Gradient Clipping Implementation

def clip_gradients(gradients, max_norm):
    """
    Clips gradients to prevent exploding gradients while addressing vanishing gradients
    """
    # Calculate total norm of all gradients
    total_norm = np.sqrt(sum(np.sum(grad ** 2) for grad in gradients))
    
    # Calculate scaling factor
    clip_coef = max_norm / (total_norm + 1e-6)
    clip_coef = min(1.0, clip_coef)
    
    # Scale gradients
    clipped_gradients = [grad * clip_coef for grad in gradients]
    
    return clipped_gradients, total_norm

# Example usage
gradients = [np.random.randn(10, 10) * np.array([10**i for i in range(10)]) 
            for _ in range(3)]

max_norm = 5.0
clipped_grads, norm = clip_gradients(gradients, max_norm)

print(f"Original norm: {norm:.4f}")
print(f"After clipping:")
for i, grad in enumerate(clipped_grads):
    print(f"Layer {i+1} max gradient: {np.max(np.abs(grad)):.4f}")

LSTM Cell for Long-Term Dependencies

class LSTMCell:
    def __init__(self, input_size, hidden_size):
        # Initialize weights for input, forget, output gates and cell state
        self.Wf = xavier_init((input_size + hidden_size, hidden_size))
        self.Wi = xavier_init((input_size + hidden_size, hidden_size))
        self.Wc = xavier_init((input_size + hidden_size, hidden_size))
        self.Wo = xavier_init((input_size + hidden_size, hidden_size))
        
    def forward(self, x, prev_h, prev_c):
        # Concatenate input and previous hidden state
        combined = np.concatenate((x, prev_h), axis=1)
        
        # Gates computation
        forget_gate = sigmoid(np.dot(combined, self.Wf))
        input_gate = sigmoid(np.dot(combined, self.Wi))
        candidate = np.tanh(np.dot(combined, self.Wc))
        output_gate = sigmoid(np.dot(combined, self.Wo))
        
        # Update cell state
        next_c = forget_gate * prev_c + input_gate * candidate
        
        # Compute hidden state
        next_h = output_gate * np.tanh(next_c)
        
        return next_h, next_c

# Example usage
lstm = LSTMCell(10, 20)
x = np.random.randn(1, 10)
h = np.zeros((1, 20))
c = np.zeros((1, 20))

new_h, new_c = lstm.forward(x, h, c)
print(f"Hidden state shape: {new_h.shape}")
print(f"Cell state shape: {new_c.shape}")

Practical Example - Time Series Prediction with Gradient Monitoring

A complete implementation of a deep network for time series prediction, demonstrating gradient monitoring across layers to detect vanishing gradients in practice.

import numpy as np

class TimeSeriesPredictor:
    def __init__(self, input_size, hidden_sizes, output_size):
        self.layers = []
        self.gradients = []
        
        # Initialize layers with Xavier initialization
        layer_sizes = [input_size] + hidden_sizes + [output_size]
        for i in range(len(layer_sizes)-1):
            self.layers.append({
                'W': xavier_init((layer_sizes[i], layer_sizes[i+1])),
                'b': np.zeros(layer_sizes[i+1])
            })
    
    def forward(self, X):
        self.activations = [X]
        current_input = X
        
        # Forward pass with gradient monitoring
        for layer in self.layers[:-1]:
            z = np.dot(current_input, layer['W']) + layer['b']
            a = np.tanh(z)  # Using tanh to demonstrate vanishing gradients
            self.activations.append(a)
            current_input = a
            
        # Linear output layer
        output = np.dot(current_input, self.layers[-1]['W']) + self.layers[-1]['b']
        self.activations.append(output)
        return output
    
    def compute_gradients(self, y_true):
        batch_size = y_true.shape[0]
        grad_magnitudes = []
        
        # Compute gradients through backpropagation
        grad = (self.activations[-1] - y_true) / batch_size
        for i in reversed(range(len(self.layers))):
            grad_W = np.dot(self.activations[i].T, grad)
            grad_magnitudes.append(np.mean(np.abs(grad_W)))
            
            if i > 0:  # Not input layer
                grad = np.dot(grad, self.layers[i]['W'].T)
                grad = grad * (1 - self.activations[i] ** 2)  # tanh derivative
                
        return grad_magnitudes[::-1]

# Example usage with synthetic data
X = np.random.randn(100, 5)  # 100 samples, 5 features
y = np.sum(X, axis=1, keepdims=True)  # Simple sum target

model = TimeSeriesPredictor(5, [32, 16, 8], 1)
predictions = model.forward(X)
gradients = model.compute_gradients(y)

print("Gradient magnitudes per layer:")
for i, grad_mag in enumerate(gradients):
    print(f"Layer {i+1}: {grad_mag:.8f}")

Real-world Application - Stock Price Prediction with Gradient Monitoring

class StockPricePredictor:
    def __init__(self):
        self.model = TimeSeriesPredictor(10, [64, 32, 16], 1)
        self.gradient_history = []
        
    def preprocess_data(self, prices, window_size=10):
        X, y = [], []
        for i in range(len(prices) - window_size):
            X.append(prices[i:i+window_size])
            y.append(prices[i+window_size])
        return np.array(X), np.array(y).reshape(-1, 1)
    
    def train_epoch(self, X, y, learning_rate=0.001):
        predictions = self.model.forward(X)
        gradients = self.model.compute_gradients(y)
        self.gradient_history.append(gradients)
        return np.mean((predictions - y) ** 2)

# Generate synthetic stock data
np.random.seed(42)
days = 1000
stock_prices = np.cumsum(np.random.randn(days) * 0.02)  # Random walk
stock_prices += 100  # Start price at 100

# Prepare data
predictor = StockPricePredictor()
X, y = predictor.preprocess_data(stock_prices)

# Train for a few epochs
n_epochs = 5
for epoch in range(n_epochs):
    loss = predictor.train_epoch(X, y)
    
    # Analyze gradients
    current_gradients = predictor.gradient_history[-1]
    print(f"\nEpoch {epoch+1}:")
    print(f"Loss: {loss:.6f}")
    print("Layer gradients:")
    for i, grad in enumerate(current_gradients):
        print(f"Layer {i+1}: {grad:.8f}")

Results Analysis for Stock Price Prediction

import matplotlib.pyplot as plt

def analyze_gradient_evolution(gradient_history):
    n_epochs = len(gradient_history)
    n_layers = len(gradient_history[0])
    
    plt.figure(figsize=(10, 6))
    for layer in range(n_layers):
        layer_grads = [epoch_grads[layer] for epoch_grads in gradient_history]
        plt.plot(layer_grads, label=f'Layer {layer+1}')
    
    plt.xlabel('Epoch')
    plt.ylabel('Gradient Magnitude')
    plt.title('Gradient Evolution During Training')
    plt.legend()
    plt.yscale('log')  # Log scale to better visualize vanishing gradients
    
    # Calculate gradient ratios
    final_gradients = gradient_history[-1]
    print("\nFinal gradient ratios (relative to output layer):")
    output_grad = final_gradients[-1]
    for i, grad in enumerate(final_gradients):
        ratio = grad / output_grad
        print(f"Layer {i+1}: {ratio:.6f}")

# Analyze results
analyze_gradient_evolution(predictor.gradient_history)

# Calculate prediction accuracy metrics
test_predictions = predictor.model.forward(X[-100:])
test_actual = y[-100:]
mse = np.mean((test_predictions - test_actual) ** 2)
mae = np.mean(np.abs(test_predictions - test_actual))

print(f"\nTest Results:")
print(f"MSE: {mse:.6f}")
print(f"MAE: {mae:.6f}")

GRU Implementation for Gradient Flow Comparison

The Gated Recurrent Unit (GRU) offers an alternative to LSTM with fewer parameters while still maintaining good gradient flow. This implementation demonstrates its effectiveness against vanishing gradients.

class GRUCell:
    def __init__(self, input_size, hidden_size):
        # Reset gate parameters
        self.Wr = xavier_init((input_size + hidden_size, hidden_size))
        # Update gate parameters
        self.Wz = xavier_init((input_size + hidden_size, hidden_size))
        # Candidate hidden state parameters
        self.Wh = xavier_init((input_size + hidden_size, hidden_size))
        
        self.gradients = []
        
    def forward(self, x, prev_h):
        # Concatenate input and previous hidden state
        combined = np.concatenate((x, prev_h), axis=1)
        
        # Compute gates
        reset_gate = sigmoid(np.dot(combined, self.Wr))
        update_gate = sigmoid(np.dot(combined, self.Wz))
        
        # Compute candidate hidden state
        reset_hidden = reset_gate * prev_h
        combined_reset = np.concatenate((x, reset_hidden), axis=1)
        candidate = np.tanh(np.dot(combined_reset, self.Wh))
        
        # Compute new hidden state
        new_h = (1 - update_gate) * prev_h + update_gate * candidate
        
        # Store intermediate values for gradient tracking
        self.cache = (x, prev_h, reset_gate, update_gate, candidate)
        
        return new_h, [np.mean(np.abs(g)) for g in [reset_gate, update_gate, candidate]]

# Example usage and gradient analysis
gru = GRUCell(10, 20)
x = np.random.randn(1, 10)
h = np.zeros((1, 20))

new_h, gate_gradients = gru.forward(x, h)
print("Gate activation statistics:")
print(f"Reset gate mean gradient: {gate_gradients[0]:.6f}")
print(f"Update gate mean gradient: {gate_gradients[1]:.6f}")
print(f"Candidate mean gradient: {gate_gradients[2]:.6f}")

Comparative Analysis: LSTM vs GRU vs Vanilla RNN

class GradientComparisonExperiment:
    def __init__(self, sequence_length, input_size, hidden_size):
        self.lstm = LSTMCell(input_size, hidden_size)
        self.gru = GRUCell(input_size, hidden_size)
        self.vanilla_rnn = np.random.randn(input_size + hidden_size, hidden_size) * 0.01
        
        self.sequence_length = sequence_length
        self.hidden_size = hidden_size
        
    def run_comparison(self, input_sequence):
        results = {'lstm': [], 'gru': [], 'vanilla': []}
        
        # Initial states
        h_lstm = np.zeros((1, self.hidden_size))
        c_lstm = np.zeros((1, self.hidden_size))
        h_gru = np.zeros((1, self.hidden_size))
        h_vanilla = np.zeros((1, self.hidden_size))
        
        # Process sequence
        for t in range(len(input_sequence)):
            x_t = input_sequence[t:t+1]
            
            # LSTM forward
            h_lstm, c_lstm = self.lstm.forward(x_t, h_lstm, c_lstm)
            results['lstm'].append(np.mean(np.abs(h_lstm)))
            
            # GRU forward
            h_gru, gru_grads = self.gru.forward(x_t, h_gru)
            results['gru'].append(np.mean(np.abs(h_gru)))
            
            # Vanilla RNN forward
            combined = np.concatenate((x_t, h_vanilla), axis=1)
            h_vanilla = np.tanh(np.dot(combined, self.vanilla_rnn))
            results['vanilla'].append(np.mean(np.abs(h_vanilla)))
            
        return results

# Run experiment
seq_len, input_size, hidden_size = 100, 10, 20
input_sequence = np.random.randn(seq_len, 1, input_size)
experiment = GradientComparisonExperiment(seq_len, input_size, hidden_size)
results = experiment.run_comparison(input_sequence)

# Print analysis
for model in results:
    gradients = results[model]
    print(f"\n{model.upper()} Statistics:")
    print(f"Mean gradient magnitude: {np.mean(gradients):.6f}")
    print(f"Gradient magnitude std: {np.std(gradients):.6f}")
    print(f"Max/Min ratio: {max(gradients)/min(gradients):.6f}")

Additional Resources

More in Backpropagation and optimization

  1. Response:7 min

    In linear regression, the Sum of Squared Residuals (SSR) measures the total deviation between predicted and actual values. It serves as our cost function...

  2. Optimizing Neural Network Training with Learning Rate Control10 min

    Learning rate control is crucial for optimizing neural network performance. It involves adjusting the step size during gradient descent to find the right...

  3. RMSProp Optimization in Deep Learning with Python17 min

    RMSProp (Root Mean Square Propagation) is an optimization algorithm used in training deep neural networks. It addresses the issue of diminishing learning...

  4. Optimizers in Deep Learning Navigating the Path to Precision11 min

    Gradient descent forms the foundation of optimization in deep learning, acting as an iterative algorithm that minimizes the loss function by updating...

  5. Saddle Points in Deep Learning Optimization10 min

    Saddle points are critical points in the loss landscape of neural networks where the gradient is zero, but they are neither local minima nor maxima. They...

  6. Optimizers for Deep Learning Momentum, Nesterov, Adagrad, RMSProp, Adam9 min

    Optimizers play a crucial role in training deep neural networks. They are algorithms or methods used to adjust the attributes of the neural network, such as...