Deep RL Advances: Double DQN, Dueling DQN, PER, and Continuous Control¶

Goal: Understand and implement key improvements over standard DQN (Double, Dueling, Prioritized Experience Replay) and continuous control with DDPG from scratch in pure NumPy, with an outlook on SAC.

Prerequisites: Basics of RL, Q-learning, and standard DQN. Theory Link: Theory Markdown

1. Problem Setup — WHY¶

• Overestimation Bias: Standard DQN applies a max operator to noisy Q-value estimates during the Bellman update, leading to systemic overestimation of action values.
• Coupled State/Action Values: Flat Q-values entangle the general value of being in a state with the specific advantage of taking an action, which slows down learning across actions.
• Uniform Sampling: Standard replay buffers sample all transitions equally, wasting computation on transitions that offer little new information (low TD error).
• Continuous Control: DQN relies on $\arg\max_a Q(s,a)$, which is computationally intractable for continuous action spaces.

2. Mathematical Core — WHAT¶

Double DQN:

• Problem: $\mathbb{E}[\max_a Q(s,a)] \ge \max_a \mathbb{E}[Q(s,a)]$
• Fix: Decouple action selection (using online network $\theta$) and action evaluation (using target network $\theta^-$).
• Target: $y = r + \gamma Q_{\theta^-}(s', \arg\max_{a'} Q_\theta(s', a'))$

Dueling DQN:

• Decomposes Q-value: $Q(s,a;\theta) = V(s;\theta) + A(s,a;\theta) - \frac{1}{|A|} \sum_{a'} A(s,a';\theta)$
• Allows the network to learn which states are valuable independent of the action chosen.

Prioritized Experience Replay (PER):

• Priority: $p_i = |\delta_i| + \epsilon$
• Sampling Probability: $P(i) = \frac{p_i^\alpha}{\sum_k p_k^\alpha}$
• Importance Sampling Weight: $w_i = (\frac{1}{N \cdot P(i)})^\beta$

DDPG (implemented below):

• Deterministic policy (Actor): $a = \mu_\theta(s)$
• Critic: $Q_w(s,a)$
• Actor Gradient: $\nabla_\theta J \approx \mathbb{E}[\nabla_a Q_w(s,a)|_{a=\mu_\theta(s)} \cdot \nabla_\theta \mu_\theta(s)]$

Outlook: SAC (described only — not implemented in this notebook):

• Maximum Entropy Objective: $J(\pi) = \mathbb{E}[\sum \gamma^t (r_t + \alpha \mathcal{H}(\pi(\cdot|s_t)))]$
• Soft Bellman Backup: $Q(s,a) = r + \gamma \mathbb{E}[V(s')]$, where $V(s) = \mathbb{E}_{a}[Q(s,a) - \alpha \log \pi(a|s)]$

3. Solution Method — HOW¶

• Double DQN: Modify the target calculation to use two networks.
• Dueling DQN: Change the final layers of the NeuralQNetwork to split into $V$ and $A$ streams.
• PER: Implement a SumTree for $O(\log N)$ priority updates and sampling.
• DDPG: Implement simple Actor and Critic MLPs with manual backprop for continuous action control, using Ornstein-Uhlenbeck or Gaussian noise for exploration.

4. Implementation — BUILD¶

In [1]:

import random

import matplotlib.pyplot as plt
import numpy as np

%load_ext autoreload
%autoreload 2

SEED = 42
random.seed(SEED)
rng = np.random.default_rng(SEED)

plt.rcParams["figure.dpi"] = 90

import random import matplotlib.pyplot as plt import numpy as np %load_ext autoreload %autoreload 2 SEED = 42 random.seed(SEED) rng = np.random.default_rng(SEED) plt.rcParams["figure.dpi"] = 90

In [2]:

class GridWorld:
    def __init__(self):
        self.size = 4
        self.start = (0, 0)
        self.goal = (3, 3)
        self.obstacles = [(1, 1), (2, 1)]
        self.state = self.start
        
    def reset(self):
        self.state = self.start
        return self.state
        
    def step(self, action):
        # Actions: 0:up, 1:right, 2:down, 3:left
        x, y = self.state
        if action == 0: y = max(0, y - 1)
        elif action == 1: x = min(self.size - 1, x + 1)
        elif action == 2: y = min(self.size - 1, y + 1)
        elif action == 3: x = max(0, x - 1)
            
        if (x, y) not in self.obstacles:
            self.state = (x, y)
            
        done = (self.state == self.goal)
        reward = 1.0 if done else -0.01
        return self.state, reward, done

def state_to_onehot(state, size=4):
    idx = state[1] * size + state[0]
    vec = np.zeros(size * size)
    vec[idx] = 1.0
    return vec

class GridWorld: def __init__(self): self.size = 4 self.start = (0, 0) self.goal = (3, 3) self.obstacles = [(1, 1), (2, 1)] self.state = self.start def reset(self): self.state = self.start return self.state def step(self, action): # Actions: 0:up, 1:right, 2:down, 3:left x, y = self.state if action == 0: y = max(0, y - 1) elif action == 1: x = min(self.size - 1, x + 1) elif action == 2: y = min(self.size - 1, y + 1) elif action == 3: x = max(0, x - 1) if (x, y) not in self.obstacles: self.state = (x, y) done = (self.state == self.goal) reward = 1.0 if done else -0.01 return self.state, reward, done def state_to_onehot(state, size=4): idx = state[1] * size + state[0] vec = np.zeros(size * size) vec[idx] = 1.0 return vec

In [3]:

class NeuralQNetwork:
    def __init__(self, input_dim, hidden_dim, output_dim, lr=0.01):
        self.W1 = rng.standard_normal((input_dim, hidden_dim)) * np.sqrt(2 / input_dim)
        self.b1 = np.zeros(hidden_dim)
        self.W2 = rng.standard_normal((hidden_dim, output_dim)) * np.sqrt(2 / hidden_dim)
        self.b2 = np.zeros(output_dim)
        self.lr = lr

    def forward(self, x):
        self.x = np.atleast_2d(x)
        self.z1 = self.x @ self.W1 + self.b1
        self.a1 = np.maximum(0, self.z1)
        self.z2 = self.a1 @ self.W2 + self.b2
        return self.z2

    def backward(self, grad_z2):
        grad_W2 = self.a1.T @ grad_z2
        grad_b2 = np.sum(grad_z2, axis=0)
        grad_a1 = grad_z2 @ self.W2.T
        grad_z1 = grad_a1 * (self.z1 > 0)
        grad_W1 = self.x.T @ grad_z1
        grad_b1 = np.sum(grad_z1, axis=0)

        self.W2 -= self.lr * grad_W2
        self.b2 -= self.lr * grad_b2
        self.W1 -= self.lr * grad_W1
        self.b1 -= self.lr * grad_b1

    def copy_weights_from(self, other):
        self.W1, self.b1 = other.W1.copy(), other.b1.copy()
        self.W2, self.b2 = other.W2.copy(), other.b2.copy()

class DuelingNetwork(NeuralQNetwork):
    def __init__(self, input_dim, hidden_dim, output_dim, lr=0.01):
        super().__init__(input_dim, hidden_dim, output_dim, lr)
        self.W_V = rng.standard_normal((hidden_dim, 1)) * np.sqrt(2 / hidden_dim)
        self.b_V = np.zeros(1)
        self.W_A = rng.standard_normal((hidden_dim, output_dim)) * np.sqrt(2 / hidden_dim)
        self.b_A = np.zeros(output_dim)
        self.output_dim = output_dim

    def forward(self, x):
        self.x = np.atleast_2d(x)
        self.z1 = self.x @ self.W1 + self.b1
        self.a1 = np.maximum(0, self.z1)
        self.V = self.a1 @ self.W_V + self.b_V
        self.A = self.a1 @ self.W_A + self.b_A
        self.mean_A = np.mean(self.A, axis=1, keepdims=True)
        self.Q = self.V + self.A - self.mean_A
        return self.Q

    def backward(self, grad_Q):
        # grad_Q shape: (batch, output_dim)
        grad_V = np.sum(grad_Q, axis=1, keepdims=True)  # (batch, 1)
        grad_A = grad_Q - np.mean(grad_Q, axis=1, keepdims=True)  # (batch, output_dim)

        grad_W_V = self.a1.T @ grad_V
        grad_b_V = np.sum(grad_V, axis=0)
        grad_W_A = self.a1.T @ grad_A
        grad_b_A = np.sum(grad_A, axis=0)

        grad_a1 = grad_V @ self.W_V.T + grad_A @ self.W_A.T
        grad_z1 = grad_a1 * (self.z1 > 0)
        grad_W1 = self.x.T @ grad_z1
        grad_b1 = np.sum(grad_z1, axis=0)

        self.W_V -= self.lr * grad_W_V
        self.b_V -= self.lr * grad_b_V
        self.W_A -= self.lr * grad_W_A
        self.b_A -= self.lr * grad_b_A
        self.W1 -= self.lr * grad_W1
        self.b1 -= self.lr * grad_b1
        
    def copy_weights_from(self, other):
        self.W1, self.b1 = other.W1.copy(), other.b1.copy()
        self.W_V, self.b_V = other.W_V.copy(), other.b_V.copy()
        self.W_A, self.b_A = other.W_A.copy(), other.b_A.copy()

class NeuralQNetwork: def __init__(self, input_dim, hidden_dim, output_dim, lr=0.01): self.W1 = rng.standard_normal((input_dim, hidden_dim)) * np.sqrt(2 / input_dim) self.b1 = np.zeros(hidden_dim) self.W2 = rng.standard_normal((hidden_dim, output_dim)) * np.sqrt(2 / hidden_dim) self.b2 = np.zeros(output_dim) self.lr = lr def forward(self, x): self.x = np.atleast_2d(x) self.z1 = self.x @ self.W1 + self.b1 self.a1 = np.maximum(0, self.z1) self.z2 = self.a1 @ self.W2 + self.b2 return self.z2 def backward(self, grad_z2): grad_W2 = self.a1.T @ grad_z2 grad_b2 = np.sum(grad_z2, axis=0) grad_a1 = grad_z2 @ self.W2.T grad_z1 = grad_a1 * (self.z1 > 0) grad_W1 = self.x.T @ grad_z1 grad_b1 = np.sum(grad_z1, axis=0) self.W2 -= self.lr * grad_W2 self.b2 -= self.lr * grad_b2 self.W1 -= self.lr * grad_W1 self.b1 -= self.lr * grad_b1 def copy_weights_from(self, other): self.W1, self.b1 = other.W1.copy(), other.b1.copy() self.W2, self.b2 = other.W2.copy(), other.b2.copy() class DuelingNetwork(NeuralQNetwork): def __init__(self, input_dim, hidden_dim, output_dim, lr=0.01): super().__init__(input_dim, hidden_dim, output_dim, lr) self.W_V = rng.standard_normal((hidden_dim, 1)) * np.sqrt(2 / hidden_dim) self.b_V = np.zeros(1) self.W_A = rng.standard_normal((hidden_dim, output_dim)) * np.sqrt(2 / hidden_dim) self.b_A = np.zeros(output_dim) self.output_dim = output_dim def forward(self, x): self.x = np.atleast_2d(x) self.z1 = self.x @ self.W1 + self.b1 self.a1 = np.maximum(0, self.z1) self.V = self.a1 @ self.W_V + self.b_V self.A = self.a1 @ self.W_A + self.b_A self.mean_A = np.mean(self.A, axis=1, keepdims=True) self.Q = self.V + self.A - self.mean_A return self.Q def backward(self, grad_Q): # grad_Q shape: (batch, output_dim) grad_V = np.sum(grad_Q, axis=1, keepdims=True) # (batch, 1) grad_A = grad_Q - np.mean(grad_Q, axis=1, keepdims=True) # (batch, output_dim) grad_W_V = self.a1.T @ grad_V grad_b_V = np.sum(grad_V, axis=0) grad_W_A = self.a1.T @ grad_A grad_b_A = np.sum(grad_A, axis=0) grad_a1 = grad_V @ self.W_V.T + grad_A @ self.W_A.T grad_z1 = grad_a1 * (self.z1 > 0) grad_W1 = self.x.T @ grad_z1 grad_b1 = np.sum(grad_z1, axis=0) self.W_V -= self.lr * grad_W_V self.b_V -= self.lr * grad_b_V self.W_A -= self.lr * grad_W_A self.b_A -= self.lr * grad_b_A self.W1 -= self.lr * grad_W1 self.b1 -= self.lr * grad_b1 def copy_weights_from(self, other): self.W1, self.b1 = other.W1.copy(), other.b1.copy() self.W_V, self.b_V = other.W_V.copy(), other.b_V.copy() self.W_A, self.b_A = other.W_A.copy(), other.b_A.copy()

In [4]:

class ReplayBuffer:
    def __init__(self, capacity):
        self.capacity = capacity
        self.buffer = []
        self.pos = 0

    def push(self, transition):
        if len(self.buffer) < self.capacity:
            self.buffer.append(transition)
        else:
            self.buffer[self.pos] = transition
        self.pos = (self.pos + 1) % self.capacity

    def sample(self, batch_size):
        indices = rng.choice(len(self.buffer), batch_size, replace=False)
        batch = [self.buffer[i] for i in indices]
        return map(np.array, zip(*batch))

class SumTree:
    def __init__(self, capacity):
        self.capacity = capacity
        self.tree = np.zeros(2 * capacity - 1)
        self.data = np.zeros(capacity, dtype=object)
        self.write = 0
        self.n_entries = 0
        
    def total(self):
        return self.tree[0]
        
    def update(self, idx, p):
        change = p - self.tree[idx]
        self.tree[idx] = p
        while idx != 0:
            idx = (idx - 1) // 2
            self.tree[idx] += change
            
    def add(self, p, data):
        idx = self.write + self.capacity - 1
        self.data[self.write] = data
        self.update(idx, p)
        self.write = (self.write + 1) % self.capacity
        if self.n_entries < self.capacity:
            self.n_entries += 1
            
    def get(self, s):
        idx = 0
        while True:
            left = 2 * idx + 1
            right = left + 1
            if left >= len(self.tree):
                break
            if s <= self.tree[left]:
                idx = left
            else:
                s -= self.tree[left]
                idx = right
        data_idx = idx - self.capacity + 1
        return idx, self.tree[idx], self.data[data_idx]

class PrioritizedReplayBuffer:
    def __init__(self, capacity, alpha=0.6):
        self.tree = SumTree(capacity)
        self.alpha = alpha
        self.capacity = capacity
        self.max_priority = 1.0
        
    def push(self, transition):
        priority = self.max_priority ** self.alpha
        self.tree.add(priority, transition)
        
    def sample(self, batch_size, beta=0.4):
        batch = []
        indices = []
        weights = []
        segment = self.tree.total() / batch_size
        
        for i in range(batch_size):
            a = segment * i
            b = segment * (i + 1)
            s = rng.uniform(a, b)
            idx, p, data = self.tree.get(s)
            batch.append(data)
            indices.append(idx)
            prob = p / self.tree.total()
            weights.append((self.tree.n_entries * prob) ** (-beta))
            
        weights = np.array(weights) / max(weights)
        return map(np.array, zip(*batch)), indices, weights
        
    def update_priorities(self, indices, td_errors):
        for idx, err in zip(indices, td_errors):
            p = (abs(err) + 1e-5) ** self.alpha
            self.tree.update(idx, p)
            self.max_priority = max(self.max_priority, p)

class ReplayBuffer: def __init__(self, capacity): self.capacity = capacity self.buffer = [] self.pos = 0 def push(self, transition): if len(self.buffer) < self.capacity: self.buffer.append(transition) else: self.buffer[self.pos] = transition self.pos = (self.pos + 1) % self.capacity def sample(self, batch_size): indices = rng.choice(len(self.buffer), batch_size, replace=False) batch = [self.buffer[i] for i in indices] return map(np.array, zip(*batch)) class SumTree: def __init__(self, capacity): self.capacity = capacity self.tree = np.zeros(2 * capacity - 1) self.data = np.zeros(capacity, dtype=object) self.write = 0 self.n_entries = 0 def total(self): return self.tree[0] def update(self, idx, p): change = p - self.tree[idx] self.tree[idx] = p while idx != 0: idx = (idx - 1) // 2 self.tree[idx] += change def add(self, p, data): idx = self.write + self.capacity - 1 self.data[self.write] = data self.update(idx, p) self.write = (self.write + 1) % self.capacity if self.n_entries < self.capacity: self.n_entries += 1 def get(self, s): idx = 0 while True: left = 2 * idx + 1 right = left + 1 if left >= len(self.tree): break if s <= self.tree[left]: idx = left else: s -= self.tree[left] idx = right data_idx = idx - self.capacity + 1 return idx, self.tree[idx], self.data[data_idx] class PrioritizedReplayBuffer: def __init__(self, capacity, alpha=0.6): self.tree = SumTree(capacity) self.alpha = alpha self.capacity = capacity self.max_priority = 1.0 def push(self, transition): priority = self.max_priority ** self.alpha self.tree.add(priority, transition) def sample(self, batch_size, beta=0.4): batch = [] indices = [] weights = [] segment = self.tree.total() / batch_size for i in range(batch_size): a = segment * i b = segment * (i + 1) s = rng.uniform(a, b) idx, p, data = self.tree.get(s) batch.append(data) indices.append(idx) prob = p / self.tree.total() weights.append((self.tree.n_entries * prob) ** (-beta)) weights = np.array(weights) / max(weights) return map(np.array, zip(*batch)), indices, weights def update_priorities(self, indices, td_errors): for idx, err in zip(indices, td_errors): p = (abs(err) + 1e-5) ** self.alpha self.tree.update(idx, p) self.max_priority = max(self.max_priority, p)

In [5]:

def train_dqn(env, agent_type='dqn', episodes=200, use_per=False):
    input_dim = env.size * env.size
    output_dim = 4
    if agent_type == 'dueling':
        q_net = DuelingNetwork(input_dim, 32, output_dim, lr=0.01)
        target_net = DuelingNetwork(input_dim, 32, output_dim, lr=0.01)
    else:
        q_net = NeuralQNetwork(input_dim, 32, output_dim, lr=0.01)
        target_net = NeuralQNetwork(input_dim, 32, output_dim, lr=0.01)
    
    target_net.copy_weights_from(q_net)
    
    if use_per:
        buffer = PrioritizedReplayBuffer(1000)
    else:
        buffer = ReplayBuffer(1000)
        
    epsilon = 1.0
    gamma = 0.95
    batch_size = 32
    
    returns = []
    q_values_record = []
    
    for ep in range(episodes):
        state = env.reset()
        s_vec = state_to_onehot(state)
        total_reward = 0
        done = False
        step = 0
        
        # Track start state Q value for overestimation analysis
        start_q = np.max(q_net.forward(s_vec)[0])
        q_values_record.append(start_q)
        
        while not done and step < 50:
            if rng.random() < epsilon:
                action = rng.integers(0, 4)
            else:
                action = np.argmax(q_net.forward(s_vec)[0])
                
            next_state, reward, done = env.step(action)
            ns_vec = state_to_onehot(next_state)
            
            buffer.push((s_vec, action, reward, ns_vec, done))
            
            if (use_per and buffer.tree.n_entries >= batch_size) or (not use_per and len(buffer.buffer) >= batch_size):
                if use_per:
                    (b_s, b_a, b_r, b_ns, b_d), indices, weights = buffer.sample(batch_size)
                else:
                    b_s, b_a, b_r, b_ns, b_d = buffer.sample(batch_size)
                    weights = np.ones(batch_size)
                    
                q_curr = q_net.forward(b_s)
                
                if agent_type == 'double_dqn':
                    # Double DQN: action selected by online net, evaluated by target net
                    next_actions = np.argmax(q_net.forward(b_ns), axis=1)
                    q_next = target_net.forward(b_ns)
                    targets = b_r + gamma * q_next[np.arange(batch_size), next_actions] * (1 - b_d)
                else:
                    # Standard DQN / Dueling (unless combined with Double)
                    q_next = target_net.forward(b_ns)
                    targets = b_r + gamma * np.max(q_next, axis=1) * (1 - b_d)
                    
                # Compute TD errors for PER
                td_errors = targets - q_curr[np.arange(batch_size), b_a]
                if use_per:
                    buffer.update_priorities(indices, td_errors)
                
                grad = np.zeros_like(q_curr)
                # Include importance sampling weights
                grad[np.arange(batch_size), b_a] = -td_errors * weights
                q_net.backward(grad / batch_size)
                
            s_vec = ns_vec
            total_reward += reward
            step += 1
            
        epsilon = max(0.01, epsilon * 0.95)
        returns.append(total_reward)
        if ep % 20 == 0:
            target_net.copy_weights_from(q_net)
            
    return returns, q_values_record

def train_dqn(env, agent_type='dqn', episodes=200, use_per=False): input_dim = env.size * env.size output_dim = 4 if agent_type == 'dueling': q_net = DuelingNetwork(input_dim, 32, output_dim, lr=0.01) target_net = DuelingNetwork(input_dim, 32, output_dim, lr=0.01) else: q_net = NeuralQNetwork(input_dim, 32, output_dim, lr=0.01) target_net = NeuralQNetwork(input_dim, 32, output_dim, lr=0.01) target_net.copy_weights_from(q_net) if use_per: buffer = PrioritizedReplayBuffer(1000) else: buffer = ReplayBuffer(1000) epsilon = 1.0 gamma = 0.95 batch_size = 32 returns = [] q_values_record = [] for ep in range(episodes): state = env.reset() s_vec = state_to_onehot(state) total_reward = 0 done = False step = 0 # Track start state Q value for overestimation analysis start_q = np.max(q_net.forward(s_vec)[0]) q_values_record.append(start_q) while not done and step < 50: if rng.random() < epsilon: action = rng.integers(0, 4) else: action = np.argmax(q_net.forward(s_vec)[0]) next_state, reward, done = env.step(action) ns_vec = state_to_onehot(next_state) buffer.push((s_vec, action, reward, ns_vec, done)) if (use_per and buffer.tree.n_entries >= batch_size) or (not use_per and len(buffer.buffer) >= batch_size): if use_per: (b_s, b_a, b_r, b_ns, b_d), indices, weights = buffer.sample(batch_size) else: b_s, b_a, b_r, b_ns, b_d = buffer.sample(batch_size) weights = np.ones(batch_size) q_curr = q_net.forward(b_s) if agent_type == 'double_dqn': # Double DQN: action selected by online net, evaluated by target net next_actions = np.argmax(q_net.forward(b_ns), axis=1) q_next = target_net.forward(b_ns) targets = b_r + gamma * q_next[np.arange(batch_size), next_actions] * (1 - b_d) else: # Standard DQN / Dueling (unless combined with Double) q_next = target_net.forward(b_ns) targets = b_r + gamma * np.max(q_next, axis=1) * (1 - b_d) # Compute TD errors for PER td_errors = targets - q_curr[np.arange(batch_size), b_a] if use_per: buffer.update_priorities(indices, td_errors) grad = np.zeros_like(q_curr) # Include importance sampling weights grad[np.arange(batch_size), b_a] = -td_errors * weights q_net.backward(grad / batch_size) s_vec = ns_vec total_reward += reward step += 1 epsilon = max(0.01, epsilon * 0.95) returns.append(total_reward) if ep % 20 == 0: target_net.copy_weights_from(q_net) return returns, q_values_record

In [6]:

class Pendulum:
    def __init__(self):
        self.max_speed = 8.0
        self.max_torque = 2.0
        self.dt = 0.05
        self.g = 10.0
        self.m = 1.0
        self.l = 1.0
        self.state = np.array([np.pi, 0.0]) # [theta, theta_dot]
        
    def reset(self):
        self.state = rng.uniform(low=[-np.pi, -1], high=[np.pi, 1])
        return self._get_obs()
        
    def _get_obs(self):
        theta, thetadot = self.state
        return np.array([np.cos(theta), np.sin(theta), thetadot])
        
    def step(self, action):
        torque = np.clip(action, -self.max_torque, self.max_torque)[0]
        th, thdot = self.state
        
        # Reward: penalize being away from upright, high speed, and high effort
        cost = ((th + np.pi) % (2 * np.pi) - np.pi)**2 + 0.1 * thdot**2 + 0.001 * torque**2
        reward = -cost
        
        new_thdot = thdot + (-3 * self.g / (2 * self.l) * np.sin(th + np.pi) + 3.0 / (self.m * self.l**2) * torque) * self.dt
        new_thdot = np.clip(new_thdot, -self.max_speed, self.max_speed)
        new_th = th + new_thdot * self.dt
        
        self.state = np.array([new_th, new_thdot])
        return self._get_obs(), reward, False

class Pendulum: def __init__(self): self.max_speed = 8.0 self.max_torque = 2.0 self.dt = 0.05 self.g = 10.0 self.m = 1.0 self.l = 1.0 self.state = np.array([np.pi, 0.0]) # [theta, theta_dot] def reset(self): self.state = rng.uniform(low=[-np.pi, -1], high=[np.pi, 1]) return self._get_obs() def _get_obs(self): theta, thetadot = self.state return np.array([np.cos(theta), np.sin(theta), thetadot]) def step(self, action): torque = np.clip(action, -self.max_torque, self.max_torque)[0] th, thdot = self.state # Reward: penalize being away from upright, high speed, and high effort cost = ((th + np.pi) % (2 * np.pi) - np.pi)**2 + 0.1 * thdot**2 + 0.001 * torque**2 reward = -cost new_thdot = thdot + (-3 * self.g / (2 * self.l) * np.sin(th + np.pi) + 3.0 / (self.m * self.l**2) * torque) * self.dt new_thdot = np.clip(new_thdot, -self.max_speed, self.max_speed) new_th = th + new_thdot * self.dt self.state = np.array([new_th, new_thdot]) return self._get_obs(), reward, False

In [7]:

class DDPGActor:
    def __init__(self, state_dim, hidden_dim, action_dim, action_bound, lr=1e-3):
        self.W1 = rng.standard_normal((state_dim, hidden_dim)) * np.sqrt(2 / state_dim)
        self.b1 = np.zeros(hidden_dim)
        self.W2 = rng.standard_normal((hidden_dim, action_dim)) * np.sqrt(2 / hidden_dim)
        self.b2 = np.zeros(action_dim)
        self.lr = lr
        self.action_bound = action_bound

    def forward(self, x):
        self.x = np.atleast_2d(x)
        self.z1 = self.x @ self.W1 + self.b1
        self.a1 = np.maximum(0, self.z1)
        self.z2 = self.a1 @ self.W2 + self.b2
        self.out = np.tanh(self.z2) * self.action_bound
        return self.out

    def backward(self, grad_action):
        # grad_action is dJ/da
        grad_z2 = grad_action * self.action_bound * (1 - np.tanh(self.z2)**2)
        grad_W2 = self.a1.T @ grad_z2
        grad_b2 = np.sum(grad_z2, axis=0)
        grad_a1 = grad_z2 @ self.W2.T
        grad_z1 = grad_a1 * (self.z1 > 0)
        grad_W1 = self.x.T @ grad_z1
        grad_b1 = np.sum(grad_z1, axis=0)

        self.W2 += self.lr * grad_W2  # ascent to maximize Q
        self.b2 += self.lr * grad_b2
        self.W1 += self.lr * grad_W1
        self.b1 += self.lr * grad_b1
        
    def copy_weights_from(self, other, tau=1.0):
        self.W1 = tau * other.W1 + (1 - tau) * self.W1
        self.b1 = tau * other.b1 + (1 - tau) * self.b1
        self.W2 = tau * other.W2 + (1 - tau) * self.W2
        self.b2 = tau * other.b2 + (1 - tau) * self.b2

class DDPGCritic:
    def __init__(self, state_dim, action_dim, hidden_dim, lr=1e-3):
        self.W1 = rng.standard_normal((state_dim + action_dim, hidden_dim)) * np.sqrt(2 / (state_dim + action_dim))
        self.b1 = np.zeros(hidden_dim)
        self.W2 = rng.standard_normal((hidden_dim, 1)) * np.sqrt(2 / hidden_dim)
        self.b2 = np.zeros(1)
        self.lr = lr
        self.state_dim = state_dim

    def forward(self, state, action):
        self.x = np.concatenate([np.atleast_2d(state), np.atleast_2d(action)], axis=1)
        self.z1 = self.x @ self.W1 + self.b1
        self.a1 = np.maximum(0, self.z1)
        self.z2 = self.a1 @ self.W2 + self.b2
        return self.z2

    def backward(self, td_error):
        # td_error = Q - target
        grad_z2 = td_error
        grad_W2 = self.a1.T @ grad_z2
        grad_b2 = np.sum(grad_z2, axis=0)
        grad_a1 = grad_z2 @ self.W2.T
        grad_z1 = grad_a1 * (self.z1 > 0)
        grad_W1 = self.x.T @ grad_z1
        grad_b1 = np.sum(grad_z1, axis=0)

        self.W2 -= self.lr * grad_W2
        self.b2 -= self.lr * grad_b2
        self.W1 -= self.lr * grad_W1
        self.b1 -= self.lr * grad_b1
        
        # Return gradient wrt action for actor update
        grad_x = grad_z1 @ self.W1.T
        return grad_x[:, self.state_dim:]
        
    def copy_weights_from(self, other, tau=1.0):
        self.W1 = tau * other.W1 + (1 - tau) * self.W1
        self.b1 = tau * other.b1 + (1 - tau) * self.b1
        self.W2 = tau * other.W2 + (1 - tau) * self.W2
        self.b2 = tau * other.b2 + (1 - tau) * self.b2

class DDPGActor: def __init__(self, state_dim, hidden_dim, action_dim, action_bound, lr=1e-3): self.W1 = rng.standard_normal((state_dim, hidden_dim)) * np.sqrt(2 / state_dim) self.b1 = np.zeros(hidden_dim) self.W2 = rng.standard_normal((hidden_dim, action_dim)) * np.sqrt(2 / hidden_dim) self.b2 = np.zeros(action_dim) self.lr = lr self.action_bound = action_bound def forward(self, x): self.x = np.atleast_2d(x) self.z1 = self.x @ self.W1 + self.b1 self.a1 = np.maximum(0, self.z1) self.z2 = self.a1 @ self.W2 + self.b2 self.out = np.tanh(self.z2) * self.action_bound return self.out def backward(self, grad_action): # grad_action is dJ/da grad_z2 = grad_action * self.action_bound * (1 - np.tanh(self.z2)**2) grad_W2 = self.a1.T @ grad_z2 grad_b2 = np.sum(grad_z2, axis=0) grad_a1 = grad_z2 @ self.W2.T grad_z1 = grad_a1 * (self.z1 > 0) grad_W1 = self.x.T @ grad_z1 grad_b1 = np.sum(grad_z1, axis=0) self.W2 += self.lr * grad_W2 # ascent to maximize Q self.b2 += self.lr * grad_b2 self.W1 += self.lr * grad_W1 self.b1 += self.lr * grad_b1 def copy_weights_from(self, other, tau=1.0): self.W1 = tau * other.W1 + (1 - tau) * self.W1 self.b1 = tau * other.b1 + (1 - tau) * self.b1 self.W2 = tau * other.W2 + (1 - tau) * self.W2 self.b2 = tau * other.b2 + (1 - tau) * self.b2 class DDPGCritic: def __init__(self, state_dim, action_dim, hidden_dim, lr=1e-3): self.W1 = rng.standard_normal((state_dim + action_dim, hidden_dim)) * np.sqrt(2 / (state_dim + action_dim)) self.b1 = np.zeros(hidden_dim) self.W2 = rng.standard_normal((hidden_dim, 1)) * np.sqrt(2 / hidden_dim) self.b2 = np.zeros(1) self.lr = lr self.state_dim = state_dim def forward(self, state, action): self.x = np.concatenate([np.atleast_2d(state), np.atleast_2d(action)], axis=1) self.z1 = self.x @ self.W1 + self.b1 self.a1 = np.maximum(0, self.z1) self.z2 = self.a1 @ self.W2 + self.b2 return self.z2 def backward(self, td_error): # td_error = Q - target grad_z2 = td_error grad_W2 = self.a1.T @ grad_z2 grad_b2 = np.sum(grad_z2, axis=0) grad_a1 = grad_z2 @ self.W2.T grad_z1 = grad_a1 * (self.z1 > 0) grad_W1 = self.x.T @ grad_z1 grad_b1 = np.sum(grad_z1, axis=0) self.W2 -= self.lr * grad_W2 self.b2 -= self.lr * grad_b2 self.W1 -= self.lr * grad_W1 self.b1 -= self.lr * grad_b1 # Return gradient wrt action for actor update grad_x = grad_z1 @ self.W1.T return grad_x[:, self.state_dim:] def copy_weights_from(self, other, tau=1.0): self.W1 = tau * other.W1 + (1 - tau) * self.W1 self.b1 = tau * other.b1 + (1 - tau) * self.b1 self.W2 = tau * other.W2 + (1 - tau) * self.W2 self.b2 = tau * other.b2 + (1 - tau) * self.b2

In [8]:

def train_ddpg(env, episodes=50, tau=0.01):
    state_dim = 3
    action_dim = 1
    action_bound = 2.0
    
    actor = DDPGActor(state_dim, 32, action_dim, action_bound, lr=1e-3)
    actor_target = DDPGActor(state_dim, 32, action_dim, action_bound)
    actor_target.copy_weights_from(actor)
    
    critic = DDPGCritic(state_dim, action_dim, 32, lr=5e-3)
    critic_target = DDPGCritic(state_dim, action_dim, 32)
    critic_target.copy_weights_from(critic)
    
    buffer = ReplayBuffer(5000)
    batch_size = 64
    gamma = 0.99
    
    returns = []
    
    for ep in range(episodes):
        state = env.reset()
        total_reward = 0
        
        for step in range(100):
            # Action selection with exploration noise
            action = actor.forward(state)[0] + rng.normal(0, 0.2, size=action_dim)
            action = np.clip(action, -action_bound, action_bound)
            
            next_state, reward, done = env.step(action)
            buffer.push((state, action, reward, next_state, done))
            
            if len(buffer.buffer) > batch_size:
                b_s, b_a, b_r, b_ns, b_d = buffer.sample(batch_size)
                
                # Critic update
                next_actions = actor_target.forward(b_ns)
                q_next = critic_target.forward(b_ns, next_actions)
                b_r_col = b_r.reshape(-1, 1)
                b_d_col = b_d.reshape(-1, 1)
                targets = b_r_col + gamma * q_next * (1 - b_d_col)
                
                q_curr = critic.forward(b_s, b_a)
                td_error = q_curr - targets
                critic.backward(td_error / batch_size)
                
                # Actor update
                actions_pred = actor.forward(b_s)
                critic.forward(b_s, actions_pred)
                # We want to maximize Q, so we take gradient of Q wrt action
                grad_action = critic.backward(np.zeros_like(q_curr)) # just to get grad_x
                # In simple manual backprop, the actor needs dQ/da
                actor.backward(grad_action / batch_size)
                
                # Soft target updates
                actor_target.copy_weights_from(actor, tau)
                critic_target.copy_weights_from(critic, tau)
                
            state = next_state
            total_reward += reward
            
        returns.append(total_reward)
        
    return returns

def train_ddpg(env, episodes=50, tau=0.01): state_dim = 3 action_dim = 1 action_bound = 2.0 actor = DDPGActor(state_dim, 32, action_dim, action_bound, lr=1e-3) actor_target = DDPGActor(state_dim, 32, action_dim, action_bound) actor_target.copy_weights_from(actor) critic = DDPGCritic(state_dim, action_dim, 32, lr=5e-3) critic_target = DDPGCritic(state_dim, action_dim, 32) critic_target.copy_weights_from(critic) buffer = ReplayBuffer(5000) batch_size = 64 gamma = 0.99 returns = [] for ep in range(episodes): state = env.reset() total_reward = 0 for step in range(100): # Action selection with exploration noise action = actor.forward(state)[0] + rng.normal(0, 0.2, size=action_dim) action = np.clip(action, -action_bound, action_bound) next_state, reward, done = env.step(action) buffer.push((state, action, reward, next_state, done)) if len(buffer.buffer) > batch_size: b_s, b_a, b_r, b_ns, b_d = buffer.sample(batch_size) # Critic update next_actions = actor_target.forward(b_ns) q_next = critic_target.forward(b_ns, next_actions) b_r_col = b_r.reshape(-1, 1) b_d_col = b_d.reshape(-1, 1) targets = b_r_col + gamma * q_next * (1 - b_d_col) q_curr = critic.forward(b_s, b_a) td_error = q_curr - targets critic.backward(td_error / batch_size) # Actor update actions_pred = actor.forward(b_s) critic.forward(b_s, actions_pred) # We want to maximize Q, so we take gradient of Q wrt action grad_action = critic.backward(np.zeros_like(q_curr)) # just to get grad_x # In simple manual backprop, the actor needs dQ/da actor.backward(grad_action / batch_size) # Soft target updates actor_target.copy_weights_from(actor, tau) critic_target.copy_weights_from(critic, tau) state = next_state total_reward += reward returns.append(total_reward) return returns

5. Library Comparison / Experiments¶

No standard library ships these DQN variants as pure-NumPy references, so we compare the variants against each other and check Q-value estimates against the known optimal start-state value. Let's run comparisons between DQN, Double DQN, and Dueling DQN.

In [9]:

env = GridWorld()
print("Training standard DQN...")
dqn_returns, dqn_q = train_dqn(env, 'dqn', episodes=200)
print("Training Double DQN...")
ddqn_returns, ddqn_q = train_dqn(env, 'double_dqn', episodes=200)
print("Training Dueling DQN...")
duel_returns, duel_q = train_dqn(env, 'dueling', episodes=200)
print("Training PER (Standard DQN with PER)...")
per_returns, per_q = train_dqn(env, 'dqn', episodes=200, use_per=True)

def moving_avg(a, n=10):
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n - 1:] / n

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(moving_avg(dqn_returns), label='DQN')
plt.plot(moving_avg(ddqn_returns), label='Double DQN')
plt.plot(moving_avg(duel_returns), label='Dueling DQN')
plt.plot(moving_avg(per_returns), label='PER')
plt.xlabel('Episodes')
plt.ylabel('Smoothed Returns')
plt.title('Learning Curves on GridWorld')
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(dqn_q, label='DQN Q-value (Start state)', alpha=0.7)
plt.plot(ddqn_q, label='Double DQN Q-value', alpha=0.7)
plt.axhline(0.95**6, color='r', linestyle='--', label='Approx True Value')
plt.xlabel('Episodes')
plt.ylabel('Max Q-value at Start')
plt.title('Overestimation Bias Check')
plt.legend()
plt.tight_layout()
plt.show()

env = GridWorld() print("Training standard DQN...") dqn_returns, dqn_q = train_dqn(env, 'dqn', episodes=200) print("Training Double DQN...") ddqn_returns, ddqn_q = train_dqn(env, 'double_dqn', episodes=200) print("Training Dueling DQN...") duel_returns, duel_q = train_dqn(env, 'dueling', episodes=200) print("Training PER (Standard DQN with PER)...") per_returns, per_q = train_dqn(env, 'dqn', episodes=200, use_per=True) def moving_avg(a, n=10): ret = np.cumsum(a, dtype=float) ret[n:] = ret[n:] - ret[:-n] return ret[n - 1:] / n plt.figure(figsize=(12, 5)) plt.subplot(1, 2, 1) plt.plot(moving_avg(dqn_returns), label='DQN') plt.plot(moving_avg(ddqn_returns), label='Double DQN') plt.plot(moving_avg(duel_returns), label='Dueling DQN') plt.plot(moving_avg(per_returns), label='PER') plt.xlabel('Episodes') plt.ylabel('Smoothed Returns') plt.title('Learning Curves on GridWorld') plt.legend() plt.subplot(1, 2, 2) plt.plot(dqn_q, label='DQN Q-value (Start state)', alpha=0.7) plt.plot(ddqn_q, label='Double DQN Q-value', alpha=0.7) plt.axhline(0.95**6, color='r', linestyle='--', label='Approx True Value') plt.xlabel('Episodes') plt.ylabel('Max Q-value at Start') plt.title('Overestimation Bias Check') plt.legend() plt.tight_layout() plt.show()

Training standard DQN...

Training Double DQN...

Training Dueling DQN...

Training PER (Standard DQN with PER)...

6. Experiments and Failures — VERIFY¶

Continuous Control with DDPG¶

Let's test DDPG on the Pendulum environment.

In [10]:

pendulum = Pendulum()
print("Training DDPG on Pendulum (with target soft update)...")
ddpg_returns = train_ddpg(pendulum, episodes=100, tau=0.01)

print("Training DDPG on Pendulum (FAIL: without soft update, tau=1.0)...")
ddpg_fail_returns = train_ddpg(pendulum, episodes=100, tau=1.0)

plt.figure(figsize=(6, 4))
plt.plot(moving_avg(ddpg_returns), label='DDPG (tau=0.01)')
plt.plot(moving_avg(ddpg_fail_returns), label='DDPG (tau=1.0 - Unstable)')
plt.xlabel('Episodes')
plt.ylabel('Smoothed Returns')
plt.title('DDPG Training on Pendulum')
plt.legend()
plt.show()

pendulum = Pendulum() print("Training DDPG on Pendulum (with target soft update)...") ddpg_returns = train_ddpg(pendulum, episodes=100, tau=0.01) print("Training DDPG on Pendulum (FAIL: without soft update, tau=1.0)...") ddpg_fail_returns = train_ddpg(pendulum, episodes=100, tau=1.0) plt.figure(figsize=(6, 4)) plt.plot(moving_avg(ddpg_returns), label='DDPG (tau=0.01)') plt.plot(moving_avg(ddpg_fail_returns), label='DDPG (tau=1.0 - Unstable)') plt.xlabel('Episodes') plt.ylabel('Smoothed Returns') plt.title('DDPG Training on Pendulum') plt.legend() plt.show()

Training DDPG on Pendulum (with target soft update)...

Training DDPG on Pendulum (FAIL: without soft update, tau=1.0)...

/tmp/ipykernel_243673/4247286845.py:71: RuntimeWarning: overflow encountered in matmul
  grad_x = grad_z1 @ self.W1.T
/tmp/ipykernel_243673/4247286845.py:52: RuntimeWarning: overflow encountered in matmul
  self.z2 = self.a1 @ self.W2 + self.b2
/tmp/ipykernel_243673/581904816.py:43: RuntimeWarning: invalid value encountered in subtract
  td_error = q_curr - targets
/tmp/ipykernel_243673/4247286845.py:58: RuntimeWarning: invalid value encountered in matmul
  grad_W2 = self.a1.T @ grad_z2
/tmp/ipykernel_243673/4247286845.py:60: RuntimeWarning: overflow encountered in matmul
  grad_a1 = grad_z2 @ self.W2.T
/tmp/ipykernel_243673/4247286845.py:61: RuntimeWarning: invalid value encountered in multiply
  grad_z1 = grad_a1 * (self.z1 > 0)

In [11]:

# Deterministic Assertions
test_sumtree = SumTree(4)
test_sumtree.add(1.0, "A")
test_sumtree.add(2.0, "B")
assert test_sumtree.total() == 3.0, "SumTree total should be 3.0"
idx, p, data = test_sumtree.get(1.5)
assert data == "B" and p == 2.0, "SumTree get failing to find correct interval"
print("Assertions passed!")

# Deterministic Assertions test_sumtree = SumTree(4) test_sumtree.add(1.0, "A") test_sumtree.add(2.0, "B") assert test_sumtree.total() == 3.0, "SumTree total should be 3.0" idx, p, data = test_sumtree.get(1.5) assert data == "B" and p == 2.0, "SumTree get failing to find correct interval" print("Assertions passed!")

Assertions passed!

7. Connections & Takeaways¶

• Double DQN fixes the overestimation bias of standard Q-learning by evaluating the greedy action with a separate target network. This makes value estimates much closer to the true expected return.
• Dueling DQN separates state-value $V$ from advantage $A$. In states where actions don't matter much, the network can learn $V(s)$ quickly and effectively generalize across actions.
• PER accelerates learning by focusing on transitions where the TD error is high (meaning the agent was "surprised"). The SumTree ensures we can sample proportionally without an $O(N)$ penalty per step.
• DDPG brings DQN ideas to continuous action spaces using an Actor-Critic architecture and deterministic policies. Soft updates to the target networks (tau) are absolutely critical to maintaining stability.