If you can load the two rasters into GDAL and extract their geotransforms, this is easy to do without warping the rasters in any way (such as resampling):
def pixToGeoToPix(x1_index, y1_index, gt1, gt2):
'''
x1_index: x index of pixel in first image
y1_index: y index of pixel in first image
gt1: geotransform tuple of first image
gt2: geotransform tuple of second image
'''
geo_x = gt1[0] + gt1[1]*x1_index
geo_y = gt1[3] + gt1[5]*y1_index
x2_index = int((geo_x-gt2[0])/gt2[1])
y2_index = int((gt2[3]-geo_y)/gt2[5])
return x2_index, y2_index
This function uses the geotransforms of both images to map a pixel index of one image to the equivalent pixel index in the second image by using the pixel's geo-coordinate to translate the pixel index between images. Once you have the matching indices you can easily use numpy to compare the two pixels.
Note that this doesn't necessarily require GDAL. The only information you need is the top left x/y coordinate and x/y spatial resolution of both images, and there are many ways to extract those values.