multiple random cordinates in matrix, optimization. python

the code below creates a matrix of zeros (line 06), then takes random "start points" (line 36) that move randomly in x and y (line 10), each start point have their own number(1,2,3,4,5,...,number_fo_start_points), they can't go over the path of other start_point (line 29), the list of start_points enter a while that has a limit of how many "steps" each start_poitn does (line 38) (the limit is defined by other external conditions).

My problem is that I need huge matrixes, about 1.000.000 x 1.000.000. but this code will take for ever doing that, is there a way to optimize this code?? is there a way to move all the start_points at the same time? do I need threading? can I change the same matrix with two different threads?

01    import random
02    import numpy
03    length_matrix = 1000
04    width_matrix = 1000
05    number_of_start_points = 400
06    matriz = numpy.zeros([length_matrix, width_matrix])
07    def direction(): 
08        direction = random.randint(1, 4)
09        return direction
10    def move(x, y):  
11        direction = direction()
12        if direction == 1:
13           return x + 1, y
14        if direction == 2:
15           return x - 1, y
16        if direction == 3:
17           return x, y + 1
18        if direction == 4:
19           return x, y - 1
20            
21    class start_point: 
22        def __init__(self):
23            self.xrand = random.randint(0, LARGOIMG)
24            self.yrand = random.randint(0, ANCHOIMG)
25    
26        def move(self, x, y, start_point_index): 
27            try:
28               position = move(x, y) 
29               if matrix[position] == 0 or matrix[position] == start_point_index + 1:  
30                  matrix[position] = start_point_index + 1
31                  return position 
32               return x, y  
33            except:
34               return x, y 
35    start_points = [ ]
36    for i in range(number_of_start_points):
37        nombres.append(start_point())
38    limit = (((length_matrix * width_matrix * 100000) / (360 * 540))*12)/number_of_start_points
39    while flag < limite:
40          for i in start_points:
41             new_position = (i.move(i.xrand, i.yrand, star_point.index(i)))
42             i.xrand, i.yrand = new_position
43          flag = flag + 1
44    numpy.save("matrix001",matrix)