To speed this up further, you should use the numpy interface of PIL (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 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])  .

To speed this up further, you should use the numpy interface of PIL (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 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]).

For some more readability and the possibility to use a different norm in the future (albeit at the cost of about 30% speed), you could make the norm to use a parameter and use np.linalg.norm, as suggested in the comments by @GarethReese:

def fitness(optimal_data, new, norm=np.linalg.norm):
    return norm(optimal_data - new, axis=-1).sum()