To speed this up further, you should use the [`numpy` interface of `PIL`][1] (if you are not yet using `PIL`, you should, for this reason):

    from PIL import Image
    import numpy as np
    
    # `int` important because otherwise it might wrap around when subtracting
    optimal_data = np.asarray(Image.open("optimal.png"), dtype=int) 
    new = np.random.randint(0, 256, optimal_data.shape)
    
    def fitness(optimal_data, new):
        return np.sqrt(((optimal_data - new)**2).sum(axis=-1)).sum()

This takes only 258 ms ± 2.21 ms for a 2424 x 2424 pixel image on my machine, while the [function by @TimCPogue][2] takes 9.93 s ± 465 ms with the same images.

Note that the array has the shape `(width, height, channels)`, where `channels` is usually `4` (red, green, blue, alpha), not `3` like your code assumes. If you want to disregard differences in alpha, either set the alpha channel of the new image to the one of the optimal data (`new[:,:,-1] = optimal_data[:,:,-1]`), or slice in the fitness (`optimal_data[...,:-1] - new[...,:-1]`)
.

  [1]: https://stackoverflow.com/a/44606972/4042267
  [2]: https://codereview.stackexchange.com/a/79028/98493