Weight Initialization Strategies

# Weight Initialization Strategies

## Introduction & Motivation

Weight Initialization Strategies dictate how trainable parameters in deep neural networks are assigned initial numerical values prior to gradient-based optimization. In deep learning architectures, proper initialization is not merely a convergence acceleration technique; it is a fundamental mathematical prerequisite for establishing stable signal propagation across deep computation graphs. Naive initialization methods, such as zero-initialization or unscaled random Gaussian initialization, cause catastrophic gradient degradation—manifesting as vanishing or exploding gradients—which prevents deep feedforward networks, convolutional networks, and transformers from training altogether.

The core motivation behind rigorous weight initialization is variance preservation across layers. As activations flow forward through matrix multiplications and non-linear activation functions, the variance of activation signals can exponentially decay to zero or explode to infinity depending on the spectral norm and variance of the weight matrices. Similarly, during backpropagation, the variance of error gradients flowing backward through the transpose of weight matrices follows identical exponential kinetics. Effective weight initialization strategies carefully configure the initial probability distribution of weights so that both activation variance during the forward pass and gradient variance during the backward pass remain bounded and stable across arbitrary network depth.

Modern deep learning frameworks depend on tailored initialization schemes designed for specific activation functions and structural topology. Key developments include Xavier (Glorot) initialization for symmetric linear or sigmoidal activations, Kaiming (He) initialization for rectified linear units (ReLU) and its variants, Layer-Sequential Unit-Variance (LSUV) initialization for data-driven empirical calibration, and Orthogonal initialization for recurrent structures and deep residual paths. Understanding the mathematical derivations, failure modes, and practical implementation trade-offs of these initialization strategies is essential for building robust, scalable neural networks.

---

## Core Concepts & Theory

### 1. The Physics of Signal Flow and Gradient Breakdown

Consider an $$L$$-layer deep neural network where each layer $$l \in \{1, 2, \dots, L\}$$ computes a linear transformation $$z^{(l)} = W^{(l)} a^{(l-1)} + b^{(l)}$$ followed by an element-wise activation $$a^{(l)} = g(z^{(l)})$$. If weights are initialized independently with mean zero and variance $$ ext{Var}(W)$$, the variance of linear outputs $$z^{(l)}$$ expands as a function of the layer input dimension $$n_{ ext{in}}^{(l)}$$.

When weights are initialized with excessive variance, the magnitude of activation vectors grows exponentially with network depth $$L$$. In saturated activation regimes (such as $$ anh$$ or sigmoid), large activation magnitudes push neurons into flat saturation regions where local derivatives approach zero, resulting in vanishing gradients. Conversely, when weights are initialized with insufficient variance, activation signals attenuate exponentially at each layer, reducing deep representations to zero and preventing useful feature learning.

### 2. Xavier (Glorot) Initialization

Introduced by Xavier Glorot and Yoshua Bengio (2010), Xavier initialization designs weight distribution parameters to equalize the variance of activations between input and output layers under the assumption of linear or symmetric zero-centered activations (such as $$ anh$$).

To preserve activation variance forward ($$ ext{Var}(a^{(l)}) = ext{Var}(a^{(l-1)})$$) and gradient variance backward ($$ ext{Var}(\frac{\partial \mathcal{L}}{\partial z^{(l-1)}}) = ext{Var}(\frac{\partial \mathcal{L}}{\partial z^{(l)}})$$), Xavier initialization balances the input fan ($$n_{ ext{in}}$$) and output fan ($$n_{ ext{out}}$$). The variance requirement yields:

$$ ext{Var}(W) = \frac{2}{n_{ ext{in}} + n_{ ext{out}}}$$

For a uniform distribution $$ ext{Uniform}(-a, a)$$, whose variance is $$\frac{a^2}{3}$$, setting $$\frac{a^2}{3} = \frac{2}{n_{ ext{in}} + n_{ ext{out}}}$$ gives the bounds $$a = \sqrt{\frac{6}{n_{ ext{in}} + n_{ ext{out}}}}$$.

### 3. Kaiming (He) Initialization

Kaiming He et al. (2015) identified that Xavier initialization fails in networks utilizing Rectified Linear Units (ReLU). Because ReLU sets all negative linear outputs to zero ($$g(z) = \max(0, z)$$), it zero-out exactly half of the activation space for symmetric zero-mean inputs. Consequently, the variance of activations is halved at every single layer:

$$ ext{Var}(a^{(l)}) = \frac{1}{2} ext{Var}(z^{(l)})$$

To compensate for this 50% variance reduction per layer, Kaiming initialization doubles the target weight variance relative to $$n_{ ext{in}}$$:

$$ ext{Var}(W) = \frac{2}{n_{ ext{in}}}$$

For He Normal initialization, weights are sampled from $$\mathcal{N}\left(0, \sqrt{\frac{2}{n_{ ext{in}}}} ight)$$. For He Uniform initialization, weights are sampled from $$ ext{Uniform}\left(-\sqrt{\frac{6}{n_{ ext{in}}}}, \sqrt{\frac{6}{n_{ ext{in}}}} ight)$$. For LeakyReLU activations with negative slope $$\alpha$$, the required variance generalizes to $$ ext{Var}(W) = \frac{2}{1 + \alpha^2}$$.

### 4. Variance Scaling and Initialization Failure Modes

Variance Scaling provides a unified framework generalizing Glorot and He initialization by parameterizing scale factor $$k$$, distribution type (normal, uniform, truncated normal), and fan mode ($$ ext{fan\_in}$$, $$ ext{fan\_out}$$, $$ ext{fan\_avg}$$).

However, standard variance scaling fails under specific conditions:
1. Unbounded Activation Drift: When non-linearities have non-zero mean outputs (e.g., standard Sigmoid or ELU without offset), activations accumulate positive bias, causing variance scaling assumptions to collapse.
2. Deep Residual Accumulation: In residual networks ($$a^{(l)} = a^{(l-1)} + F(a^{(l-1)}, W^{(l)})$$), adding the identity connection causes activation variance to grow linearly with depth $$L$$, leading to gradient explosion at $$L > 100$$ unless residual branches are scaled down (e.g., by $$\frac{1}{\sqrt{2L}}$$ or initialized near zero).
3. Mismatched Gain Factors: Applying Xavier initialization to deep ReLU networks causes activation variance to decay by $$\left(\frac{1}{2} ight)^L$$, causing signal extinction.

---

## Mathematical Formulation

### 1. Derivation of Forward Pass Variance Preservation

Let $$z_i^{(l)} = \sum_{j=1}^{n_{ ext{in}}} W_{ij}^{(l)} a_j^{(l-1)} + b_i^{(l)}$$. Assuming $$W_{ij}^{(l)}$$ and $$a_j^{(l-1)}$$ are independent random variables with zero mean:

$$\mathbb{E}[z_i^{(l)}] = \sum_{j=1}^{n_{ ext{in}}} \mathbb{E}[W_{ij}^{(l)}] \mathbb{E}[a_j^{(l-1)}] = 0$$

Computing the variance of $$z_i^{(l)}$$:

$$ ext{Var}(z_i^{(l)}) = \sum_{j=1}^{n_{ ext{in}}} ext{Var}\left(W_{ij}^{(l)} a_j^{(l-1)} ight)$$

Using the identity for independent zero-mean variables $$ ext{Var}(X Y) = ext{Var}(X) ext{Var}(Y)$$:

$$ ext{Var}(z_i^{(l)}) = n_{ ext{in}}^{(l)} ext{Var}\left(W^{(l)} ight) ext{Var}\left(a^{(l-1)} ight)$$

For linear/symmetric activations where $$ ext{Var}(a^{(l-1)}) = ext{Var}(z^{(l-1)})$$, maintaining constant variance across layers $$ ext{Var}(z^{(l)}) = ext{Var}(z^{(l-1)})$$ requires:

$$n_{ ext{in}}^{(l)} ext{Var}\left(W^{(l)} ight) = 1 \implies ext{Var}\left(W^{(l)} ight) = \frac{1}{n_{ ext{in}}^{(l)}}$$

### 2. Derivation of Backward Pass Gradient Variance Preservation

During backpropagation, the gradient of loss $$\mathcal{L}$$ with respect to linear inputs $$z_j^{(l-1)}$$ is:

$$\frac{\partial \mathcal{L}}{\partial z_j^{(l-1)}} = \left( \sum_{i=1}^{n_{ ext{out}}} \frac{\partial \mathcal{L}}{\partial z_i^{(l)}} W_{ij}^{(l)} ight) g'\left(z_j^{(l-1)} ight)$$

Assuming derivative $$g'(z) \approx 1$$ on average for normalized linear regimes:

$$ ext{Var}\left(\frac{\partial \mathcal{L}}{\partial z_j^{(l-1)}} ight) = n_{ ext{out}}^{(l)} ext{Var}\left(W^{(l)} ight) ext{Var}\left(\frac{\partial \mathcal{L}}{\partial z_i^{(l)}} ight)$$

Preserving gradient variance backward requires:

$$n_{ ext{out}}^{(l)} ext{Var}\left(W^{(l)} ight) = 1 \implies ext{Var}\left(W^{(l)} ight) = \frac{1}{n_{ ext{out}}^{(l)}}$$

Combining the forward condition ($$\frac{1}{n_{ ext{in}}}$$) and backward condition ($$\frac{1}{n_{ ext{out}}}$$) by harmonic mean yields the Xavier Glorot variance target:

$$ ext{Var}(W) = \frac{2}{n_{ ext{in}} + n_{ ext{out}}}$$

### 3. Derivation of Kaiming He Rectified Variance

For ReLU activations $$a = \max(0, z)$$, assuming $$z \sim \mathcal{N}(0, \sigma_z^2)$$:

$$\mathbb{E}[a] = \int_{0}^{\infty} z \frac{1}{\sqrt{2\pi}\sigma_z} e^{-\frac{z^2}{2\sigma_z^2}} dz = \frac{\sigma_z}{\sqrt{2\pi}}$$

$$\mathbb{E}[a^2] = \int_{0}^{\infty} z^2 \frac{1}{\sqrt{2\pi}\sigma_z} e^{-\frac{z^2}{2\sigma_z^2}} dz = \frac{1}{2} \sigma_z^2$$

Thus, the variance of $$a$$ is:

$$ ext{Var}(a) = \mathbb{E}[a^2] - (\mathbb{E}[a])^2 = \left(\frac{1}{2} - \frac{1}{2\pi} ight) \sigma_z^2$$

Since $$ ext{Var}(z^{(l)}) = n_{ ext{in}} ext{Var}(W^{(l)}) \mathbb{E}[(a^{(l-1)})^2] = n_{ ext{in}} ext{Var}(W^{(l)}) \frac{1}{2} ext{Var}(z^{(l-1)})$$, forcing $$ ext{Var}(z^{(l)}) = ext{Var}(z^{(l-1)})$$ yields:

$$n_{ ext{in}} ext{Var}(W^{(l)}) \frac{1}{2} = 1 \implies ext{Var}(W^{(l)}) = \frac{2}{n_{ ext{in}}}$$

---

## Advanced Theory & Extensions

### 1. Orthogonal Initialization and Isotropic Singular Values

While Gaussian variance scaling preserves variance on average, random Gaussian matrices possess singular values distributed according to the Marchenko-Pastur law. This spectral dispersion means some direction vectors are stretched while others are compressed, leading to gradual norm distortion in very deep networks ($$L > 50$$).

Orthogonal initialization computes an orthogonal matrix $$Q$$ via QR decomposition of a random Gaussian matrix:

$$W = Q \quad ext{where} \quad A = Q R, \quad A_{ij} \sim \mathcal{N}(0, 1)$$

Since orthogonal matrices satisfy $$Q^T Q = I$$, all singular values of $$W$$ are exactly equal to 1. This guarantees norm-preserving dynamical isometry, allowing gradients to propagate unattenuated through hundreds of layers without batch normalization.

### 2. Layer-Sequential Unit-Variance (LSUV) Empirical Calibration

LSUV (Mishkin & Matas, 2015) combines orthogonal initialization with iterative empirical data-driven scaling:
1. Initialize weights of convolutional/dense layers with orthogonal matrices.
2. Pass a mini-batch of real training data through the network layer by layer.
3. Compute empirical variance $$ ext{Var}(a^{(l)})$$ of activations at layer $$l$$.
4. Scale weights $$W^{(l)} \leftarrow \frac{W^{(l)}}{\sqrt{ ext{Var}(a^{(l)})}}$$ iteratively until $$| ext{Var}(a^{(l)}) - 1.0| < \epsilon$$.

LSUV handles non-standard activation functions, arbitrary complex topologies, and pre-activation structures where analytical variance formulas are difficult to derive.

### 3. ReZero and Fixup Initialization for Deep Residual Networks

In ultra-deep Residual Networks (ResNets), standard Kaiming initialization causes activation variance at block $$L$$ to scale as $$\mathcal{O}(L)$$, requiring Batch Normalization to prevent explosion. ReZero and Fixup remove normalization layers by modifying initialization:
- Fixup: Scales the weights of the final layer in each residual branch to zero ($$W_{ ext{final}} = 0$$) or by $$\frac{1}{L^{1/(2m)}}$$.
- ReZero: Introduces a learnable scalar parameter $$\alpha_i$$ initialized to zero for each residual branch: $$a^{(l)} = a^{(l-1)} + \alpha_l F(a^{(l-1)})$$. At step 0, the network functions as an exact identity map with identity signal propagation.

---

## Computational Considerations

---

## Practical Implementation Strategies

### 1. Per-Layer Activation-Aware Initialization Protocol

1. Convolutional & Dense Layers with ReLU/GELU/SiLU: Always select Kaiming (He) Normal initialization with $$ ext{fan\_in}$$ scaling mode.
2. Convolutional & Dense Layers with Tanh/Sigmoid: Select Xavier (Glorot) Normal or Uniform initialization.
3. Recurrent Layers (LSTM/GRU gates): Initialize recurrent weight matrices $$W_{hh}$$ with Orthogonal initialization and input weights $$W_{xh}$$ with Kaiming Normal.
4. Bias Vector Initialization: Initialize bias vectors to constant zero ($$b = 0$$). Exceptions:
- LSTM Forget Gate Bias: Initialize to $$+1.0$$ or $$+2.0$$ to prevent early memory loss.
- Output Class Imbalance: Initialize output layer bias to $$\log(\frac{p}{1-p})$$ where $$p$$ is prior class probability.

### 2. Gain Factor Scaling for Non-Linearities

Modern framework implementations (such as PyTorch `torch.nn.init.calculate_gain`) incorporate activation gain $$g$$:

$$W \sim \mathcal{N}\left(0, \frac{g^2}{ ext{fan\_in}} ight)$$

• Linear / Identity: $$g = 1.0$$
• Tanh: $$g = \frac{5}{3} \approx 1.667$$
• ReLU: $$g = \sqrt{2} \approx 1.414$$
• LeakyReLU(\alpha): $$g = \sqrt{\frac{2}{1 + \alpha^2}}$$
• GELU / SiLU: $$g \approx 1.414$$

---

## Benchmark Datasets & Evaluation

Initialization performance is evaluated across standard benchmark datasets by measuring convergence rate, initial loss value, gradient norm stability across depth, and final classification accuracy:

• CIFAR-10 / CIFAR-100: Evaluating 20-layer to 110-layer ResNets without Batch Normalization under Xavier vs. He vs. Fixup initialization.
• ImageNet (ILSVRC): Assessing training stability during early epochs of deep Vision Transformers (ViT) and ConvNeXt architectures.
• MNIST / Fashion-MNIST: Diagnostic benchmark for testing extreme depth (e.g., 50-layer MLP) signal propagation without normalizers.

---

## Key Challenges & Limitations

### 1. Activation Function Mismatch and Variance Collapse

Applying Xavier initialization ($$g=1.0$$) to a 30-layer ReLU network causes activation variance at layer 30 to collapse by a factor of $$\left(\frac{1}{2} ight)^{30} \approx 9.3 imes 10^{-10}$$. The network fails to train completely because weights in deeper layers receive effectively zero gradient updates.

### 2. Variance Scaling Failure in Residual Networks

In architectures with skip connections ($$y = x + f(x)$$), standard variance scaling assumes $$f(x)$$ operates independently. Adding $$x$$ increases activation variance at each residual block:

$$ ext{Var}(y) = ext{Var}(x) + ext{Var}(f(x))$$

For a 100-block ResNet, total variance grows to $$100 imes ext{Var}(x)$$, inducing gradient explosion unless residual branch output weights are explicitly downscaled.

### 3. Batch Normalization Dependency Masking

The widespread adoption of Batch Normalization (BN) and Layer Normalization (LN) can mask poor weight initialization because BN re-scales activations to unit variance at every layer. However, poor initialization with BN still impairs optimization by causing ill-conditioned initial gradients and slower convergence.

---

## Hyperparameter Tuning

• Fan Mode: Choose $$ ext{fan\_in}$$ to preserve forward activation variance; choose $$ ext{fan\_out}$$ to preserve backward gradient variance; choose $$ ext{fan\_avg}$$ for balanced compromise (Glorot default).
• Distribution Type: Normal distribution avoids sharp boundary truncations; Truncated Normal (within $$\pm 2\sigma$$) avoids rare extreme outlier values that provoke early gradient spikes.
• LeakyReLU Negative Slope (\alpha): Must adjust Kaiming gain $$g = \sqrt{\frac{2}{1+\alpha^2}}$$. For $$\alpha = 0.2$$, $$g = \sqrt{\frac{2}{1.04}} \approx 1.386$$.

---

## Real-World Applications & Case Studies

1. Ultra-Deep Computer Vision (ResNet-101 / ConvNeXt): Utilizing Kaiming Normal initialization paired with zero-initialized final $$1 imes 1$$ conv blocks in residual branches ensures rapid, stable convergence on multi-GPU clusters.
2. Large Language Models & Transformers (GPT-4 / LLaMA): Transformer linear layers utilize scaled Xavier initialization where weights are scaled by $$\frac{1}{\sqrt{2 N_{ ext{layers}}}}$$ for residual projections, ensuring attention logit stability.
3. Deep Reinforcement Learning (Actor-Critic Networks): Policy head output weights are initialized with very small uniform bounds (e.g., $$ ext{Uniform}(-3 imes 10^{-3}, 3 imes 10^{-3})$$) to ensure initial exploration policy is near-uniform across action spaces.

---

## Integration with Other Methods

Weight Initialization operates synergistically with optimizer dynamics, learning rate warm-up, and normalization layers:
- Initialization + Learning Rate Warm-up: Linear warm-up during early steps prevents large Adam optimizer second-moment momentum vectors ($$v_t$$) from destroying carefully initialized weight representations.
- Initialization + Weight Decay (L2 Regularization): Proper initial weight scale prevents early weight decay steps from shrinking parameters into zero-activation collapse.

---

## Summary & Key Takeaways

1. Purpose: Weight initialization controls signal propagation and prevents vanishing/exploding gradients in deep networks.
2. Xavier (Glorot): Designed for linear/symmetric activations (Tanh), target variance $$ ext{Var}(W) = \frac{2}{n_{ ext{in}} + n_{ ext{out}}}$$.
3. Kaiming (He): Designed for ReLU/GELU activations, compensates for 50% activation death with target variance $$ ext{Var}(W) = \frac{2}{n_{ ext{in}}}$$.
4. Orthogonal: Preserves dynamical isometry and singular values, critical for RNNs and ultra-deep structures.
5. Residual Networks: Require branch scaling ($$rac{1}{\sqrt{2L}}$$ or zero-init) to prevent linear variance accumulation across depth.

---

## Appendix: Practical Labs

### Lab 1: Forward Activation Variance Propagation Analysis

This laboratory script simulates forward activation variance across a 30-layer deep Multi-Layer Perceptron (MLP) under three initialization strategies: Naive Gaussian, Xavier Uniform, and Kaiming Normal with ReLU activations.

import numpy as np

def run_forward_variance_simulation():
    np.random.seed(42)
    num_layers = 30
    layer_dim = 256
    batch_size = 128

    # Input data with unit variance
    x = np.random.normal(0, 1, size=(batch_size, layer_dim))

    # Test 1: Naive Small Normal Initialization (std = 0.01)
    activations_naive = [x]
    a = x
    for _ in range(num_layers):
        w = np.random.normal(0, 0.01, size=(layer_dim, layer_dim))
        z = a @ w
        a = np.maximum(0, z) # ReLU
        activations_naive.append(a)

    var_naive_end = np.var(activations_naive[-1])

    # Test 2: Xavier Uniform Initialization
    activations_xavier = [x]
    a = x
    limit = np.sqrt(6.0 / (layer_dim + layer_dim))
    for _ in range(num_layers):
        w = np.random.uniform(-limit, limit, size=(layer_dim, layer_dim))
        z = a @ w
        a = np.maximum(0, z) # ReLU
        activations_xavier.append(a)

    var_xavier_end = np.var(activations_xavier[-1])

    # Test 3: Kaiming (He) Normal Initialization
    activations_he = [x]
    a = x
    std_he = np.sqrt(2.0 / layer_dim)
    for _ in range(num_layers):
        w = np.random.normal(0, std_he, size=(layer_dim, layer_dim))
        z = a @ w
        a = np.maximum(0, z) # ReLU
        activations_he.append(a)

    var_he_end = np.var(activations_he[-1])

    assert var_naive_end < 1e-10, "Naive initialization must suffer activation collapse"
    assert var_xavier_end < var_he_end, "Xavier must produce smaller variance than He under ReLU"
    assert 0.1 < var_he_end < 10.0, "He initialization must maintain stable non-zero variance"

    print("Lab 1 Output:")
    print("Naive Init Final Layer Variance:", f"{var_naive_end:.4e}")
    print("Xavier Init Final Layer Variance:", f"{var_xavier_end:.4e}")
    print("Kaiming He Final Layer Variance:", f"{var_he_end:.4f}")

run_forward_variance_simulation()

### Lab 2: Backward Gradient Variance Tracking

This laboratory script tracks backward gradient norms across 25 layers during backpropagation to compare gradient flow stability under Xavier vs Kaiming He initialization.

import numpy as np

def run_backward_gradient_simulation():
    np.random.seed(101256)
    num_layers = 25
    dim = 200
    batch_size = 64

    # Simulate weights and forward passes under He Normal
    weights_he = []
    activations_he = []
    masks_he = []

    a = np.random.normal(0, 1, size=(batch_size, dim))
    activations_he.append(a)

    std_he = np.sqrt(2.0 / dim)
    for _ in range(num_layers):
        w = np.random.normal(0, std_he, size=(dim, dim))
        z = a @ w
        mask = (z > 0).astype(float)
        a = z * mask
        weights_he.append(w)
        masks_he.append(mask)
        activations_he.append(a)

    # Backpropagate gradient from loss
    grad = np.random.normal(0, 1, size=(batch_size, dim))
    grad_norms_he = []

    for l in reversed(range(num_layers)):
        grad = (grad * masks_he[l]) @ weights_he[l].T
        grad_norms_he.append(np.linalg.norm(grad))

    grad_norms_he.reverse()

    first_layer_grad = grad_norms_he[0]
    last_layer_grad = grad_norms_he[-1]
    ratio = first_layer_grad / last_layer_grad

    assert np.isfinite(first_layer_grad), "Gradients must remain finite"
    assert 0.01 < ratio < 100.0, "Gradient ratio across 25 layers must remain bounded"

    print("
Lab 2 Output:")
    print("Layer 25 Grad Norm:", round(last_layer_grad, 2))
    print("Layer 1 Grad Norm:", round(first_layer_grad, 2))
    print("Backward Gradient Stability Ratio:", round(ratio, 4))

run_backward_gradient_simulation()

### Lab 3: Orthogonal Matrix QR Initialization and Spectral Property Verification

This laboratory script constructs an orthogonal initialization using QR decomposition and verifies that all singular values of the weight matrix are identically equal to 1.0.

import numpy as np

def verify_orthogonal_initialization():
    np.random.seed(42)
    fan_in = 150
    fan_out = 150

    # Generate Gaussian matrix and perform QR decomposition
    gaussian_matrix = np.random.normal(0, 1, size=(fan_in, fan_out))
    q_matrix, r_matrix = np.linalg.qr(gaussian_matrix)

    # Correct for sign flips in R to ensure uniform Haar distribution
    d = np.diag(r_matrix)
    ph = np.sign(d)
    q_matrix = q_matrix * ph

    # Compute singular values via SVD
    singular_values = np.linalg.svd(q_matrix, compute_uv=False)
    
    # Verify orthogonality: Q^T * Q = I
    identity_check = q_matrix.T @ q_matrix
    eye = np.eye(fan_in)
    max_ortho_error = np.max(np.abs(identity_check - eye))

    assert max_ortho_error < 1e-12, "Orthogonal matrix must satisfy Q^T Q = I"
    assert np.allclose(singular_values, 1.0, atol=1e-12), "All singular values must equal 1.0"

    print("
Lab 3 Output:")
    print("Max Orthogonality Error (Q^T Q - I):", f"{max_ortho_error:.4e}")
    print("Min Singular Value:", round(float(np.min(singular_values)), 6))
    print("Max Singular Value:", round(float(np.max(singular_values)), 6))

verify_orthogonal_initialization()

### Lab 4: Data-Driven LSUV Empirical Calibration

This laboratory script implements Layer-Sequential Unit-Variance (LSUV) empirical calibration on a 5-layer network with custom non-linear activations to bring layer variance to exactly 1.0.

import numpy as np

class SimpleDenseLayer:
    def __init__(self, in_dim, out_dim):
        # Initial random matrix
        a = np.random.normal(0, 1, size=(in_dim, out_dim))
        q, _ = np.linalg.qr(a)
        self.w = q[:in_dim, :out_dim]

    def forward(self, x):
        # Custom activation: LeakyReLU with alpha = 0.1
        z = x @ self.w
        return np.where(z > 0, z, 0.1 * z)

def run_lsuv_calibration():
    np.random.seed(42)
    batch_size = 256
    dim = 128
    num_layers = 5

    calibration_data = np.random.normal(0, 1, size=(batch_size, dim))

    layers = [SimpleDenseLayer(dim, dim) for _ in range(num_layers)]

    print("
Lab 4 Output:")
    # Perform LSUV empirical calibration layer by layer
    for idx, layer in enumerate(layers):
        # Pass calibration data up to current layer
        h = calibration_data
        for prev_layer in layers[:idx]:
            h = prev_layer.forward(h)

        # Iteratively scale weights until variance is 1.0
        for iteration in range(10):
            output = layer.forward(h)
            current_var = np.var(output)
            if abs(current_var - 1.0) < 1e-4:
                break
            layer.w = layer.w / np.sqrt(current_var)

        final_var = np.var(layer.forward(h))
        assert abs(final_var - 1.0) < 1e-2, f"Layer {idx+1} LSUV variance calibration failed"
        print(f"Layer {idx+1} Calibrated Activation Variance: {final_var:.4f}")

run_lsuv_calibration()