SmoothQuant: Mathematical Foundations, Per-Channel Outlier Migration, and Hardware-Efficient W8A8 Inference in Large Language Models

Serving large language models (LLMs) in production environments presents two distinct hardware bottlenecks. During the autoregressive generation (decode) phase with small batch sizes, inference is memory-bandwidth bound, as billions of parameters must be streamed from High Bandwidth Memory (HBM) to on-chip SRAM for every generated token. Conversely, during the prompt processing (prefill) phase and during high-concurrency batched serving, inference is compute bound, limited by the raw floating-point operations per second (FLOPS) of tensor execution units.

While weight-only quantization methods (such as GPTQ and AWQ) compress model storage and reduce memory bus traffic by storing weights in INT4 or INT8, they require dequantizing weights back to 16-bit floating-point (FP16 or BF16) before matrix multiplication. Consequently, weight-only schemes do not accelerate compute-bound workloads, do not reduce activation or Key-Value (KV) cache memory footprints, and fail to leverage integer arithmetic units on modern GPUs.

Achieving true 8-bit weight and 8-bit activation (W8A8) inference enables the direct use of hardware INT8 Tensor Core General Matrix Multiplications (GEMMs), offering up to double the computational throughput and halving memory usage compared to FP16. However, standard W8A8 post-training quantization has historically triggered severe accuracy degradation in models exceeding 6.7 billion parameters.

To resolve this limitation, researchers from MIT and NVIDIA developed SmoothQuant, a training-free post-training quantization (PTQ) framework. SmoothQuant introduces a mathematically equivalent transformation that migrates quantization difficulty from activations to weights offline, enabling lossless W8A8 quantization across frontier open architectures including OPT-175B, BLOOM-176B, GLM-130B, MT-NLG 530B, and Llama models.

The Activation Outlier Problem in Scaled Transformers

Understanding why standard 8-bit quantization fails in large language models requires examining the distribution of activations across model layers.

Distributional Asymmetry: Flat Weights vs. Spiky Activations

In standard Transformer architectures, the weights $\mathbf{W} \in \mathbb{R}^{C_i \times C_o}$ exhibit a relatively uniform and bounded numerical distribution across input channels $C_i$ and output channels $C_o$. The maximum absolute weight values rarely exceed moderate bounds, making weights straightforward to quantize to INT8 or INT4 with negligible quantization error.

In contrast, as demonstrated by Dettmers et al. in LLM.int8(), scaling language models beyond 6.7 billion parameters leads to the emergent phenomenon of systematic activation outliers. These outliers exhibit specific characteristics:

1. Extreme Magnitude: Outlier activations reach magnitudes up to 100 times larger than the median activation values (often exceeding $|X| > 100$ while typical features remain below $|X| < 1$).
2. Channel Localization: Outliers do not appear randomly across all dimensions; rather, they concentrate persistently within a tiny fraction (less than 0.1% to 1%) of the hidden feature channels across nearly all tokens.
3. Low Intra-Channel Variance: While variance across different channels for a single token is massive, the variance within an individual outlier channel across different sequence tokens is remarkably small.

Quantization Granularities and Hardware Constraints

Standard uniform symmetric $N$-bit quantization maps a real-valued tensor $\mathbf{X}$ to an integer grid:

$$\bar{\mathbf{X}}^{\text{INT8}} = \text{clamp}\left( \left\lfloor \frac{\mathbf{X}}{\Delta} \right\rceil, -2^{N-1}, 2^{N-1}-1 \right)$$

where the quantization step size $\Delta$ is governed by the maximum absolute value within the quantization group:

$$\Delta = \frac{\max(|\mathbf{X}|)}{2^{N-1} - 1}$$

When performing per-tensor or per-token activation quantization, the extreme magnitude $m = \max(|\mathbf{X}|)$ of the outlier channel dictates the step size $\Delta$. For a non-outlier channel $j$ with maximum value $m_j \ll m$, its effective dynamic range is compressed into an extremely narrow band of discrete integer bins:

$$\text{Effective Levels}_j \approx 2^N \cdot \frac{m_j}{m}$$

When $m_j / m \approx 0.01$, non-outlier features are assigned only 2 to 3 discrete integer values, discarding crucial representational information and causing catastrophic perplexity spikes.

Per-Tensor Quantization Dynamic Range Failure:

Channel 1 (Normal):  [-0.8, +0.6]  --> Quantized into ~2 integer bins (Massive Rounding Error)
Channel 2 (Normal):  [-0.4, +0.5]  --> Quantized into ~1 integer bin  (Information Destroyed)
Channel 3 (Outlier): [-98.0, +95.0] --> Sets global step size Delta = 98.0 / 127 = 0.7717

Why Per-Channel Activation Quantization Fails on Hardware

Mathematically, performing per-channel activation quantization (assigning an independent step size $\Delta_{X, j}$ to each input channel $j$) resolves the precision collapse. However, modern GPU tensor accelerators (such as NVIDIA Tensor Cores) execute matrix multiplication as a hardware-pipelined outer-product reduction:

$$\mathbf{Y} = \mathbf{X} \mathbf{W}, \quad \mathbf{X} \in \mathbb{R}^{T \times C_i}, \quad \mathbf{W} \in \mathbb{R}^{C_i \times C_o}$$

In hardware INT8 GEMM micro-architectures (e.g., CUTLASS / TensorRT-LLM), scaling factors can only be applied along the outer dimensions: the sequence dimension $T$ for activations (per-token scaling) and the output channel dimension $C_o$ for weights (per-channel weight scaling):

$$\mathbf{Y} = \text{diag}(\mathbf{\Delta}_X^{\text{FP16}}) \cdot \left( \bar{\mathbf{X}}^{\text{INT8}} \cdot \bar{\mathbf{W}}^{\text{INT8}} \right) \cdot \text{diag}(\mathbf{\Delta}_W^{\text{FP16}})$$

Applying per-channel scaling along the inner reduction dimension $C_i$ would require inserting non-uniform scaling operations inside the high-throughput systolic multiply-accumulate (MMA) pipeline, which hardware execution units do not support.

Earlier attempts to bypass this issue, such as LLM.int8(), extracted outlier channels into a separate FP16 matrix multiplication while computing the remaining 99.9% in INT8. However, this dynamic decomposition introduces branch divergent memory copies and kernel dispatch overheads, resulting in a 20% to 30% latency slowdown relative to native FP16 execution.

Mathematical Foundations of SmoothQuant

SmoothQuant circumvents the hardware limitation by recognizing that linear transformations allow exact algebraic factorizations across intermediate matrix boundaries.

Equivalent Linear Transformation

Consider a linear layer computing $\mathbf{Y} = \mathbf{X} \mathbf{W}$. We can introduce a diagonal smoothing matrix $\mathbf{S} = \text{diag}(\mathbf{s})$, where $\mathbf{s} \in \mathbb{R}^{C_i}$ and $s_j > 0$ represents the per-channel scaling factor for channel $j$:

$$\mathbf{Y} = \mathbf{X} \mathbf{W} = \left( \mathbf{X} \mathbf{S}^{-1} \right) \left( \mathbf{S} \mathbf{W} \right) = \hat{\mathbf{X}} \hat{\mathbf{W}}$$

In component notation, for each token $t$, input channel $j$, and output channel $k$:

$$\hat{X}_{t, j} = \frac{X_{t, j}}{s_j}, \quad \hat{W}_{j, k} = s_j \cdot W_{j, k}$$

Because multiplying activation channel $j$ by $s_j^{-1}$ is exactly offset by multiplying row $j$ of the weight matrix by $s_j$, the mathematical output $\mathbf{Y}$ remains identical.

Activation X (Spiky)       Weight W (Flat)
[ ...  98.0  ... ]    x    [ ...   0.02  ... ]
       |                          |
   Divide by s_j              Multiply by s_j (e.g., s_j = 10.0)
       v                          v
Smoothed X_hat             Adjusted W_hat
[ ...   9.8  ... ]    x    [ ...   0.20  ... ]
(Easy to Quantize)         (Easy to Quantize)

Derivation of the Migration Strength Formulation

To minimize total quantization error, the smoothing vector $\mathbf{s}$ must balance the dynamic range between activations and weights.

If we choose $s_j = \max(|X_j|)$, where $\max(|X_j|)$ is the maximum activation magnitude of channel $j$ across all tokens, the smoothed activation channel satisfies:

$$\max(|\hat{X}_j|) = \frac{\max(|X_j|)}{\max(|X_j|)} = 1.0$$

While this completely flattens the activation matrix, it transfers the entire outlier magnitude directly into row $j$ of the weight matrix:

$$\max(|\hat{W}_j|) = \max(|X_j|) \cdot \max(|W_j|)$$

This causes large quantization errors in the weight tensor. Conversely, setting $s_j = 1 / \max(|W_j|)$ leaves activations unaltered while smoothing weights.

To find the optimal operating point, SmoothQuant defines the per-channel smoothing factor via a migration strength hyperparameter $\alpha \in [0, 1]$:

$$s_j = \frac{\max(|X_j|)^\alpha}{\max(|W_j|)^{1-\alpha}}$$

Substituting $s_j$ into the smoothed activations and weights yields their scaled channel maxima:

$$\max(|\hat{X}_j|) = \frac{\max(|X_j|)}{s_j} = \max(|X_j|)^{1-\alpha} \cdot \max(|W_j|)^{1-\alpha} = \left( \max(|X_j|) \cdot \max(|W_j|) \right)^{1-\alpha}$$

$$\max(|\hat{W}_j|) = s_j \cdot \max(|W_j|) = \max(|X_j|)^\alpha \cdot \max(|W_j|)^\alpha = \left( \max(|X_j|) \cdot \max(|W_j|) \right)^\alpha$$

Properties of the Formulation

1. Symmetric Difficulty Split ($\alpha = 0.5$): When $\alpha = 0.5$, we have $\max(|\hat{X}_j|) = \max(|\hat{W}_j|) = \sqrt{\max(|X_j|) \max(|W_j|)}$. The quantization difficulty is split equally between activations and weights. Both tensors share identical channel maximum magnitudes, maximizing the effective bits for both operands.
2. Robustness Across Architectures: For standard architectures such as OPT, BLOOM, and Llama, $\alpha = 0.5$ preserves floating-point accuracy across all zero-shot benchmarks.
3. Asymmetric Tuning for Heavy Outliers ($\alpha \approx 0.75\text{--}0.85$): In architectures where activation outliers are particularly severe (such as GLM-130B, where outlier channels comprise up to 30% of certain projection inputs), setting $\alpha = 0.75$ migrates more dynamic range into the weights, where INT8 quantization error remains well tolerated.

Offline Layer Fusion: Eliminating Runtime Scaling Overhead

A critical design requirement for production serving is that smoothing transformations must not introduce runtime kernel latency or additional memory allocations. SmoothQuant achieves this by absorbing $\mathbf{S}^{-1}$ into preceding layers offline and quantizing the scaled weights $\hat{\mathbf{W}}$ ahead of time.

Standard Transformer Block:
[LayerNorm / RMSNorm] ---> [Linear Q, K, V Projections] ---> [Attention Core]

SmoothQuant Offline Fusion:
[Fold s^-1 into LayerNorm Weights/Bias] ---> [INT8 Linear GEMM (W_hat)]
(No runtime scaling kernel required)

1. Fusing into LayerNorm and RMSNorm

In modern Transformer blocks, linear layers (such as the Query, Key, Value projections in self-attention, and the Gate/Up projections in the Feed-Forward Network) are immediately preceded by normalization layers.

For standard LayerNorm:

$$\mathbf{X} = \text{LayerNorm}(\mathbf{Z}) = \frac{\mathbf{Z} - \mu}{\sigma} \odot \mathbf{\gamma} + \mathbf{\beta}$$

To compute the smoothed activation $\hat{\mathbf{X}} = \mathbf{X} \cdot \text{diag}(\mathbf{s})^{-1}$, we scale the learnable affine parameters offline:

$$\hat{\gamma}_j = \frac{\gamma_j}{s_j}, \quad \hat{\beta}_j = \frac{\beta_j}{s_j}$$

At inference time, executing LayerNorm with $\hat{\mathbf{\gamma}}$ and $\hat{\mathbf{\beta}}$ directly produces $\hat{\mathbf{X}}$ with zero additional floating-point operations.

For RMSNorm (common in Llama, Mistral, and modern open-weight models):

$$\mathbf{X} = \text{RMSNorm}(\mathbf{Z}) = \frac{\mathbf{Z}}{\text{RMS}(\mathbf{Z})} \odot \mathbf{\gamma} \implies \hat{\gamma}_j = \frac{\gamma_j}{s_j}$$

2. Fusing into Preceding Linear Layers

When an activation is produced by a preceding linear projection (such as the Self-Attention output projection feeding into the FFN layer norm, or the MLP down-projection), the smoothing factor $\mathbf{S}^{-1}$ is folded directly into the columns of the previous weight matrix and its bias:

$$\mathbf{X} = \mathbf{X}_{\text{in}} \mathbf{W}_{\text{prev}} + \mathbf{b}_{\text{prev}}$$

$$\hat{\mathbf{X}} = \mathbf{X} \cdot \text{diag}(\mathbf{s})^{-1} = \mathbf{X}_{\text{in}} \left( \mathbf{W}_{\text{prev}} \cdot \text{diag}(\mathbf{s})^{-1} \right) + \left( \mathbf{b}_{\text{prev}} \cdot \text{diag}(\mathbf{s})^{-1} \right)$$

$$\hat{\mathbf{W}}_{\text{prev}, :, j} = \frac{\mathbf{W}_{\text{prev}, :, j}}{s_j}, \quad \hat{b}_{\text{prev}, j} = \frac{b_{\text{prev}, j}}{s_j}$$

3. Handling Residual Connections

For residual additions where activations from two branches merge, a lightweight fused scaling vector is applied to the residual branch before addition, maintaining numerical consistency across residual streams.

SmoothQuant Precision Topologies: O1, O2, and O3

SmoothQuant defines three progressive efficiency levels based on the quantization granularity applied to weights and activations:

| Mode | Weight Quantization | Activation Quantization | Hardware Kernel Mapping | Target Use Case | | :--- | :--- | :--- | :--- | :--- | | O1 | Per-Channel (INT8) | Per-Token Dynamic (INT8) | Standard INT8 GEMM (cuBLAS / CUTLASS) | Maximum zero-shot accuracy preservation across all models | | O2 | Per-Tensor (INT8) | Per-Token Dynamic (INT8) | Standard INT8 GEMM | Intermediate efficiency with simplified weight scaling | | O3 | Per-Tensor (INT8) | Per-Tensor Static (INT8) | Pure fixed-point INT8 MMA (No dynamic scale reduction) | Maximum throughput and lowest latency on dedicated accelerators |

Detailed Execution Mechanics

O1: Per-Token Dynamic + Per-Channel Weight

In O1, weights are statically quantized per-channel offline: $\Delta_{W, k} = \max(|W_{:, k}|) / 127$. At runtime, after LayerNorm emits smoothed activations $\hat{\mathbf{X}}$, a dynamic reduction kernel computes the maximum per token $t$: $\Delta_{X, t} = \max(|\hat{X}_{t, :}|) / 127$. The INT8 GEMM kernel multiplies $\bar{\mathbf{X}}^{\text{INT8}}$ and $\bar{\mathbf{W}}^{\text{INT8}}$, and the output is scaled in FP16 epilogue code:

$$Y_{t, k} = \left( \sum_{j=1}^{C_i} \bar{X}_{t, j}^{\text{INT8}} \bar{W}_{j, k}^{\text{INT8}} \right) \cdot \Delta_{X, t} \cdot \Delta_{W, k}$$

O3: Per-Tensor Static + Per-Tensor Weight

In O3, both activation step sizes $\Delta_X$ and weight step sizes $\Delta_W$ are precomputed statically during calibration. This eliminates the dynamic per-token activation reduction kernel, allowing the GPU to execute pure fixed-point INT8 arithmetic with a single scalar output scaling factor $\Delta = \Delta_X \cdot \Delta_W$.

Transformer Block Mapping and Batched Matrix Multiplications

SmoothQuant quantizes all high-compute matrix operations within the Transformer architecture:

+-------------------------------------------------------------------+
|                        Transformer Block                          |
|                                                                   |
| Input Activation Z                                                |
|        |                                                          |
|        v                                                          |
| [ Smoothed RMSNorm (gamma_hat = gamma / s_attn) ]                 |
|        |                                                          |
|        +-------------------+-------------------+                  |
|        | (INT8 W8A8)       | (INT8 W8A8)       | (INT8 W8A8)      |
|        v                   v                   v                  |
| [ Q Projection ]    [ K Projection ]    [ V Projection ]          |
|        |                   |                   |                  |
|        +---------> [ INT8 BMM1 ] <-------------+                  |
|                          |                                        |
|                          v                                        |
|                   [ FP16 Softmax ]                                |
|                          |                                        |
|                          +---------------------> [ INT8 BMM2 ]    |
|                                                        |          |
|                                                        v          |
|                                               [ INT8 Out-Proj ]   |
|                                                        |          |
|                                                        v          |
|                                                 Residual Add      |
+-------------------------------------------------------------------+

1. Self-Attention Projections: The Q, K, and V linear layers share the input smoothing factor $\mathbf{s}_{\text{attn}}$ derived from the LayerNorm output.
2. Attention Batched Matrix Multiplications (BMMs):

• BMM1 ($Q \cdot K^T$): Query and Key tensors are quantized to INT8 before attention score computation.
• Softmax: Softmax is executed in FP16 to maintain probability distribution precision and prevent numerical underflow.
• BMM2 ($\text{Score} \cdot V$): Attention probabilities are dynamically requantized to INT8 and multiplied by the INT8 Value tensor.

1. Feed-Forward Network (FFN):

• For SwiGLU / Gated architectures (Llama, Mistral), the Gate and Up projections share the smoothing scale $\mathbf{s}_{\text{ffn}}$ from the post-attention normalization layer.
• The Down projection absorbs the smoothing scale $\mathbf{s}_{\text{down}}$ computed across intermediate activation states.

Empirical Benchmarks and Systems Evaluation

SmoothQuant was extensively evaluated across multiple open model families up to 530 billion parameters, benchmarked against FP16 baselines and prior quantization methods.

Zero-Shot Accuracy on OPT-175B

On the 175-billion parameter OPT model, standard naive W8A8 quantization and existing post-training methods fail catastrophically due to outlier corruption. SmoothQuant maintains baseline FP16 accuracy across zero-shot natural language understanding benchmarks and WikiText perplexity:

| Quantization Method | Precision | Average Zero-Shot Accuracy | WikiText Perplexity ($\downarrow$) | | :--- | :--- | :--- | :--- | | FP16 Baseline | FP16 | 66.9% | 10.99 | | Naive W8A8 | INT8 | 35.5% | 93,080.00 | | ZeroQuant | INT8 | 35.8% | 84,648.00 | | LLM.int8() | Mixed INT8/FP16 | 66.7% | 11.10 | | Outlier Suppression | INT8 | 36.0% | 96,151.00 | | SmoothQuant-O1 | INT8 (W8A8) | 66.5% | 11.11 | | SmoothQuant-O2 | INT8 (W8A8) | 66.4% | 11.14 | | SmoothQuant-O3 | INT8 (W8A8) | 66.8% | 11.17 |

Accuracy Across Frontier LLM Architectures

SmoothQuant generalizes across diverse architectures without per-model manual heuristic tuning:

| Model | Parameter Count | FP16 Accuracy | SmoothQuant-O1 | SmoothQuant-O2 | SmoothQuant-O3 | | :--- | :--- | :--- | :--- | :--- | :--- | | OPT-175B | 175B | 71.6% | 71.2% | 71.1% | 71.1% | | BLOOM-176B | 176B | 68.2% | 68.3% | 68.4% | 67.4% | | GLM-130B | 130B | 73.8% | 73.7% | 72.5% | 72.8% |

(Evaluated on standard benchmark suites: WinoGrande, HellaSwag, PIQA, LAMBADA, MMLU, and MNLI).

Hardware Serving Throughput and Memory Scaling

When integrated into high-performance serving runtimes such as FasterTransformer and TensorRT-LLM:

1. Hardware Memory Halved: SmoothQuant reduces the memory footprint of weights and activations by up to 50%, enabling a 175B model to run on 4x 80GB A100 GPUs instead of the 8x 80GB GPUs required for FP16.
2. Single-Node 530B Serving: SmoothQuant enabled serving the MT-NLG 530B parameter model within a single 8-GPU node (8x 80GB A100s), where FP16 required two full DGX nodes (16 GPUs).
3. Execution Speedup: On NVIDIA Ampere (A100) and Hopper (H100) architectures, SmoothQuant delivers up to 1.56x end-to-end inference speedup over FP16 baselines, overcoming the 20-30% latency penalty introduced by mixed-precision decomposition frameworks.

Calibration and Implementation Guidelines

Applying SmoothQuant in practice involves a simple two-phase workflow:

Step 1: Activation Calibration

To determine the per-channel activation scales $\max(|X_j|)$, run forward passes over a modest calibration dataset (typically 512 random sequences drawn from general corpora such as The Pile or C4).

import torch

@torch.no_grad()
def get_act_scales(model, dataloader, num_samples=512):
    model.eval()
    act_scales = {}

    def stat_tensor(name, tensor):
        hidden_dim = tensor.shape[-1]
        tensor = tensor.view(-1, hidden_dim).abs().detach()
        comming_max = torch.max(tensor, dim=0)[0].float().cpu()
        if name in act_scales:
            act_scales[name] = torch.max(act_scales[name], comming_max)
        else:
            act_scales[name] = comming_max

    def hook_fn(name):
        return lambda module, inp, out: stat_tensor(name, inp[0])

    hooks = []
    for name, module in model.named_modules():
        if isinstance(module, torch.nn.Linear):
            hooks.append(module.register_forward_hook(hook_fn(name)))

    for i, batch in enumerate(dataloader):
        if i >= num_samples:
            break
        model(batch.to(model.device))

    for h in hooks:
        h.remove()
    return act_scales

Step 2: Smoothing and Layer Transformation

Given calibration scales act_scales and model weights:

@torch.no_grad()
def smooth_ln_fcs(ln, fcs, act_scales, alpha=0.5):
    if not isinstance(fcs, list):
        fcs = [fcs]

    # Compute weight maximums across input channels
    weight_scales = torch.cat([fc.weight.abs().max(dim=0, keepdim=True)[0] for fc in fcs], dim=0)
    weight_scales = weight_scales.max(dim=0)[0].clamp(min=1e-5)

    act_scales = act_scales.to(device=fcs[0].weight.device).clamp(min=1e-5)

    # Compute smoothing scale vector s
    scales = (act_scales.pow(alpha) / weight_scales.pow(1.0 - alpha)).clamp(min=1e-5)

    # Fold s^-1 into LayerNorm parameters
    ln.weight.div_(scales)
    if hasattr(ln, "bias") and ln.bias is not None:
        ln.bias.div_(scales)

    # Scale weight matrices by s
    for fc in fcs:
        fc.weight.mul_(scales.view(1, -1))

Following transformation, all linear layers can be exported directly to standard INT8 format for deployment in TensorRT-LLM, vLLM, or FasterTransformer engines.

Sources

• Xiao, G., Lin, J., Seznec, M., Wu, H., Demouth, J., & Han, S. (2023). SmoothQuant: Accurate and Efficient Post-Training Quantization for Large Language Models. International Conference on Machine Learning (ICML 2023).
• Dettmers, T., Lewis, M., Belkada, Y., & Zettlemoyer, L. (2022). LLM.int8(): 8-bit Matrix Multiplication for Transformers at Scale. Advances in Neural Information Processing Systems (NeurIPS 2022).
• Yao, Z., Aminabadi, R. Y., Zhang, M., Wu, X., Li, C., & He, Y. (2022). ZeroQuant: Efficient and Affordable Post-Training Quantization for Large-Scale Transformers. Advances in Neural Information Processing Systems (NeurIPS 2022).
• Frantar, E., Saleh, S., Zhang, D., & Alistarh, D. (2022). GPTQ: Accurate Post-Training Quantization for Generative Pre-Trained Transformers. International Conference on Learning Representations (ICLR 2023).
• Lin, J., Tang, J., Tang, H., Yang, S., Chen, W.-M., Wang, W.-C., & Han, S. (2023). AWQ: Activation-aware Weight Quantization for LLM Compression and Acceleration. Machine Learning and Systems (MLSys 2024).
• Wei, X., Zhang, Y., Zhang, X., Gong, R., Zhang, S., Zhang, Q., Yu, F., & Liu, X. (2022). Outlier Suppression: Pushing the Limit of Low-Bit Transformer Language Models. Advances in Neural Information Processing Systems (NeurIPS 2022).
• Bondarenko, Y., Nagel, M., & Blankevoort, T. (2021). Understanding and Overcoming the Challenges of Efficient Transformer Quantization. Empirical Methods in Natural Language Processing (EMNLP 2021).