How to quickly assign an individual material to color thousands of cubes at the same time with python API in Blender 2.8?
One mesh, object linked materials.

You are trying to make same linked mesh have different materials on a different object. Using object linked materials can use only the one mesh, and N x N materials and objects.

import bpy
context = bpy.context
scene = context.scene
count = 32

bpy.ops.mesh.primitive_cube_add(location=(0, 0, 0))
ob = context.object
dx, dy = ob.dimensions.xy
bpy.ops.object.material_slot_add()
ob.material_slots[0].link = 'OBJECT'

mat = bpy.data.materials.new("Mat")
ob.material_slots[0].material = mat
for i in range(1, count * count):
    o = ob.copy()
    o.material_slots[0].material = mat.copy()
    # set some mat props
    o.location = (dx * (i % count), dy * (i // count), 0)
    scene.collection.objects.link(o)


One material could be used here instead by mapping texture coordinates via the object.



Method 2: 1 material, N x N objects and meshes.

Grid, UVMap image texture, edge split, solidify, separate by loose parts

Using the same methods as outlined to you here https://blender.stackexchange.com/a/131765/15543

Material is an image texture mapped to the UV. This way we only have one material. The UV can be mapped to 1 pixel squares.  

 
Result of script. Image is 128x128, grid was split 32 x 32


So each block has a top with 4x4 pixels.

import bpy
context = bpy.context

bpy.ops.mesh.primitive_grid_add(
        size=10,
        x_subdivisions=33,
        y_subdivisions=33,
        calc_uvs=True)

ob = context.object
me = ob.data
me.materials.append(bpy.data.materials["Material"])
sfy = ob.modifiers.new(name="Solidify", type='SOLIDIFY')
sfy.thickness = 1 / 3
bpy.ops.object.mode_set(mode='EDIT')
bpy.ops.mesh.edge_split()
bpy.ops.mesh.separate(type='LOOSE')
bpy.ops.object.mode_set()

bpy.ops.object.origin_set(type='ORIGIN_GEOMETRY')


Creating 128 x 128 objects and as many materials too is going to be slooooooooooooooooow.

Can be optimized for speed somewhat. Not sure it's worth it for zero feedback.



Method 3: 1 object, 1 mesh, N x N materials.

Single cube, applied array modifiers, materials assigned to each part

Similarly to this answer https://blender.stackexchange.com/a/61254/15543 can create an array mesh very quickly.

After the array modifiers are applied the faces are still ordered such that each 6 belongs to a cube.  A new material is added for each and assigned.  The row and column of the parts can be obtained using integer div and modulus and hence can be mapped to pixel colour.

For example sake I've only made a 32 x 32 array. Note the offset is at 1.5, change to 1.0 to have touching cubes. 

To make each a single object, separate by loose parts as above. This is the rate-determining-step.  As one mesh, individual parts can be animated via hooks, armature, shapekeys.

import bpy
from bpy import context
# add a cube
bpy.ops.mesh.primitive_cube_add()
ob = context.active_object
me = ob.data
# add two array modifiers
count = 32
arrays = {"x": (1.5, 0, 0),
          "y": (0, 1.5, 0)}

for axis, displace in arrays.items():          
    mod = ob.modifiers.new(axis, 'ARRAY')
    mod.count = count
    mod.relative_offset_displace = displace
    bpy.ops.object.modifier_apply(modifier=mod.name)

def chunks(l, n):
    for i in range(0, len(l), n):
        yield l[i:i+n]

material_index = 0
for cube in chunks(me.polygons, 6):
    mat = bpy.data.materials.new("Mat%d" % material_index)

    mat.diffuse_color[0] = material_index / 100
    me.materials.append(mat)
    for f in cube:

        f.material_index = material_index
    material_index += 1


Alternatively could use a UV map here also, for 1 object, 1 mesh (with NxN parts) and 1 material.