BYOL

BYOL (Bootstrap Your Own Latent; Grill et al., NeurIPS 2020) is a self-supervised method for learning image representations without labels and without negative pairs. Two networks see two augmented views of the same image. An online network, followed by a small predictor, is trained to predict what a target network outputs for the other view. The target is not trained by gradient descent: its weights are an exponential moving average (EMA) of the online weights, and the target branch sits behind a stop-gradient. The backbone in the experiments discussed here is a ResNet-50 convolutional network, and later a Vision Transformer.

What makes BYOL interesting is what it leaves out. Contrastive methods such as SimCLR and MoCo push the embeddings of different images apart with an InfoNCE-type loss; without that push, a network that maps every image to the same vector would minimize a pure agreement loss. Tian, Chen and Ganguli (2021) call the capacity of BYOL and SimSiam to learn strong representations from positive pairs alone “remarkable”, and set out to explain why they do not collapse. This article covers that question in depth: the architecture, the evidence on which component prevents collapse, SimSiam (Chen and He, 2021), which showed that the moving average can be dropped, Tian et al.’s linear theory, and a deterministic integration of that theory on image patches.

Sources. The BYOL paper itself was not read for this article: the NeurIPS proceedings site did not respond when this batch was planned. BYOL is described here only through papers that discuss and reproduce it (SimSiam; Tian, Chen and Ganguli; MoCo v3) and through published articles on this site. Every statement about BYOL’s own paper is attributed to the source that reports it.

Problem Addressed

Two-view self-supervised learning trains an encoder so that two augmentations of one image get similar representations. The loss has a trivial minimum: a constant output. Methods published before BYOL ruled it out in different ways, which Self-Supervised Learning surveys. Contrastive learning uses negatives, from a large batch (SimCLR) or from a queue maintained by a momentum encoder (MoCo). Clustering methods such as SwAV assign views to balanced clusters. Chen and He (2021, §1–2) summarize the alternatives this way and note that contrastive methods benefit from many negatives, kept in a memory bank, a queue or a large batch, and that clustering methods likewise need a memory bank, large batches or a queue.

BYOL, in Chen and He’s description, relies only on positive pairs and does not collapse provided a momentum encoder is used. Tian et al. (2021, §1) list what this removes: the large batches (4,096 in SimCLR) or memory queue (MoCo) that supply negatives. Doing away with negatives also removes the question of which images count as different, a question that Contrastive Learning and InfoNCE Loss discuss for the contrastive family.

Architecture

BYOL is a Siamese network with what Chen and He call “indirect” weight sharing: the two branches have the same architecture, and the target’s weights follow the online weights through the moving average rather than being shared directly (Chen & He, 2021, footnote 2). Chen and He’s Figure 3 draws the two branches as encoder, predictor and similarity on the online side, and a momentum encoder updated by a moving average on the other, with gradients flowing only through the online side.

Encoder and Projector

Each branch has a backbone followed by a projection MLP. In Tian et al.’s ImageNet reimplementation, which they say follows BYOL strictly in architecture, augmentations and linear evaluation, a ResNet-50 backbone produces the features used by the linear probe and is followed by a projector and a predictor (Tian et al., 2021, §4).

Predictor

The predictor is an extra MLP on the online branch only. Tian et al. describe BYOL’s default as a two-layer predictor with batch normalization and a ReLU, a 4,096-dimensional hidden layer and 256-dimensional input and output. It maps the online projection of view 1 to a prediction of the target projection of view 2. The target branch has no predictor, so the two branches are deliberately asymmetric.

Target Network

The target has the same backbone and projector. Its weights $\xi$ follow the online weights $\theta$ by an exponential moving average after every step,

$$\xi \leftarrow \gamma\, \xi + (1 - \gamma)\, \theta ,$$

with $\gamma = 0.996$ in BYOL according to Tian et al. (§2). The notation here is ours. In their reproduction recipes, Chen and He used momentum coefficients of 0.99, 0.99 and 0.996 for 100, 200–400 and 800 epochs, following BYOL’s 100-, 300- and 1,000-epoch recipes (Chen & He, 2021, footnote 7). Chen and He describe MoCo as turning one branch into a momentum encoder to keep its queue of negatives consistent; BYOL keeps the momentum encoder but has no queue and no negatives.

Training

Objective

Let $z_2$ be the target’s projection of view 2 and $p_1$ the online prediction from view 1. Chen and He note that BYOL’s loss is the mean squared error of $\ell_2$-normalized vectors, which equals the negative cosine similarity up to a scale of 2 (Chen & He, 2021, §3). With $\bar{u} = u/\lVert u \rVert_2$,

$$\mathcal{L}(p_1, z_2) = \big\lVert \bar{p}_1 - \mathrm{sg}(\bar{z}_2) \big\rVert_2^2 = 2 - 2\,\frac{p_1^\top z_2}{\lVert p_1\rVert_2 \lVert z_2\rVert_2},$$

where $\mathrm{sg}$ marks the stop-gradient: the target branch gets no gradient. The loss is symmetrized by also predicting view 1’s target from view 2. Chen and He define their symmetrized loss “following” BYOL, and Tian et al.’s ImageNet setting that “follows BYOL more closely” uses a symmetrized loss, a batch size of 4,096, the LARS optimizer and 300 epochs (Tian et al., 2021, §4). Tian et al. list color jittering, blur and solarization among the augmentation strategies their reimplementation took from BYOL, and DINO records that DINO adopted BYOL’s augmentations.

Which Component Prevents Collapse

Three components separate BYOL from a plain agreement loss: the predictor, the stop-gradient and the slow target. The readable evidence on their roles comes from SimSiam, from Tian et al.’s ablations, and from what both papers report of BYOL’s own ablations.

BYOL’s own account, as others report it. Chen and He (2021, §2) write that BYOL hypothesized the momentum encoder to be important for avoiding collapse, and that BYOL reported failure without it (0.3% accuracy, BYOL’s Table 5). Chen and He add that BYOL’s arXiv v3 reports 66.9% with 300-epoch pre-training when the momentum encoder is removed and the predictor’s learning rate is raised tenfold (footnote 3). Tian et al. likewise report that BYOL needed the EMA at first but later confirmed it is not necessary if the predictor has a higher relative learning rate (§1). They cite BYOL’s Table 15 for unstable results without weight decay (§1) and its Table 22 for the predictor learning rate (§3.2).

SimSiam. Chen and He removed the moving average entirely. SimSiam shares the encoder weights directly between the branches, keeps a prediction MLP $h$ on one side and a stop-gradient on the other, and minimizes the symmetrized negative cosine similarity:

$$\mathcal{L} = \tfrac12 D\big(p_1, \mathrm{sg}(z_2)\big) + \tfrac12 D\big(p_2, \mathrm{sg}(z_1)\big), \qquad D(p, z) = -\frac{p}{\lVert p\rVert_2} \cdot \frac{z}{\lVert z\rVert_2}.$$

In their own words, SimSiam can be thought of as “BYOL without the momentum encoder”. Their ablations on ImageNet (ResNet-50, 100-epoch pre-training, linear evaluation) separate the three components:

SimSiam variant (Chen & He, 2021)	Linear acc. (%)
Baseline (predictor, stop-gradient)	67.7 ± 0.1 (5 trials)
Without stop-gradient	0.1
Without predictor	0.1
Predictor fixed at random initialization	1.5
Predictor learning rate not decayed	68.1

Without the stop-gradient, the optimizer quickly reached the loss’s minimum of −1 and the per-channel standard deviation of the normalized output fell to zero, the signature of a constant output; with it, the standard deviation stayed near $1/\sqrt{d}$, the value for outputs spread over the unit sphere (Figure 2). Chen and He stress that the collapse is not prevented by architecture alone: the predictor, batch normalization and $\ell_2$ normalization were all kept when the stop-gradient was removed. Without the predictor, they observe, the symmetrized loss with stop-gradient has a gradient in the same direction as the loss without stop-gradient, scaled by 1/2, so removing the predictor amounts to removing the stop-gradient. A random fixed predictor failed differently: the loss stayed high and training did not converge, which they distinguish from collapse. Their conclusion is that the moving average may help accuracy with the right momentum coefficient but is not what prevents collapse; the stop-gradient is critical, and its role was obscured in BYOL because a momentum encoder always comes with a stop-gradient (§2, §4.1, §4.2).

Chen and He’s explanation is a hypothesis, not a proof: SimSiam behaves like an alternating, expectation–maximization-like optimization over the network weights and a second set of variables, one representation per image (§5). The stop-gradient follows from holding one set fixed while the other is updated, and the predictor, they suggest, approximates the expectation over augmentations that one sampled view leaves out. Supporting this, replacing the predictor by a per-image moving average of past outputs gave 55.0% without a predictor, where the plain predictor-free version failed. They state that the hypothesis describes what is optimized but does not explain why collapse is prevented, which remains an empirical observation (§5.3).

Self-Supervised Learning keeps a short summary of this argument among its other collapse-avoidance strategies; DINO shows a third strategy, centering and sharpening a momentum teacher, and reports that adding a BYOL-style predictor to DINO changed little.

Tian, Chen and Ganguli’s Linear Theory

Tian, Chen and Ganguli (2021) analyse a bias-free linear version of BYOL in which the online network, the predictor and the target are matrices $W$, $W_p$ and $W_a$. With $f_1 = W x_1$ and $f_{2a} = W_a x_2$ for two augmented views of an input $x$, the online network and the predictor minimize

$$J(W, W_p) = \tfrac12\, \mathbb{E}_{x_1, x_2}\Big[\big\lVert W_p f_1 - \mathrm{StopGrad}(f_{2a}) \big\rVert_2^2\Big],$$

while $W_a$ follows $W$ by an EMA. In the limit of large batches and small learning rates, training becomes a gradient flow (their Lemma 1, Eqs. 2–4):

$$\begin{aligned} \dot W_p &= \alpha_p\big(-W_p W (X + X') + W_a X\big) W^\top - \eta W_p, \\ \dot W &= W_p^\top\big(-W_p W (X + X') + W_a X\big) - \eta W, \\ \dot W_a &= \beta\,(-W_a + W). \end{aligned}$$

Here $X = \mathbb{E}[\bar{x}\bar{x}^\top]$ is the correlation of the average augmented view $\bar x$ of each input, $X'$ is the covariance of the augmented views around that average, averaged over inputs, $\alpha_p$ and $\beta$ are the predictor’s and the EMA’s rates relative to the online network, and $\eta$ is weight decay. SimSiam is the case $W_a = W$, with the stop-gradient still in place. The time $t$ is gradient-flow time: one step of learning rate $\alpha$ advances it by $\alpha$, so the per-step EMA rate is $\alpha\beta = 1 - \gamma$.

From this system they derive the main results:

• Balancing (Theorem 1). Whatever the target does, $W W^\top = \alpha_p^{-1} W_p^\top W_p + e^{-2\eta t} C$, with $C$ fixed by the initialization. What the predictor learns, the online network learns too, so the match between branches cannot be carried by the predictor alone; weight decay erases the memory of the initialization $C$.
• The stop-gradient is essential (Theorem 2). With $W_a = W$ and no stop-gradient, $\mathrm{vec}(W)$ obeys $\frac{d}{dt}\mathrm{vec}(W) = -H(t)\,\mathrm{vec}(W)$ with $H(t) = X' \otimes (W_p^\top W_p + I) + X \otimes \tilde W_p^\top \tilde W_p + \eta I$, $\tilde W_p = W_p - I$. $H$ is positive semi-definite, and if its smallest eigenvalue stays bounded away from zero, $W(t) \to 0$: collapse. Their supplement shows that removing the predictor ($W_p = I$) with the stop-gradient kept gives $\dot W = -(X' + \eta I) W$, so $W \to 0$ as well.
• Eigenspace alignment (Theorem 3). Under three simplifying assumptions (a target proportional to the online network, $W_a = \tau W$; isotropic data, $X = I$ and $X' = \sigma^2 I$; a symmetric predictor), the commutator of $W_p$ with $F = W X W^\top$, the correlation of the predictor’s inputs, decays to zero whenever the smallest eigenvalue of a matrix $K(t)$ stays above a positive constant. The predictor’s eigenvectors align with those of its input correlation, which they observed in BYOL training of a ResNet-18 on STL-10 (Figure 2), where $W_p$ also became nearly symmetric.

Once aligned, each eigenmode $j$ has three scalar variables, the predictor eigenvalue $p_j$, the eigenvalue $s_j$ of $W W^\top$, and $\tau$ (Eqs. 11–13):

$$\dot p_j = \alpha_p s_j\big[\tau - (1+\sigma^2)p_j\big] - \eta p_j, \qquad \dot s_j = 2 p_j s_j \big[\tau - (1+\sigma^2)p_j\big] - 2\eta s_j, \qquad s_j \dot\tau = \beta(1-\tau)s_j - \tau \dot s_j / 2 .$$

They draw these conclusions from the per-mode system (§3.2):

• An exact integral. $s_j(t) = \alpha_p^{-1} p_j^2(t) + e^{-2\eta t} c_j$ (Eq. 14).
• No weight decay ($\eta = 0$). There is a collapsed branch $s_j = 0$ and a non-collapsed branch $p_j = \tau/(1+\sigma^2)$, which models invariance: stronger augmentation means a smaller predictor eigenvalue. Since $c_j$ is conserved, $s_j(\infty) = s_j(0) + \alpha_p^{-1}(p_j^2(\infty) - p_j^2(0))$. A fast predictor (large $\alpha_p$) shrinks the basin of collapse (Obs. 1), but too fast a predictor leaves $s_j$ where it started (Obs. 2), consistent with SimSiam’s finding that updating the predictor too often hurts.
• With weight decay ($\eta > 0$). The trajectory is drawn to the invariant parabola $s_j = \alpha_p^{-1}p_j^2$ (Obs. 3), on which $\dot p_j = p_j^2[\tau - (1+\sigma^2)p_j] - \eta p_j$. For $0 < \eta \le \tau^2/(4(1+\sigma^2))$ there are three fixed points, $p^*_{j0} = 0$ and

$$p^*_{j\pm} = \frac{\tau \pm \sqrt{\tau^2 - 4\eta(1+\sigma^2)}}{2(1+\sigma^2)} ;$$

$0$ and $p^*_{j+}$ are stable and $p^*_{j-}$ is unstable, so a mode collapses if it starts below $p^*_{j-}$. Larger $\eta$ widens that basin (Obs. 4), and above the bound only the collapsed point remains (Obs. 5).

• The EMA. Tian et al. argue that the EMA acts as an automatic curriculum (Obs. 8). Early on, $\tau$ is small and sets a low target $\tau/(1+\sigma^2)$ for $p_j$, which keeps $K(t)$ positive definite and the alignment stable. A very slow EMA slows training (Obs. 9), and a small $\tau$ raises $p^*_{j-}$, trapping more modes (Obs. 10).

A sign in the paper. The main text after Eq. 14 defines $c_j = \alpha_p^{-1}p_j^2(0) - s_j(0)$. Eq. 14 at $t = 0$, and the supplement’s derivation (Eqs. 82–83), give the opposite sign, $c_j = s_j(0) - \alpha_p^{-1}p_j^2(0)$. This article uses the supplement’s sign, the one consistent with Eq. 14.

The theory matched several of their STL-10 ablations (ResNet-18, 100 epochs, top-1 linear accuracy, Table 3):

• With an EMA, a symmetric one-layer linear predictor gave 75.09 ± 0.48 against 74.51 ± 0.47 for an unconstrained one.
• Without an EMA, the symmetric predictor crashed (36.62 ± 1.85), while the unconstrained one reached 72.85 ± 0.16.
• On ImageNet, a symmetric linear predictor in SimSiam scored 0.1% top-1, against 68.1% for an asymmetric one (§3).
• Raising the predictor’s relative learning rate above 1 rescued the symmetric predictor without an EMA (Figure 6), as the theory predicts.
• Removing weight decay lowered BYOL’s STL-10 accuracy to between 70.62 and 71.37, depending on EMA and predictor bias (Table 2).

DirectPred. The theory suggests skipping the predictor’s gradient training altogether. DirectPred estimates the correlation matrix $\hat F$ of the predictor’s inputs with a moving average, takes its eigendecomposition $\hat F = \hat U \operatorname{diag}(s_j) \hat U^\top$, and sets

$$W_p = \hat U \operatorname{diag}(p_j) \hat U^\top, \qquad p_j = \sqrt{s_j} + \epsilon \max_j s_j ,$$

the eigenbasis from Theorem 3 and the square root from the invariant parabola (Eq. 18). In their ImageNet reimplementation of BYOL (ResNet-50, Table 9), the results were as follows; the 60-epoch runs used an asymmetric loss, SGD and batch size 256, the 300-epoch runs the symmetrized, batch-4,096 LARS setting above (§4):

Predictor	60 epochs, top-1	300 epochs, top-1	300 epochs, top-5
BYOL’s default two-layer predictor	64.7	72.5	90.8
Learned linear predictor	59.4	69.9	89.6
DirectPred	64.4	72.4	91.0

Tian et al. state that the 72.5% two-layer result reproduces the accuracy reported in BYOL.

Collapse and Survival in the Linear Model

The script below integrates Tian et al.’s dynamics without their isotropy assumption. Image patches supply the input statistics, and the code checks the integration against every closed form the theory gives. It uses NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2 with one BLAS thread. It runs in about 10 seconds and gives identical output on two runs. Nothing is trained on samples. Each run integrates the population gradient flow of Lemma 1 deterministically with scipy.integrate.solve_ivp.

Data. The script cuts non-overlapping $8 \times 8$ patches from three bundled scikit-image photographs: camera, and astronaut and coffee converted to grey (11,942 patches). It subtracts each patch’s own mean brightness and sets $X$ to the patches’ $64 \times 64$ second-moment matrix, scaled so its largest eigenvalue is $\lambda_1 = 1$. The removed brightness direction leaves one eigenvalue at zero.

Augmentation and parameters. Augmentation is modelled as additive isotropic noise, $x' = x + \varepsilon$ with $\varepsilon \sim \mathcal{N}(0, \sigma^2 I)$. This gives Tian et al.’s $X$ as above and $X' = \sigma^2 I$, with $\sigma^2 = 0.03$. The other parameters are $\alpha_p = 1$ and $\beta = 0.4$, the EMA rate in the caption of Tian et al.’s Figure 5. Weight decay is $\eta \in \{0, 10^{-4}, 10^{-3}, 10^{-2}\}$.

Per-mode equations. If $W$, $W_p$ and $W_a$ are diagonal in $X$’s eigenbasis, every product in Lemma 1 stays diagonal. Each mode $k$, with eigenvalue $\lambda_k$ of $X$, then evolves by the scalar version of Lemma 1, with $w$, $p$ and $a$ the diagonal entries of $W$, $W_p$ and $W_a$. Tian et al.’s supplement gives this reduction without weight decay for non-isotropic data (Appendix D.4, Eqs. 102–104), with the eigenvalues $\lambda_s = \lambda_k + \sigma^2$ and $\lambda_d = \lambda_k$ of their single-view and two-view correlations $X + X'$ and $X$; the weight-decay terms are our addition:

$$G = -p\,w\,(\lambda_k + \sigma^2) + a \lambda_k, \qquad \dot p = \alpha_p G w - \eta p, \qquad \dot w = p\, G - \eta w, \qquad \dot a = \beta (w - a).$$

With $\lambda_k = 1$, $s = w^2$ and $\tau = a/w$, these are exactly Tian et al.’s Eqs. 11–13; the script checks this. For other $\lambda_k$, the non-collapsed branch at $\eta = 0$ is $p_\infty = \lambda_k/(\lambda_k + \sigma^2)$ when $\tau = 1$, the supplement’s $w_p = \lambda_d \lambda_s^{-1}$: the predictor eigenvalue is the mode’s share of signal in signal plus augmentation noise. With weight decay, the fixed points on the invariant parabola become

$$p^*_{\pm}(\lambda_k) = \frac{\lambda_k \pm \sqrt{\lambda_k^2 - 4\eta(\lambda_k + \sigma^2)}}{2(\lambda_k + \sigma^2)} .$$

A mode can therefore survive only if $\lambda_k^2 \ge 4\eta(\lambda_k + \sigma^2)$. Weight decay keeps the strong directions of the patch distribution and removes the weak ones. The fixed points $p^*_\pm(\lambda_k)$ and this survival condition are our derivation; Tian et al. derive the weight-decay fixed points only in the isotropic case $\lambda_k = 1$.

Ablations (per mode). Without a predictor ($p \equiv 1$) and with $a = w$, the flow is $\dot w = -(\sigma^2 + \eta) w$, Tian et al.’s Remark, with closed form $w_0 e^{-(\sigma^2+\eta)t}$. Without a predictor but with the EMA, $(w, a)$ follow a linear system with matrix $\begin{pmatrix} -(\lambda + \sigma^2 + \eta) & \lambda \\ \beta & -\beta \end{pmatrix}$. Its determinant $\beta(\sigma^2 + \eta)$ is positive and its trace is negative, so it decays as well (our derivation); the script compares it with the matrix exponential. Without the stop-gradient (SimSiam with gradients through both branches), each mode follows Theorem 2’s form, $\dot w = -w[(p^2+1)\sigma^2 + (p-1)^2\lambda] - \eta w$. Since $(p^2 + 1)\sigma^2 \ge \sigma^2$, this gives $|w(t)| \le |w_0| e^{-(\sigma^2+\eta)t}$.

Conventions.

• $s = w^2$ is the eigenvalue of $W W^\top$ for the mode, as in Tian et al.
• A mode counts as surviving if $s(T) > 10^{-6}$.
• Times are gradient-flow times, not epochs, with $T = 2 \times 10^5$ for the per-mode runs.
• Two starting values of the target are printed, because the sources leave the choice open. “Target from 0” is $W_a(0) = 0$, that is $\tau(0) = 0$, the start Tian et al. use for Eq. 15 and in their simulations. “Target = online” starts the target as a copy of the online network.

Parts of the script.

• A checks the per-mode equations against Eqs. 11–13, and the full matrix flow against the per-mode flow from an aligned start.
• B runs $\eta = 0$ from $w_0 = 0.3$, $p_0 = 0.1$, a start off the parabola with $c = 0.08$. It runs SimSiam, BYOL with both target starts, and the three ablations.
• C sweeps $\eta$ from $w_0 = p_0 = 0.1$, on the parabola. It compares survivors with the fixed-point count (modes with $\lambda_k^2 \ge 4\eta(\lambda_k+\sigma^2)$) and with SimSiam’s exact prediction. On the parabola with $\tau \equiv 1$ the flow is one-dimensional, so SimSiam keeps exactly the modes whose fixed points exist and whose start $p_0$ lies above $p^*_-$.
• D integrates the full $64 \times 64$ BYOL flow from a random, unaligned start, with $\eta = 10^{-3}$, to $t = 2 \times 10^4$. It reports four quantities:
  • the predictor’s asymmetry $\lVert W_p - W_p^\top\rVert_F / \lVert W_p\rVert_F$;
  • the normalized commutator $\lVert [F, W_p]\rVert_F / (\lVert F\rVert_F \lVert W_p\rVert_F)$;
  • the error of Theorem 1’s invariant;
  • at the end, Tian et al.’s Figure 2 alignment measure $u^\top W_p u / \lVert W_p u\rVert$ over the surviving eigenvectors $u$ of $F$. It also reports how far the surviving input directions are from eigenvectors of $X$, and how far their eigenvalues are from the per-mode prediction $p^*_+(\lambda)^2/\alpha_p$.

import os
for v in ("OPENBLAS_NUM_THREADS", "OMP_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v] = "1"                       # one BLAS thread: reproducible sums
import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import expm
from skimage import data, color

# --- data: 8x8 patches of three bundled images, each patch's mean removed ---
ims = [data.camera() / 255.0, color.rgb2gray(data.astronaut()), color.rgb2gray(data.coffee())]
P = 8
pt = np.array([im[i:i+P, j:j+P].ravel() for im in ims
               for i in range(0, im.shape[0] - P + 1, P) for j in range(0, im.shape[1] - P + 1, P)])
pt = pt - pt.mean(1, keepdims=True)
X = pt.T @ pt / len(pt)
lam, U = np.linalg.eigh(X)
lam, U = lam[::-1].copy(), U[:, ::-1].copy()
X, lam = X / lam[0], lam / lam[0]
n = len(lam)
S2, AP, BETA = 0.03, 1.0, 0.4                 # augmentation variance, predictor lr ratio, EMA rate
T = 2e5
print(f"patches {len(pt)}, dim {n}; lambda_1 {lam[0]:.6f}, lambda_2 {lam[1]:.6f}, lambda_63 {lam[-2]:.3e}, "
      f"|lambda_64| {abs(lam[-1]):.1e}; X symmetric {np.allclose(X, X.T)}; sorted {np.all(np.diff(lam) <= 0)}")

# --- per-mode dynamics: Lemma 1 restricted to matrices diagonal in X's eigenbasis ---
def byol(t, y, eta):          # y = [w, p, a] for all modes
    w, p, a = y.reshape(3, n)
    G = -p * w * (lam + S2) + a * lam
    return np.concatenate([p * G - eta * w, AP * G * w - eta * p, BETA * (w - a)])

def simsiam(t, y, eta):       # y = [w, p], target = online (a = w), stop-gradient kept
    w, p = y.reshape(2, n)
    G = -p * w * (lam + S2) + w * lam
    return np.concatenate([p * G - eta * w, AP * G * w - eta * p])

def nopred_simsiam(t, y, eta):    # Wp = I, stop-gradient
    return -(S2 + eta) * y

def nopred_byol(t, y, eta):       # Wp = I, EMA target
    w, a = y.reshape(2, n)
    return np.concatenate([-w * (lam + S2) + a * lam - eta * w, BETA * (w - a)])

def nosg_simsiam(t, y, eta):      # a = w and the gradient flows through both branches
    w, p = y.reshape(2, n)
    dw = -w * ((p**2 + 1) * S2 + (p - 1)**2 * lam) - eta * w
    dp = AP * w**2 * (lam - p * (lam + S2)) - eta * p
    return np.concatenate([dw, dp])

TE = np.concatenate([[0.0, 50.0, 100.0], np.logspace(0, np.log10(T), 300)])
TE = np.unique(np.clip(TE, 0, T))
def run(f, y0, eta):
    s = solve_ivp(f, (0, T), y0, args=(eta,), method="LSODA", rtol=1e-10, atol=1e-14, t_eval=TE)
    assert s.success, s.message
    return s.t, s.y

def pstar(eta):
    disc = lam**2 - 4 * eta * (lam + S2)
    ok = disc >= 0
    r = np.sqrt(np.clip(disc, 0, None))
    return ok, np.where(ok, (lam - r) / (2 * (lam + S2)), np.inf), np.where(ok, (lam + r) / (2 * (lam + S2)), np.nan)

# --- A. exactness checks ---
print("\nA. exactness checks")
rng = np.random.default_rng(0)
w, p, a = rng.uniform(0.05, 1, 1000), rng.uniform(-1, 1, 1000), rng.uniform(-1, 1, 1000)
eta = 1e-3
G = -p * w * (1 + S2) + a * 1.0
dw, dp, da = p * G - eta * w, AP * G * w - eta * p, BETA * (w - a)
s_, tau = w**2, a / w
ds, dtau = 2 * w * dw, da / w - a * dw / w**2
e11 = AP * s_ * (tau - (1 + S2) * p) - eta * p
e12 = 2 * p * s_ * (tau - (1 + S2) * p) - 2 * eta * s_
e13 = BETA * (1 - tau) * s_ - tau * ds / 2
print(f"A1 per-mode RHS at lambda=1 vs Tian et al. Eqs. 11-13, max abs diff: "
      f"{max(np.abs(dp - e11).max(), np.abs(ds - e12).max(), np.abs(s_ * dtau - e13).max()):.2e}")

def matrix_rhs(t, y, eta):
    W, Wp, Wa = y.reshape(3, n, n)
    G = -Wp @ W @ (X + S2 * np.eye(n)) + Wa @ X
    return np.concatenate([(Wp.T @ G - eta * W).ravel(), (AP * G @ W.T - eta * Wp).ravel(), (BETA * (W - Wa)).ravel()])

def run_matrix(y0, eta, tmax, te):
    s = solve_ivp(matrix_rhs, (0, tmax), y0, args=(eta,), method="DOP853", rtol=1e-10, atol=1e-13, t_eval=te)
    assert s.success, s.message
    return s.t, s.y

ta = [0, 50, 100, 200]
W0, Wp0 = U @ np.diag(np.full(n, 0.3)) @ U.T, U @ np.diag(np.full(n, 0.1)) @ U.T
_, ym = run_matrix(np.concatenate([W0.ravel(), Wp0.ravel(), np.zeros(n * n)]), 1e-3, 200, ta)
sol = solve_ivp(byol, (0, 200), np.concatenate([np.full(n, 0.3), np.full(n, 0.1), np.zeros(n)]), args=(1e-3,),
                method="LSODA", rtol=1e-10, atol=1e-14, t_eval=ta)
err = 0.0
for k in range(len(ta)):
    Ms = [U.T @ M.reshape(n, n) @ U for M in ym[:, k].reshape(3, n * n)]
    for M, v in zip(Ms, sol.y[:, k].reshape(3, n)):
        err = max(err, np.abs(M - np.diag(v)).max())
print(f"A2 matrix ODE from an aligned start vs per-mode ODE (eta=1e-3, t<=200), max abs diff: {err:.2e}")

# --- B. no weight decay, start off the parabola: w0 = 0.3, p0 = 0.1, c = w0^2 - p0^2/ap ---
print("\nB. eta = 0, w0 = 0.3, p0 = 0.1, c = 0.08, T = 2e5")
w0, p0 = 0.3, 0.1
c = w0**2 - p0**2 / AP
pinf = lam / (lam + S2); sinf = c + pinf**2 / AP
one = np.ones(n)
print(f"{'method':24s} kept  min s(T)   max s(T)   |p-p_inf|  |s-s_inf|  invariant  closed-form")
for name, f, y0 in [("SimSiam", simsiam, np.concatenate([w0*one, p0*one])),
                    ("BYOL, target from 0", byol, np.concatenate([w0*one, p0*one, 0*one])),
                    ("BYOL, target = online", byol, np.concatenate([w0*one, p0*one, w0*one])),
                    ("SimSiam, no predictor", nopred_simsiam, w0*one),
                    ("BYOL, no predictor", nopred_byol, np.concatenate([w0*one, 0*one])),
                    ("SimSiam, no stop-grad", nosg_simsiam, np.concatenate([w0*one, p0*one]))]:
    t, y = run(f, y0, 0.0)
    w = y[:n]; sT = w[:, -1]**2
    kept = int((sT > 1e-6).sum())
    row = f"{name:24s} {kept:4d}  {sT.min():.2e}  {sT.max():.2e}"
    if f in (simsiam, byol):
        pp = y[n:2*n]
        row += f"  {np.abs(pp[:, -1] - pinf).max():.2e}  {np.abs(sT - sinf).max():.2e}  {np.abs(w**2 - pp**2/AP - c).max():.2e}  min(s-c) {(sT - c).min():.2e}"
    elif f is nopred_simsiam:
        i = np.searchsorted(t, [50, 100])
        row += f"  {'':9s}  {'':9s}  {'':9s}  {np.abs(w[:, i] / (w0*np.exp(-S2*t[i])) - 1).max():.2e}"
    elif f is nopred_byol:
        i = np.searchsorted(t, [50, 100]); e = 0.0
        for k in range(n):
            A = np.array([[-(lam[k] + S2), lam[k]], [BETA, -BETA]])
            for j in i:
                e = max(e, abs(w[k, j] / (expm(A * t[j]) @ [w0, 0.0])[0] - 1))
        row += f"  {'':9s}  {'':9s}  {'':9s}  {e:.2e}"
    else:
        bound = w0 * np.exp(-S2 * t)
        use = bound >= 1e-8                      # only where the bound has not underflowed
        ratio = np.abs(w[:, use]) / bound[use]
        row += f"  {'':9s}  {'':9s}  {'':9s}  bound ratio max {ratio.max():.8f} (t <= {t[use].max():.0f})"
    print(row)

# --- C. weight decay, start on the parabola: w0 = p0 = 0.1 (c = 0) ---
print("\nC. weight decay, w0 = p0 = 0.1 (c = 0), T = 2e5")
print(f"{'eta':>7s} {'method':22s} surv fp  pred top-k in-fp =pred  max rel|p-p*+|  invariant  max s collapsed  min s surv")
for eta in (1e-4, 1e-3, 1e-2):
    ok, pm, pp_ = pstar(eta)
    pred = ok & (0.1 > pm)
    for name, f, y0 in [("SimSiam", simsiam, np.concatenate([0.1*one, 0.1*one])),
                        ("BYOL, target from 0", byol, np.concatenate([0.1*one, 0.1*one, 0*one])),
                        ("BYOL, target = online", byol, np.concatenate([0.1*one, 0.1*one, 0.1*one]))]:
        t, y = run(f, y0, eta)
        w, p = y[:n], y[n:2*n]
        sT = w[:, -1]**2; surv = sT > 1e-6
        k = int(surv.sum())
        relp = np.abs(p[surv, -1] / pp_[surv] - 1).max() if k else np.nan
        print(f"{eta:7.0e} {name:22s} {k:4d} {int(ok.sum()):3d} {int(pred.sum()):4d} "
              f"{str(np.array_equal(np.where(surv)[0], np.arange(k))):5s} {str(bool(np.all(ok[surv]))):5s} {str(np.array_equal(surv, pred)):5s} "
              f"{relp:14.2e}  {np.abs(w**2 - p**2/AP).max():.2e}  "
              f"{(sT[~surv].max() if (~surv).any() else np.nan):15.2e}  {(sT[surv].min() if k else np.nan):.2e}")

# --- D. full 64 x 64 matrix dynamics (BYOL, eta = 1e-3): does the predictor align? ---
print("\nD. full matrix BYOL, eta = 1e-3, random start (seed 0), target from 0")
eta = 1e-3
rng = np.random.default_rng(0)
W0 = 0.1 * rng.standard_normal((n, n)) / np.sqrt(n)
Wp0 = 0.1 * rng.standard_normal((n, n)) / np.sqrt(n)
C = W0 @ W0.T - Wp0.T @ Wp0 / AP
td = [0, 10, 100, 1000, 1e4, 2e4]
t, y = run_matrix(np.concatenate([W0.ravel(), Wp0.ravel(), np.zeros(n * n)]), eta, 2e4, td)
print(f"{'t':>7s}  asymmetry  commutator  invariant  directions(s>1e-6)")
for k in range(len(t)):
    W, Wp, Wa = y[:, k].reshape(3, n, n)
    F = W @ X @ W.T
    asym = np.linalg.norm(Wp - Wp.T) / np.linalg.norm(Wp)
    comm = np.linalg.norm(F @ Wp - Wp @ F) / (np.linalg.norm(F) * np.linalg.norm(Wp))
    inv = np.linalg.norm(W @ W.T - Wp.T @ Wp / AP - np.exp(-2 * eta * t[k]) * C) / np.linalg.norm(C)
    print(f"{t[k]:7.0f}  {asym:.2e}   {comm:.2e}    {inv:.2e}   {int((np.linalg.eigvalsh(W @ W.T) > 1e-6).sum())}")
Ul, sv, Vt = np.linalg.svd(W)
m = int((sv**2 > 1e-6).sum())
fe, fu = np.linalg.eigh(F); fu = fu[:, ::-1][:, :m]
corr = min(u @ Wp @ u / np.linalg.norm(Wp @ u) for u in fu.T)
V = Vt[:m]
lh = np.array([v @ X @ v for v in V])
resid = max(np.linalg.norm(X @ v - (v @ X @ v) * v) for v in V)
ok_, pm_, pp_ = (lambda d: (d >= 0, None, (lh + np.sqrt(np.clip(d, 0, None))) / (2 * (lh + S2))))(lh**2 - 4 * eta * (lh + S2))
pp_ = np.where(ok_, pp_, np.nan)
rel = np.abs(sv[:m]**2 / (pp_**2 / AP) - 1).max()
print(f"surviving directions {m} (fixed points exist for {int(pstar(eta)[0].sum())} modes)")
print(f"min Fig.-2 correlation u'Wp u/|Wp u| over surviving eigenvectors of F: {corr:.6f}")
print(f"max |Xv - (v'Xv)v| over surviving input directions v: {resid:.2e}")
print(f"max relative diff of surviving eigenvalues of WW' from p*+(v'Xv)^2: {rel:.2e}")
print("surviving v'Xv:", np.array2string(np.sort(lh)[::-1], precision=6), "| lambda_1..m:", np.array2string(lam[:m], precision=6))

Output:

patches 11942, dim 64; lambda_1 1.000000, lambda_2 0.769866, lambda_63 2.506e-03, |lambda_64| 3.4e-18; X symmetric True; sorted True

A. exactness checks
A1 per-mode RHS at lambda=1 vs Tian et al. Eqs. 11-13, max abs diff: 4.44e-16
A2 matrix ODE from an aligned start vs per-mode ODE (eta=1e-3, t<=200), max abs diff: 2.28e-10

B. eta = 0, w0 = 0.3, p0 = 0.1, c = 0.08, T = 2e5
method                   kept  min s(T)   max s(T)   |p-p_inf|  |s-s_inf|  invariant  closed-form
SimSiam                    64  8.00e-02  1.02e+00  1.11e-16  2.53e-11  2.53e-11  min(s-c) 6.72e-12
BYOL, target from 0        64  8.00e-02  1.02e+00  1.11e-16  5.46e-11  5.46e-11  min(s-c) 5.30e-12
BYOL, target = online      64  8.00e-02  1.02e+00  1.11e-16  5.72e-11  5.72e-11  min(s-c) 6.68e-12
SimSiam, no predictor       0  1.40e-63  1.40e-63                                   2.68e-10
BYOL, no predictor          0  7.85e-100  4.22e-55                                   4.12e-10
SimSiam, no stop-grad       0  2.71e-152  8.05e-63                                   bound ratio max 1.00000000 (t <= 560)

C. weight decay, w0 = p0 = 0.1 (c = 0), T = 2e5
    eta method                 surv fp  pred top-k in-fp =pred  max rel|p-p*+|  invariant  max s collapsed  min s surv
  1e-04 SimSiam                  54  54   54 True  True  True        2.55e-05  1.78e-15         2.04e-10  4.25e-03
  1e-04 BYOL, target from 0      54  54   54 True  True  True        2.56e-05  6.07e-18         2.06e-10  4.25e-03
  1e-04 BYOL, target = online    54  54   54 True  True  True        2.56e-05  5.55e-17         2.06e-10  4.25e-03
  1e-03 SimSiam                  23  25   23 True  True  True        2.22e-16  7.77e-16         3.16e-40  5.45e-02
  1e-03 BYOL, target from 0      23  25   23 True  True  True        2.22e-16  2.22e-16         4.24e-35  5.45e-02
  1e-03 BYOL, target = online    23  25   23 True  True  True        2.22e-16  1.11e-16         4.70e-36  5.45e-02
  1e-02 SimSiam                   5   8    5 True  True  True        2.22e-16  7.77e-16         6.75e-62  5.99e-01
  1e-02 BYOL, target from 0       5   8    5 True  True  True        2.22e-16  2.22e-16         1.01e-51  5.99e-01
  1e-02 BYOL, target = online     5   8    5 True  True  True        2.22e-16  6.78e-21         4.22e-49  5.99e-01

D. full matrix BYOL, eta = 1e-3, random start (seed 0), target from 0
      t  asymmetry  commutator  invariant  directions(s>1e-6)
      0  1.39e+00   1.65e-01    0.00e+00   63
     10  1.37e+00   1.37e-01    1.45e-13   64
    100  4.59e-01   7.42e-03    2.41e-12   63
   1000  1.06e-01   8.17e-04    4.27e-13   63
  10000  1.16e-05   7.79e-11    1.13e-14   19
  20000  5.26e-10   1.65e-15    1.19e-14   19
surviving directions 19 (fixed points exist for 25 modes)
min Fig.-2 correlation u'Wp u/|Wp u| over surviving eigenvectors of F: 1.000000
max |Xv - (v'Xv)v| over surviving input directions v: 9.95e-06
max relative diff of surviving eigenvalues of WW' from p*+(v'Xv)^2: 9.55e-07
surviving v'Xv: [1.       0.769866 0.244514 0.221871 0.160078 0.100166 0.085491 0.076024
 0.059542 0.054627 0.042698 0.040966 0.037535 0.033692 0.031167 0.027705
 0.025537 0.023417 0.020045] | lambda_1..m: [1.       0.769866 0.244514 0.221871 0.160078 0.100166 0.085491 0.076024
 0.059542 0.054627 0.042698 0.040966 0.037535 0.033692 0.031167 0.027705
 0.025537 0.023417 0.020045]

Exactness. Part A shows that the per-mode system is Tian et al.’s Eqs. 11–13 at $\lambda = 1$, up to rounding. It also shows that the full matrix flow, started aligned, stays on the per-mode solution to within $2.3 \times 10^{-10}$. In part B, the predictor-free runs match their closed forms to a few parts in $10^{10}$, and the invariant of Eq. 14 holds to $6 \times 10^{-11}$. In part C, with $c = 0$, the integration preserves $s = p^2$ to $2 \times 10^{-15}$. In part D, Theorem 1’s balancing invariant holds at every printed time to $2.4 \times 10^{-12}$, although that run is neither isotropic nor symmetric. The theorem does not need those assumptions.

Collapse without the predictor or the stop-gradient (part B). Without weight decay, SimSiam and BYOL keep all 64 modes under both target starts. Every $s(T)$ is at least $c$: the printed $\min(s - c)$ is positive in every row. This is Eq. 14 at work. With $\eta = 0$ the conserved $c$ keeps every mode off the collapsed branch, and each predictor eigenvalue ends at $\lambda_k/(\lambda_k + \sigma^2)$, as derived above, and each $s(T)$ equals $c + p_\infty^2/\alpha_p$ to within $6 \times 10^{-11}$. Even the zero-eigenvalue brightness direction keeps $s = c$ while its predictor goes to zero. That mode learned nothing, so “kept” here means only that it did not collapse. All three ablations collapse every mode, to $s(T)$ of $10^{-55}$ or less, from the same start and with the same data. The no-stop-gradient run never exceeds its bound $w_0 e^{-\sigma^2 t}$; its printed ratio is 1 at most, checked over the times where the bound is still above $10^{-8}$. Removing the predictor collapses BYOL even with its slow target; in our reading, the target then only averages a quantity that is already shrinking, as the decaying linear system derived above shows.

Weight decay keeps the strong modes (part C). SimSiam’s survivors are exactly the predicted sets: the 54, 23 and 5 strongest patch directions at $\eta = 10^{-4}$, $10^{-3}$ and $10^{-2}$. Every method’s survivors lie inside the set whose fixed points exist (54, 25 and 8 modes), and the counts fall as $\eta$ rises, as Obs. 4 and 5 predict. At $\eta = 10^{-3}$ and $10^{-2}$, every survivor sits at $p^*_+(\lambda_k)$ to rounding. At $\eta = 10^{-4}$, the check that survivors are within $10^{-6}$ (relative) of $p^*_+$ failed: the largest relative gaps are $2.55 \times 10^{-5}$ for SimSiam and $2.56 \times 10^{-5}$ for both BYOL starts. In our reading this is slow convergence, not a different end point. On the parabola with $\tau \equiv 1$, $\dot p = -(\lambda + \sigma^2)\,p\,(p - p^*_-)(p - p^*_+)$, so the linearized rate at $p^*_+$ is $(\lambda+\sigma^2)\,p^*_+\,(p^*_+ - p^*_-)$. For the weakest survivors at the smallest $\eta$, $p^*_+$ is small and, near the survival threshold, close to $p^*_-$, so that rate is small, and $T = 2\times 10^5$ was not long enough. The script shows no run that tests this further.

We expected the EMA’s slow start to trap more modes (Obs. 10). In this example that did not happen: BYOL with its target starting from zero kept 54, 23 and 5 modes, the same as SimSiam. The hypothesis that it keeps no more than SimSiam held, with equality. Starting the target as a copy of the online network changed nothing either. In the per-mode system, alignment holds by construction from the start, so the EMA’s claimed benefit through $K(t)$ (Obs. 8) cannot appear here.

The surviving modes are the top eigenvectors of the patch distribution, the smooth, low-frequency patterns of a principal component analysis. The weak, high-frequency ones collapse. Their augmentation noise $\sigma^2$ is large relative to their signal $\lambda_k$, so the predictor’s target $\lambda_k/(\lambda_k+\sigma^2)$ is small and weight decay wins. In this linear model, BYOL behaves like a noise-aware principal component analysis whose cut-off is set by $\eta$ and $\sigma^2$. Tian et al.’s supplement makes the related point for additive augmentation noise without weight decay: the learned modes are those of $X(X + X')^{-1}$, and dimensions with larger noise variance are attenuated (Appendix D.2).

Alignment from a random start (part D). The full matrix flow starts with an asymmetric predictor (asymmetry 1.39) that does not commute with its input correlation (commutator 0.165). By $t = 2\times 10^4$, the predictor is symmetric to $5 \times 10^{-10}$, the commutator is $2 \times 10^{-15}$, and the Figure 2 correlation is 1.000000 for every surviving direction. Nineteen directions survive, within the 25 whose fixed points exist. They are eigenvectors of $X$ to within $10^{-5}$, they are exactly the 19 strongest ($v^\top X v$ equals $\lambda_1, \ldots, \lambda_{19}$ at the printed six decimals), and their eigenvalues match $p^*_+(\lambda)^2$ to $10^{-6}$. Every planned check on this run held. The matrix run keeps fewer modes than the per-mode prediction from $p_0 = 0.1$ (23 at this $\eta$). In our reading, its small random start puts more weak directions below $p^*_-$, and that basin decides. The flow thus finds the alignment that Theorem 3 proves only under isotropy, symmetry and a proportional target, and DirectPred’s eigenbasis rule describes where it ends up.

Relation to the paper’s figures. The example matches Tian et al.’s figures only qualitatively and reproduces none of them. Their Figures 3–5 simulate isotropic data with settings the captions do not fully give, and their Figure 2 comes from training a ResNet-18 on STL-10. What carries over is the structure: the collapsed and non-collapsed branches and the basin boundary $p^*_-$ (Figures 3–4), weight decay trading surviving modes for stability (Figure 5), and the predictor becoming symmetric and aligned (Figure 2).

Inference

After pre-training, the backbone serves as a feature extractor. In Tian et al.’s and Chen and He’s ImageNet experiments, the linear classifier is trained on the frozen backbone features. Chen and He also fine-tune the pretrained ResNet-50 end to end for detection and segmentation (§6.1). In our reading, the projector, predictor and target exist only for training. Evaluating Representations covers the linear, k-NN and fine-tuning protocols.

Performance

These are linear-probe results on ImageNet, each as reported by a readable source:

• Under a common protocol. Chen and He reproduced BYOL with a ResNet-50 and two $224 \times 224$ views, batch size 4,096 and single-crop evaluation (Table 4). It reached 66.5%, 70.6%, 73.2% and 74.3% after 100, 200, 400 and 800 epochs. SimSiam, at batch size 256, reached 68.1%, 70.0%, 70.8% and 71.3%. BYOL’s reproduction was the best of the five methods compared at 200 epochs and beyond. SimSiam was the best at 100.
• Transfer. In Chen and He’s transfer comparison (Table 5, 200-epoch pre-training, fine-tuned end to end), BYOL’s COCO detection AP on val2017 (Mask R-CNN, C4 backbone, 1× schedule) was 37.9, the same as SimCLR’s and SimSiam’s base recipe, against 39.2 for MoCo v2. On VOC 2007 test detection (Faster R-CNN, trained on VOC 2007 trainval), its AP50 was 77.1, equal to MoCo v2’s and above SimSiam’s 75.5 (base recipe).
• Vision Transformers. MoCo v3 (Chen, Xie & He, 2021, Table 4) trained ViTs for 300 epochs. BYOL reached 71.0% with ViT-S/16 and 73.9% with ViT-B/16, against 72.5% and 76.5% for MoCo v3; with ResNet-50 at 800 epochs, where the other frameworks’ numbers come from SimSiam’s improved implementation, BYOL’s 74.3% was the highest of the four. MoCo v3 also names BYOL’s 79.6% with a ResNet-200 (2× width), next to SimCLR v2’s 79.8%, as the previous best results under linear probing (§6.3).
• Further comparisons. Self-Supervised Learning lists BYOL’s ResNet-50 and large-ResNet results from the BYOL paper, together with the drop under crop-only augmentation. The DINO paper compares DINO with BYOL on ResNet-50 and ViT-S. These are not repeated here.

Strengths

• No negatives. BYOL needs neither the large batches nor the memory queue that contrastive methods use to supply negative pairs (Tian et al., 2021, §1). The SimSiam variant also works at ordinary batch sizes: 66.1–68.1% from 64 to 2,048 images (Chen & He, 2021, Table 2).
• Competitive representations. In Chen and He’s common-protocol reproduction, BYOL had the highest 800-epoch linear accuracy of the methods compared. Chen and He report that all the methods in their transfer comparison, BYOL included, surpass or match ImageNet supervised pre-training on every task they tested (Table 5).
• An analysable mechanism. Because its loss is a plain regression onto a stop-gradient target, its dynamics can be analysed in a linear model. That analysis led Tian et al. to DirectPred, which sets the predictor in closed form and matched BYOL’s default predictor on ImageNet at 300 epochs.

Limitations

• Collapse is avoided, not excluded. Collapsed solutions exist for BYOL’s loss. Without the stop-gradient or the predictor, SimSiam collapses (Chen & He, 2021), and the linear theory proves the same collapse in its setting. Chen and He call SimSiam’s non-collapsing behaviour an empirical observation.
• Sensitivity to the predictor and weight decay. In Tian et al.’s experiments, outcomes hinged on settings that a contrastive loss does not depend on: the predictor’s learning rate, its symmetry and the weight decay. Without the EMA, a symmetric predictor crashed. Without weight decay, accuracy fell. Too strong a weight decay collapses modes in the theory, as the example shows.
• Instability with ViTs. MoCo v3 found that BYOL training of ViTs can be unstable. A fixed random patch projection improved BYOL’s ViT-B/16 accuracy at 100 epochs from 69.7% to 71.0% (Chen, Xie & He, 2021, Figure 6).
• Multi-crop. As the DINO paper reports, multi-crop training lowered BYOL’s accuracy where it raised DINO’s, SwAV’s and MoCo v2’s.
• Theory beyond linear networks. Tian et al. leave the two-layer predictor and the asymmetric predictor without a formal analysis.

Historical Importance

BYOL was published at NeurIPS in 2020. Chen and He and Tian et al. both cite it as the method that learned strong representations from positive pairs alone. SimSiam (CVPR 2021) reduced it to a stop-gradient and a predictor with shared weights, and identified the stop-gradient as the component that prevents collapse. Tian, Chen and Ganguli (ICML 2021) gave the first analytic account of the dynamics, in their words, and turned it into DirectPred. Later methods kept pieces of the design:

• MoCo v3 adds a BYOL-style prediction head to a contrastive method (Chen, Xie & He, 2021; MoCo covers the version).
• DINO keeps the momentum teacher and the stop-gradient, but replaces the predictor with centering and sharpening (DINO).
• Consistency models: their authors relate their EMA target with a stop-gradient to momentum-based methods such as BYOL and MoCo (Consistency Model).

Variants

• SimSiam (Chen & He, 2021): BYOL without the momentum encoder. Weights are shared directly, a prediction MLP sits on one branch and a stop-gradient on the other, with a symmetrized negative-cosine loss and SGD, at batch size 512 by default and 256 in its comparison with other methods. It shows that the moving average is not needed to prevent collapse.
• BYOL without EMA (BYOL arXiv v3, as Chen and He report): 66.9% at 300 epochs when the target equals the online network and the predictor’s learning rate is raised tenfold.
• Linear and symmetric predictors (Tian et al., 2021): one-layer predictors, constrained to be symmetric or not. Their behaviour with and without the EMA is the evidence for the eigenspace analysis.
• DirectPred (Tian et al., 2021): the predictor is set at each step from the eigendecomposition of its input correlation, $p_j = \sqrt{s_j} + \epsilon \max_j s_j$. A variant mixing this with gradient steps (DirectPred with freq = 5) reached 80.28 ± 0.62 on STL-10 after 500 epochs (Table 8).
• Prediction heads in contrastive methods: MoCo v3 uses a BYOL-style prediction MLP. Removing it lowered ViT-B linear accuracy from 76.5% to 75.5%. The authors note that a contrastive method, unlike negative-free ones, does not need the predictor to work (Chen, Xie & He, 2021, §6.2).
• DINO (DINO): a momentum teacher with softmax outputs, centering and sharpening, and no predictor or negatives.