← Back to the problem

Sweden SAIO · Correction by IOAI Training

Colorize

This page gives the full approach. Read it once you have searched.

Metric
MSE on sampled RGB values
Evaluated on
held-out validation (official test labels are private)
Our score
MSE 527.2 (100/100 by the formula)
Expected level
full marks at MSE <= 567, zero at >= 666

Expected level reached

We solved the problem ourselves: downloaded the data, ran the experiments and measured the score below. The statement and data belong to their organisers.

The problem in brief

You receive 2,871 photos of 128x128 pixels (from the Yahoo Flickr Creative Commons collection) together with their grayscale versions, and 1,547 grayscale test photos whose colours were lost. You must predict the colour of every 131st pixel of the test images (counting pixels across all images in a fixed order) and print the R, G, B values as integers; the printed list is the submission. The hidden original colours are scored by mean squared error. The task comes from the final of the Swedish AI Olympiad 2026 (saio.kattis.com/problems/colorize); the test colours are private, so we score on 15 % of the training pairs held out from training (446 images).

How you are scored

With gig_i the predicted and tit_i the true value over all evaluated R, G and B entries,

MSE=1N∑i=1N(gi−ti)2,Score=100⋅clip⁡(666−MSE99, 0, 1).\text{MSE} = \frac{1}{N}\sum_{i=1}^{N}(g_i - t_i)^2 , \qquad \text{Score} = 100\cdot\operatorname{clip}\left(\frac{666 - \text{MSE}}{99},\,0,\,1\right).

So an MSE of 666 or more earns nothing and 567 or less earns full marks. The organisers' baseline, which copies the gray value into all three channels, scores 696 on our validation images: below zero points. The target is to cut the error of "no colour at all" by about a fifth.

First look at the data

  • The gray image is exactly the luma. On every training pair, the gray value equals Y=0.299R+0.587G+0.114BY = 0.299R + 0.587G + 0.114B of the colour image (mean absolute difference 0.2 gray levels, JPEG noise). That is the standard conversion, and it means the brightness of each pixel is known: only its colour (hue and saturation) must be guessed.
  • Colour has a strong average. The mean colour offset from gray over all training pixels is (+8.1,−1.8,−11.8)(+8.1, -1.8, -11.8) for (R, G, B): photos are on average warm, slightly yellow-orange. Adding this constant shift lowers the MSE from 696 to 616.5 (50 points). Making the shift depend on the gray level (the mean offset of all training pixels with that gray value, a 256-entry lookup table) gives 599 (67 points): dark pixels tend to be brownish, bright ones bluish-white like the sky.
  • The error is dominated by colourful images. Per-image MSE of our final model ranges from below 70 (gray buildings) to over 2,500 (a cyan illuminated clock). A few vivid, unpredictable images weigh as much as dozens of easy ones.

The key idea

Two observations shape the solution.

1. Predict only what is unknown, and keep the luma exact. Write each colour as the gray value plus two chroma components, Cr=R−YC_r = R - Y and Cb=B−YC_b = B - Y. Given YY (the input) and the two chroma values, the green channel follows from the luma equation:

R=Y+Cr,B=Y+Cb,G=Y−0.299 Cr+0.114 Cb0.587.R = Y + C_r, \qquad B = Y + C_b, \qquad G = Y - \frac{0.299\,C_r + 0.114\,C_b}{0.587}.

The network outputs two numbers per pixel instead of three, and every prediction has exactly the right brightness, so no capacity is wasted relearning it.

2. Under MSE, the best guess is the conditional mean, which is a muted colour. For a pixel whose colour is uncertain (a sweater could be red, blue or green), the prediction minimising the expected squared error is the average of the possible colours, which is close to gray-brown. A model trained with the MSE loss learns exactly that: it paints confident colours where the content is unambiguous (grass, sky, skin, wood) and stays cautious elsewhere. The colourisations look dull, but that is what the metric rewards. A model that produces vivid, "realistic" colours (for instance with an adversarial loss) would be punished hard whenever it picks the wrong hue: predicting red for a blue sweater costs far more than predicting gray. Choose the loss that is the metric.

Context is what makes the conditional mean sharp: knowing that a region is grass, sky or a face requires recognising objects, so we use a convolutional network with an ImageNet-pretrained encoder.

A first baseline and the model

The lookup table by gray level (MSE 599) is the natural non-neural baseline. The model is a U-Net:

  • Encoder: ResNet-18 with ImageNet weights, applied to the gray image copied into three channels (normalised like ImageNet), giving feature maps at 1/2, 1/4, 1/8, 1/16 and 1/32 of the resolution.
  • Decoder: five upsampling stages, each concatenating the encoder features of the same resolution (skip connections) and the gray image itself at full resolution, then two 3x3 convolutions; a final 1x1 convolution outputs CrC_r and CbC_b.
  • Training: loss = MSE between the reconstructed RGB image and the original (the metric itself), AdamW with a one-cycle schedule (peak learning rate 3×10−43\times10^{-4} for the pretrained encoder, 10−310^{-3} for the decoder), batch 32, random horizontal flips; at test time the prediction is averaged with that of the mirrored image. One epoch over the 2,425 training images takes about 10 s on an Apple laptop GPU.

The skip connections matter: colour must be placed precisely at object boundaries, and the deep features alone (4x4 pixels at the bottleneck) know what the objects are but not exactly where their edges lie.

Validation images sorted from easy to hard. The model colours grass, skin, wood and sky correctly but leaves ambiguous objects (clothes, the illuminated clock) in safe, muted tones, as the MSE loss demands.

The full solution

  1. Read the training pairs and the test images in sorted filename order.
  2. Train the U-Net for 15 epochs on 85 % of the pairs and check the validation MSE; then retrain the same configuration on all 2,871 pairs.
  3. Predict the test images (with mirror averaging), clip to [0,255][0, 255] and round to integers.
  4. Write out.py, which prints the R G B values of every 131st pixel in the official order (global counter over images, xx in the outer loop, yy in the inner loop), exactly like the organisers' baseline: 193,481 pixels, 580,443 integers.
cr, cb = chroma[:, 0], chroma[:, 1]                      # network output
rgb = torch.stack([g + cr, g - (0.299 * cr + 0.114 * cb) / 0.587, g + cb], 1)
loss = F.mse_loss(rgb, true_rgb)                          # the competition metric itself

Results

Validation MSE on 446 held-out training pairs, over all pixels (the official every-131st-pixel sample gives nearly the same value, for example 528.7 against 527.2 for the final model).

approachvalidation MSEscore
copy gray into R, G, B (organisers' baseline)695.90
constant colour shift616.550.0
mean shift per gray level (lookup table)598.967.8
U-Net with randomly initialised encoder, 30 epochs567.999.1
U-Net with ImageNet ResNet-18 encoder, 30 epochs533.1100
U-Net with ImageNet ResNet-18 encoder, 15 epochs (final)527.2100

The ladder of approaches. The green band is where points are earned; everything below MSE 567 is full marks.

The pretrained encoder is worth 35 points of MSE over training from scratch: 2,871 images are too few to learn object recognition, but ImageNet features already know grass from sky. With 30 epochs the training MSE (about 200) falls far below the validation MSE: the model starts memorising the training images, which is why the shorter schedule generalises slightly better. With a margin of 40 below the full-marks threshold, the result should be robust to the difference between our validation images and the hidden test set.

What to remember

  • Decompose the target into what is known and what is not. The gray input is the luma, so predict only the two chroma channels and reconstruct RGB exactly.
  • Train on the metric. Under squared error the optimal prediction is the conditional mean, so muted colours are not a failure but the right answer to an MSE question; vivid generative colourisation would score worse.
  • Climb the ladder of baselines. Constant shift, then a lookup table by gray level, then a network: each step tells you how much structure is left to exploit.
  • Pretrained encoders transfer to dense prediction. A U-Net with an ImageNet backbone and skip connections combines "what is this object" with "where exactly is its edge".
  • Reproduce the output format exactly. The global pixel counter across images and the xx-before-yy loop order are easy to get wrong; reuse the organisers' loop.
Full code (solution.py)
"""SAIO (Sweden) final 2026, "Colorize": predict RGB from grayscale, scored by MSE on sampled pixels.

Files: train/colored/*.jpg, train/blackwhite/*.jpg, test/*.jpg (Kattis attachments, 128x128 px);
ImageNet ResNet-18 weights as safetensors (timm/resnet18.tv_in1k) in r18.safetensors.
usage: python solution.py <pretrained 0|1> <epochs>        e.g. python solution.py 1 40
Holds out 15 % of the training pairs for validation, then retrains on everything and writes out.py
(the Kattis submission: a program that prints the predicted pixels in the official order).
"""
import os, sys, time
import numpy as np, torch, torch.nn as nn, torch.nn.functional as F
from PIL import Image
from safetensors.torch import load_file
torch.set_num_threads(2)
DEV = torch.device("mps" if torch.backends.mps.is_available() else "cpu")
PRE, EPOCHS = sys.argv[1] == "1", int(sys.argv[2])

def load(folder, mode):
    files = sorted(os.listdir(folder))
    return files, np.stack([np.asarray(Image.open(f"{folder}/{f}").convert(mode)) for f in files])
names, C = load("train/colored", "RGB")
_, G = load("train/blackwhite", "L")
test_names, T = load("test", "L")
va = np.random.default_rng(0).random(len(C)) < 0.15

def to_rgb(g, chroma):
    """The gray image is the luma Y = 0.299 R + 0.587 G + 0.114 B. Predict only Cr = R - Y and Cb = B - Y;
    then G follows exactly, so the luma of the prediction always equals the input."""
    cr, cb = chroma[:, 0], chroma[:, 1]
    return torch.stack([g + cr, g - (0.299 * cr + 0.114 * cb) / 0.587, g + cb], 1)

class UNet(nn.Module):
    """ResNet-18 encoder (ImageNet weights, optional) + light upsampling decoder with skip connections."""
    def __init__(self, pretrained):
        super().__init__()
        import torchvision
        r = torchvision.models.resnet18()
        if pretrained:
            r.load_state_dict(load_file("r18.safetensors"))
        self.stem = nn.Sequential(r.conv1, r.bn1, r.relu)                    # 64 ch, 1/2
        self.l1, self.l2, self.l3, self.l4 = nn.Sequential(r.maxpool, r.layer1), r.layer2, r.layer3, r.layer4
        def up(cin, cskip, cout):
            return nn.Sequential(nn.Conv2d(cin + cskip, cout, 3, padding=1), nn.BatchNorm2d(cout), nn.ReLU(inplace=True),
                                 nn.Conv2d(cout, cout, 3, padding=1), nn.BatchNorm2d(cout), nn.ReLU(inplace=True))
        self.d4, self.d3, self.d2 = up(512, 256, 256), up(256, 128, 128), up(128, 64, 64)
        self.d1, self.d0 = up(64, 64, 48), up(48, 1, 32)
        self.out = nn.Conv2d(32, 2, 1)
        self.register_buffer("mean", torch.tensor([0.485, 0.456, 0.406]).view(1, 3, 1, 1))
        self.register_buffer("std", torch.tensor([0.229, 0.224, 0.225]).view(1, 3, 1, 1))
    def forward(self, g):                                                      # g: (B, 1, H, W) in [0, 255]
        x = ((g / 255).repeat(1, 3, 1, 1) - self.mean) / self.std
        s0 = self.stem(x); s1 = self.l1(s0); s2 = self.l2(s1); s3 = self.l3(s2); s4 = self.l4(s3)
        cat_up = lambda a, b: torch.cat([F.interpolate(a, size=b.shape[-2:], mode="bilinear", align_corners=False), b], 1)
        y = self.d4(cat_up(s4, s3)); y = self.d3(cat_up(y, s2)); y = self.d2(cat_up(y, s1))
        y = self.d1(cat_up(y, s0)); y = self.d0(cat_up(y, g / 255))
        return 40 * self.out(y)                                               # chroma (Cr, Cb), in grey levels

def run(train_idx, epochs, bs=32):
    torch.manual_seed(0)
    model = UNet(PRE).to(DEV)
    enc = [p for n, p in model.named_parameters() if n.startswith(("stem", "l1", "l2", "l3", "l4"))]
    dec = [p for n, p in model.named_parameters() if not n.startswith(("stem", "l1", "l2", "l3", "l4"))]
    opt = torch.optim.AdamW([{"params": enc, "lr": 3e-4 if PRE else 1e-3}, {"params": dec, "lr": 1e-3}], weight_decay=1e-4)
    steps = epochs * (len(train_idx) // bs)
    sch = torch.optim.lr_scheduler.OneCycleLR(opt, [3e-4 if PRE else 1e-3, 1e-3], total_steps=steps, pct_start=0.1)
    Gt = torch.tensor(G[train_idx], dtype=torch.float32)[:, None]; Ct = torch.tensor(C[train_idx], dtype=torch.float32).permute(0, 3, 1, 2)
    t0 = time.time()
    for ep in range(epochs):
        model.train(); perm = torch.randperm(len(Gt))
        for k in range(0, len(perm) - bs + 1, bs):
            b = perm[k:k + bs]; g, c = Gt[b].to(DEV), Ct[b].to(DEV)
            if torch.rand(1) < 0.5: g, c = g.flip(3), c.flip(3)                 # mirror augmentation
            loss = F.mse_loss(to_rgb(g[:, 0], model(g)), c)                    # the metric itself
            opt.zero_grad(); loss.backward(); opt.step(); sch.step()
        if ep % 5 == 4 or ep == epochs - 1:
            print(f"  epoch {ep + 1}: train MSE {loss.item():.1f}, {time.time() - t0:.0f}s", flush=True)
    return model

@torch.no_grad()
def predict(model, gray, tta=True):
    model.eval(); out = []
    for k in range(0, len(gray), 64):
        g = torch.tensor(gray[k:k + 64], dtype=torch.float32)[:, None].to(DEV)
        ch = model(g)
        if tta: ch = (ch + model(g.flip(3)).flip(3)) / 2
        out.append(to_rgb(g[:, 0], ch).clamp(0, 255).round().permute(0, 2, 3, 1).cpu().numpy())
    return np.concatenate(out)

def sampled_mse(pred, true):
    """Official sampling: global pixel counter over images in sorted order, x outer, y inner, every 131st."""
    p, t = pred.transpose(0, 2, 1, 3).reshape(-1, 3), true.transpose(0, 2, 1, 3).reshape(-1, 3).astype(float)
    keep = (np.arange(1, len(p) + 1) % 131) == 0
    return ((p[keep] - t[keep]) ** 2).mean()

if __name__ == "__main__":
    model = run(np.where(~va)[0], EPOCHS)
    P = predict(model, G[va])
    m_all, m_s = ((P - C[va].astype(float)) ** 2).mean(), sampled_mse(P, C[va])
    print(f"pretrained={PRE} epochs={EPOCHS}: validation MSE all pixels {m_all:.1f}, sampled {m_s:.1f}, "
          f"score {100 * max(0, min(1, (666 - m_all) / 99)):.1f}/100")
    np.save(f"valpred_pre{int(PRE)}_{EPOCHS}.npy", P)
    if len(sys.argv) > 3 and sys.argv[3] == "final":                         # retrain on all pairs, write out.py
        model = run(np.arange(len(C)), EPOCHS)
        P = predict(model, T).astype(int)
        with open("out.py", "w") as f:
            f.write('print("""\n')
            counter = 0
            for img in P:
                for x in range(img.shape[1]):
                    for y in range(img.shape[0]):
                        counter += 1
                        if counter % 131 == 0:
                            f.write(f"{img[y, x, 0]} {img[y, x, 1]} {img[y, x, 2]} ")
                f.write("\n")
            f.write('""")\n')