NumPy에 루프 종속성이있는 for 벡터를 벡터화하는 중첩

Question

4 면체에 일련의 격자 점을 생성하는 다음 함수가 있습니다.

def tet_grid(n): 

    xv = np.array([ 
     [-1.,-1.,-1.], 
     [ 1.,-1.,-1.], 
     [-1., 1.,-1.], 
     [-1.,-1., 1.], 
     ]) 

    nsize = int((n+1)*(n+2)*(n+3)/6); 
    xg = np.zeros((nsize,3)) 
    p = 0 

    for i in range (0, n + 1): 
     for j in range (0, n + 1 - i): 
      for k in range (0, n + 1 - i - j): 
       l = n - i - j - k 
       xg[p,0]=(i * xv[0,0] + j * xv[1,0] + k * xv[2,0] + l * xv[3,0])/n 
       xg[p,1]=(i * xv[0,1] + j * xv[1,1] + k * xv[2,1] + l * xv[3,1])/n 
       xg[p,2]=(i * xv[0,2] + j * xv[1,2] + k * xv[2,2] + l * xv[3,2])/n 
       p = p + 1 

    return xg 

NumPy에서 쉽게 벡터화 할 수 있습니까?

return xg/n 

: 다음은 n으로 나누기가 매우 끝으로 이동 될 수 있습니다 것입니다

xg[p]=(i * xv[0] + j * xv[1] + k * xv[2] + l * xv[3])/n 

:

Answer 1

당신이 할 수있는 첫번째 간단한 일

하나에 세 가지 계산을 설정하기 위해 사용하는 방송입니다

그런 다음 4 개의 곱셈을 분리하여 결과를 별도로 저장 한 다음 끝에 결합 할 수 있습니다. 이제 우리는이 :

는

xg = np.empty((nsize,4)) # n.b. zeros not required 
p = 0 

for i in range (0, n + 1): 
    for j in range (0, n + 1 - i): 
     for k in range (0, n + 1 - i - j): 
      l = n - i - j - k 
      xg[p,0] = i 
      xg[p,1] = j 
      xg[p,2] = k 
      xg[p,3] = l 
      p = p + 1 

return (xg[:,:,None] * xv).sum(1)/n 

하단에 xg[:,:,None]의 비결은 (nsize,n) * (n,3)을 방송하는 것입니다 - 우리가 xv의 열을 곱되고있는에 의존하지 않는 (nsize,n,3)i,j,k,l 때문에에 (nsize,n)xg를 확장합니다.

우리는 루프에서 컴퓨팅 l를 생략하고 대신 한 번 바로 전에 return에서 모든 작업을 수행 할 수 있습니다 이제

xg[:,3] = n - xg[:,0:3].sum(1) 

당신이에 따라 벡터화 된 방법으로 i,j,k를 만드는 방법을 알아 내야한다 할 필요 p.

일반적으로 가장 안쪽 루프의 코드를보고 최대한 많은 루프를 벗어나 이러한 문제를 "내부에서"처리하는 것이 가장 쉽습니다. 이 작업을 반복하면 결국 루프가 생기지 않게됩니다.

Answer 2

종속 중첩 루프는 없지만 메모리가 낭비됩니다 (벡터화 문제에서 평소보다 많이 발생 함). for 루프에서 3d 상자의 작은 하위 집합을 처리하고 있습니다. (n+1)(n+2)(n+3)/6 개만 사용하여 (n+1)^3 개의 항목을 생성해도 상관 없으며 메모리에 저장하면 벡터화 된 버전이 훨씬 더 빠를 것입니다.

내 제안, 아래 설명 :

import numpy as np 

def tet_vect(n): 

    xv = np.array([ 
     [-1.,-1.,-1.], 
     [ 1.,-1.,-1.], 
     [-1., 1.,-1.], 
     [-1.,-1., 1.], 
     ]) 

    # spanning arrays of a 3d grid according to range(0,n+1) 
    ii,jj,kk = np.ogrid[:n+1,:n+1,:n+1] 
    # indices of the triples which fall inside the original for loop 
    inds = (jj < n+1-ii) & (kk < n+1-ii-jj) 
    # the [i,j,k] indices of the points that fall inside the for loop, in the same order 
    combs = np.vstack(np.where(inds)).T 
    # combs is now an (nsize,3)-shaped array 

    # compute "l" column too 
    lcol = n - combs.sum(axis=1) 
    combs = np.hstack((combs,lcol[:,None])) 
    # combs is now an (nsize,4)-shaped array 

    # all we need to do now is to take the matrix product of combs and xv, divide by n in the end 
    xg = np.matmul(combs,xv)/n 

    return xg 

핵심 단계는 3 차원 그리드에 걸치는 ii,jj,kk 인덱스를 생성한다. np.ogrid은 브로드 캐스트 작업에서 전체 그리드를 참조하는 데 사용할 수있는 스패닝 배열을 만들기 때문에이 단계는 메모리 효율적입니다. 우리는 다음 단계에서 전체 배열 (n+1)^3 크기의 배열을 생성합니다 : inds 부울 배열은 관심있는 영역 내부에있는 3d 상자의 점을 선택합니다 (배열 브로드 캐스트를 사용하여 그 점을 선택합니다). 다음 단계에서 np.where(inds)은이 큰 배열의 truthy 요소를 선택하고 더 작은 수의 nsize 요소로 끝납니다. 단일 메모리 낭비 단계는 따라서 inds이 생성됩니다.

나머지는 간단합니다. 각 행에 [i,j,k] 색인이 포함 된 추가 열을 생성해야합니다.이 작업은 배열의 열을 합산하여 수행 할 수 있습니다 (다시 벡터화 연산).행에 각각 (i,j,k,l)을 포함하는 (nsize,4) 모양의 보조 배열을 얻은 다음이 객체의 행렬 곱하기 xv을 수행해야합니다.

크기가 작은 n 크기를 테스트하면 위의 함수가 사용자와 동일한 결과를 생성한다는 것을 알 수 있습니다. n = 100에 대한 타이밍 : 원본 1.15 초, 벡터화 된 19ms.