计算速度慢：itertools.product可能是罪魁祸首吗？

数值积分所花费的时间比我预期的长。我想知道我在网格上实现迭代的方式是否可能是一个促成因素。我的代码如下所示：

import numpy as np
import itertools as it

U = np.linspace(0, 2*np.pi)
V = np.linspace(0, np.pi)

for (u, v) in it.product(U,V):
    # values = computation on each grid point, does not call any outside functions
    # solution = sum(values)
return solution

我省略了计算，因为它们很长，我的问题特别是关于我在参数空间（u，v）上实现计算的方式。我知道诸如以下的替代方法numpy.meshgrid：但是，所有这些似乎都创建了（非常大）矩阵的实例，我猜想将它们存储在内存中会使速度变慢。

有没有其他选择it.product可以加快我的程序的速度，还是我应该在其他地方寻找瓶颈？

编辑：这是有问题的for循环（以查看是否可以向量化）。

import random  
import numpy as np  
import itertools as it 

##########################################################################
# Initialize the inputs with random (to save space)
##########################################################################
mat1 = np.array([[random.random() for i in range(3)] for i in range(3)])
mat2 = np.array([[random.random() for i in range(3)] for i in range(3)]) 
a1, a2, a3 = np.array([random.random() for i in range(3)]) 
plane_normal = np.array([random.random() for i in range(3)])  
plane_point = np.array([random.random() for i in range(3)])  
d = np.dot(plane_normal, plane_point)  
truthval = True

##########################################################################
# Initialize the loop
##########################################################################
N = 100 
U = np.linspace(0, 2*np.pi, N + 1, endpoint = False) 
V = np.linspace(0, np.pi, N + 1, endpoint = False) 
U = U[1:N+1] V = V[1:N+1]

Vsum = 0
Usum = 0

##########################################################################
# The for loops starts here
##########################################################################   
for (u, v) in it.product(U,V):

    cart_point = np.array([a1*np.cos(u)*np.sin(v), 
                           a2*np.sin(u)*np.sin(v), 
                           a3*np.cos(v)])

    surf_normal = np.array(
            [2*x / a**2 for (x, a) in zip(cart_point, [a1,a2,a3])])


    differential_area = \
        np.sqrt((a1*a2*np.cos(v)*np.sin(v))**2 + \
        a3**2*np.sin(v)**4 * \
        ((a2*np.cos(u))**2 + (a1*np.sin(u))**2)) * \
        (np.pi**2 / (2*N**2)) 


    if (np.dot(plane_normal, cart_point) - d > 0) == truthval:
        perp_normal = plane_normal
        f = np.dot(np.dot(mat2, surf_normal), perp_normal)
        Vsum += f*differential_area
    else:
        perp_normal = - plane_normal
        f = np.dot(np.dot(mat2, surf_normal), perp_normal)
        Usum += f*differential_area

integral = abs(Vsum) + abs(Usum)