Yes, there is a way to improve performance... Like Dan mentioned, don't use ArcGIS (directly at least).
Here is another approach leveraging Python dictionaries. Note that this code only works for planar coordinates (UTM, State Plane, etc.) so if your data is in Geographic coordinates you need to use another distance function! BTW: This code assumes all fields are fully populated with valid values - if not you'll have to handle those exceptions.
Numpy has some nice ways to do all this stuff too - I usually prefer to roll my own though .
import math
def getDist(x1, y1, x2, y2):
return math.sqrt((x1-x2)**2 + (y1-y2)**2)
sitesDict = {r[0],r[1]:r[2:] for r in arcpy.da.SearchCursor(sitesFC, ["SHAPE@X","SHAPE@Y","RESTAURANT_ID"])}
customersDict = {r[0],r[1]:r[2:] for r in arcpy.da.SearchCursor(custFC, ["SHAPE@X","SHAPE@Y","INCOME"])}
for x1, y1 in sitesDict:
incomeList = []
for x2, y2, in customersDict:
if getDist(x1, y1, x2, y2) <= 1000: #or whatever map unit theashold...
income = customersDict[(x2, y2)][0])
if income > 0:
incomeList.append(customersDict[(x2, y2)][0])
#blah
incomeListSum = sum(incomeList)
incomeListLen = len(incomeList)
print ("Restaurant #" + str(sitesDict[(x1,y1)][0] + " has " + str(incomeListLen) + " customers - average income is " = str(incomeListSum / float(incomeListLen)))