1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

import numpy as np
 
def rmsprop_step(theta, gradient, v, alpha=0.001, rho=0.9, epsilon=1e-8):
    """
    Single RMSprop update step.
    
    Args:
        theta: Parameter vector
        gradient: Current gradient
        v: Running EMA of squared gradients (modified in-place)
        alpha: Base learning rate
        rho: Decay rate for EMA (typically 0.9)
        epsilon: Numerical stability constant
        
    Returns:
        Updated theta
    """
    # Update EMA of squared gradients
    v[:] = rho * v + (1 - rho) * gradient ** 2
    
    # Compute adapted update
    # Large v[i] -> small step for parameter i
    # Small v[i] -> large step for parameter i
    update = alpha / (np.sqrt(v) + epsilon) * gradient
    
    return theta - update
 
 
# Demonstration: RMSprop vs AdaGrad over long training
def compare_accumulator_growth(steps=10000):
    """Compare how accumulator/EMA evolves over time."""
    gradient_scale = 1.0
    
    # AdaGrad accumulator
    G_adagrad = 0.0
    adagrad_lrs = []
    
    # RMSprop EMA
    v_rmsprop = 0.0
    rho = 0.9
    rmsprop_lrs = []
    
    alpha = 1.0  # Base LR for comparison
    
    for t in range(1, steps + 1):
        g = np.random.randn() * gradient_scale
        
        # AdaGrad: accumulate
        G_adagrad += g ** 2
        adagrad_lrs.append(alpha / (np.sqrt(G_adagrad) + 1e-8))
        
        # RMSprop: EMA
        v_rmsprop = rho * v_rmsprop + (1 - rho) * g ** 2
        rmsprop_lrs.append(alpha / (np.sqrt(v_rmsprop) + 1e-8))
    
    return adagrad_lrs, rmsprop_lrs
 
adagrad_lrs, rmsprop_lrs = compare_accumulator_growth()
 
print("Effective learning rates over 10,000 steps:")
print(f"                  Step 100    Step 1000   Step 10000")
print(f"AdaGrad:          {adagrad_lrs[99]:.6f}  {adagrad_lrs[999]:.6f}  {adagrad_lrs[9999]:.6f}")
print(f"RMSprop (ρ=0.9):  {rmsprop_lrs[99]:.6f}  {rmsprop_lrs[999]:.6f}  {rmsprop_lrs[9999]:.6f}")
 
# RMSprop maintains roughly constant effective LR
# AdaGrad's effective LR decays as O(1/sqrt(t))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

import numpy as np
 
def demonstrate_ema_properties():
    """Visualize key EMA properties."""
    
    # 1. Bias toward zero (cold start problem)
    print("=== Cold Start Bias ===")
    x = 1.0  # Constant input
    rho = 0.9
    ema = 0.0
    
    for t in range(1, 21):
        ema = rho * ema + (1 - rho) * x
        print(f"Step {t:2d}: EMA = {ema:.4f}, True mean = 1.0, Bias = {1.0 - ema:.4f}")
    
    # Takes ~23 steps to get within 10% of true value!
    # This motivates "bias correction" (used in Adam)
    
    # 2. Effective window demonstration
    print("
=== Effective Window ===")
    print("Weight of past values in EMA (ρ=0.9):")
    for k in range(10):
        weight = (1 - 0.9) * (0.9 ** k)
        cumulative = 1 - (0.9 ** (k + 1))
        print(f"  x[t-{k}]: weight = {weight:.4f}, cumulative = {cumulative:.4f}")
    
    # 3. Response to sudden change
    print("
=== Response to Sudden Change ===")
    sequence = [1.0] * 50 + [10.0] * 50  # Jump from 1 to 10 at t=50
    rho = 0.9
    ema = 0.0
    
    for t, x in enumerate(sequence):
        ema = rho * ema + (1 - rho) * x
        if t in [48, 49, 50, 51, 55, 60, 70, 99]:
            print(f"Step {t}: x = {x:.1f}, EMA = {ema:.2f}")
    
    # EMA takes ~10 steps to mostly catch up to the new value
 
demonstrate_ema_properties()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

import numpy as np
 
def centered_rmsprop_step(theta, gradient, v, g_bar, 
                          alpha=0.001, rho=0.99, epsilon=1e-8):
    """
    Centered RMSprop: normalizes by variance instead of second moment.
    
    Args:
        theta: Parameters
        gradient: Current gradient
        v: EMA of squared gradients (E[g²])
        g_bar: EMA of gradients (E[g])
        alpha: Learning rate
        rho: Decay rate
        epsilon: Stability constant
        
    Returns:
        Updated theta
    """
    # Update EMAs
    v[:] = rho * v + (1 - rho) * gradient ** 2
    g_bar[:] = rho * g_bar + (1 - rho) * gradient
    
    # Variance estimate: Var[g] = E[g²] - E[g]²
    variance = v - g_bar ** 2
    
    # Ensure non-negative (numerical issues can make this slightly negative)
    variance = np.maximum(variance, 0)
    
    # Update with variance normalization
    update = alpha / (np.sqrt(variance) + epsilon) * gradient
    
    return theta - update
 
 
# Compare standard vs centered on biased gradients
def compare_variants():
    np.random.seed(42)
    theta_std = 0.0
    theta_ctr = 0.0
    
    v_std = 0.0
    v_ctr = 0.0
    g_bar = 0.0
    
    alpha = 0.1
    rho = 0.9
    
    # Gradient with large mean, small variance
    # Mean = 5, Std = 0.5
    for step in range(100):
        g = 5.0 + np.random.randn() * 0.5
        
        # Standard RMSprop
        v_std = rho * v_std + (1 - rho) * g ** 2
        update_std = alpha / (np.sqrt(v_std) + 1e-8) * g
        theta_std -= update_std
        
        # Centered RMSprop
        v_ctr = rho * v_ctr + (1 - rho) * g ** 2
        g_bar = rho * g_bar + (1 - rho) * g
        var = max(0, v_ctr - g_bar ** 2)
        update_ctr = alpha / (np.sqrt(var) + 1e-8) * g
        theta_ctr -= update_ctr
    
    print(f"After 100 steps with biased gradient (mean=5, std=0.5):")
    print(f"  Standard RMSprop:  θ = {theta_std:.2f}")
    print(f"  Centered RMSprop:  θ = {theta_ctr:.2f}")
    print(f"  Centered moved further (larger effective LR for consistent gradients)")
 
compare_variants()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151

import numpy as np
from typing import Dict, Optional
 
class RMSprop:
    """
    RMSprop (Root Mean Square Propagation) optimizer.
    
    Maintains exponential moving average of squared gradients for
    adaptive per-parameter learning rate scaling.
    
    Update rule:
        v_t = ρ * v_{t-1} + (1-ρ) * g_t²
        θ_t = θ_{t-1} - α / (√v_t + ε) * g_t
    
    Key differences from AdaGrad:
        - Uses EMA instead of sum (bounded accumulator)
        - Learning rate can increase when gradients shrink
        - Suitable for long training runs
    
    Reference:
        Hinton, Srivastava, Swersky. "Neural Networks for Machine 
        Learning" Lecture 6e (2012, unpublished).
    """
    
    def __init__(
        self,
        learning_rate: float = 0.01,
        rho: float = 0.9,
        epsilon: float = 1e-8,
        weight_decay: float = 0.0,
        momentum: float = 0.0,
        centered: bool = False,
    ):
        """
        Args:
            learning_rate: Base learning rate α
            rho: Decay rate for squared gradient EMA (typically 0.9 or 0.99)
            epsilon: Small constant for numerical stability
            weight_decay: L2 regularization coefficient
            momentum: Optional momentum term (not in original RMSprop)
            centered: If True, use centered variant (normalize by variance)
        """
        if not 0.0 < rho < 1.0:
            raise ValueError(f"rho must be in (0, 1), got {rho}")
        
        self.lr = learning_rate
        self.rho = rho
        self.eps = epsilon
        self.weight_decay = weight_decay
        self.momentum = momentum
        self.centered = centered
        
        # State
        self.v: Dict[str, np.ndarray] = {}  # EMA of squared gradients
        self.g_avg: Dict[str, np.ndarray] = {}  # EMA of gradients (centered)
        self.momentum_buffer: Dict[str, np.ndarray] = {}  # Momentum buffer
        self.steps = 0
    
    def step(
        self,
        parameters: Dict[str, np.ndarray],
        gradients: Dict[str, np.ndarray],
    ) -> Dict[str, np.ndarray]:
        """
        Perform one RMSprop step.
        
        Args:
            parameters: Dict mapping names to parameter arrays
            gradients: Dict mapping names to gradient arrays
            
        Returns:
            Updated parameters dict
        """
        self.steps += 1
        updated = {}
        
        for name, param in parameters.items():
            if name not in gradients:
                updated[name] = param
                continue
            
            grad = gradients[name].copy()
            
            # Apply weight decay
            if self.weight_decay > 0:
                grad = grad + self.weight_decay * param
            
            # Initialize state
            if name not in self.v:
                self.v[name] = np.zeros_like(param)
                if self.centered:
                    self.g_avg[name] = np.zeros_like(param)
                if self.momentum > 0:
                    self.momentum_buffer[name] = np.zeros_like(param)
            
            v = self.v[name]
            
            # Update EMA of squared gradients
            v[:] = self.rho * v + (1 - self.rho) * grad ** 2
            
            if self.centered:
                # Centered variant: normalize by variance
                g_avg = self.g_avg[name]
                g_avg[:] = self.rho * g_avg + (1 - self.rho) * grad
                avg = v - g_avg ** 2
                avg = np.maximum(avg, 0)  # Ensure non-negative
            else:
                avg = v
            
            # Compute update
            update = grad / (np.sqrt(avg) + self.eps)
            
            # Optional: apply momentum
            if self.momentum > 0:
                buf = self.momentum_buffer[name]
                buf[:] = self.momentum * buf + update
                update = buf
            
            updated[name] = param - self.lr * update
        
        return updated
    
    def get_state_summary(self) -> Dict[str, float]:
        """Get summary statistics for monitoring."""
        all_v = np.concatenate([v.flatten() for v in self.v.values()])
        all_lr = self.lr / (np.sqrt(all_v) + self.eps)
        return {
            'v_mean': float(np.mean(all_v)),
            'v_max': float(np.max(all_v)),
            'effective_lr_mean': float(np.mean(all_lr)),
            'effective_lr_min': float(np.min(all_lr)),
            'effective_lr_max': float(np.max(all_lr)),
        }
 
 
# Demonstration
if __name__ == "__main__":
    optimizer = RMSprop(learning_rate=0.01, rho=0.9)
    
    # Simple optimization: minimize (x - 3)² + (y - 2)²
    params = {"theta": np.array([0.0, 0.0])}
    target = np.array([3.0, 2.0])
    
    for step in range(200):
        grad = 2 * (params["theta"] - target)
        params = optimizer.step(params, {"theta": grad})
        
        if step % 50 == 49:
            loss = np.sum((params["theta"] - target) ** 2)
            stats = optimizer.get_state_summary()
            print(f"Step {step+1}: loss = {loss:.6f}, eff_lr = {stats['effective_lr_mean']:.6f}")