넘파이 알고 쓰자 - Linear Algebra Library 2 : 아인슈타인 표기법

안녕하세요. 지난 포스팅의 넘파이 알고 쓰자 - Linear Algebra Library 1 : Matrix&Vector Product에서는 넘파이의 선형대수 라이브러리에 구현되어 있는 다양한 행렬/벡터곱을 알아보고 dot 함수와 matmul 함수의 차이점에 대해서 분석해보았습니다. 오늘도 이어서 행렬/벡터곱과 관련된 내용이지만 이를 다른 표기법으로 작성하는 아인슈타인 표기법에 대해서 알아보도록 하겠습니다.

 

본격적으로 시작하기 전에 아인슈타인 표기법에 대해서 간단하게 알아보록 넘어가도록 하겠습니다. 아인슈타인 표기법의 정확한 정의는 아래의 링크를 참조해주시길 바랍니다. 링크

아인슈타인 표기법 - 위키백과, 우리 모두의 백과사전

몇 가지 중요한 부분을 보도록 하겠습니다. 예를 들어 아래의 합을 생각해보겠습니다.

$$y = c_{1}x + c_{2}x^{2} + \cdots c_{n}x^{n}$$

이 수식은 $\Sigma$ 기호를 이용하면 아래와 같이 간단하게 표현할 수 있습니다.

$$y = \sum_{i = 1}^{n} c_{i}x^{i}$$

여기서 아인슈타인 표기법은 이 $\Sigma$ 기호를 없앤상태로 아래와 같이 쓰게 됩니다.

$$y = c_{i}x^{i}$$

이때 중요한 점은 $x$위에 붙은 첨자 $i$는 지수승을 의미하지 않는 다는 것입니다. 아인슈타인 표기법에서는 이 중복된 첨자 $i$를 이용하여 계산을 하게 됩니다. 사실 첨자의 위치에 따라서 여벡터/벡터로 나누지만 이는 너무 복잡하기 때문에 간단하게 첨자가 붙은 것들은 그냥 전부 벡터라고 생각하셔도 됩니다. 즉, 위 수식은 $c$라는 벡터와 $x$라는 벡터를 동일한 첨자끼리 곱한다는 의미로 해석할 수 있다는 것이죠. 이 해석을 이용해서 몇 가지 변경점을 확인해보도록 하겠습니다.

 

1. 임의의 $1 \times n$ 행벡터 $v^{i}$와 $n \times 1$의 열벡터 $u_{i}$에 대한 내적 : $u_{i}v_{i} = u_{1}v_{1} + u_{2}v_{2} + \cdots u_{n}v_{n}$
2. 임의의 $m \times 1$ 열벡터 $u^{i}$와 $1 \times n$의 행벡터 $v_{j}$에 대한 외적 : $u^{i}v_{j} = A^{i}_{j}$. 이는 결과적으로 $m \times n$ 크기의 행렬 $A$를 얻는 다는 것입니다.
3. 임의의 $m \times n$ 행렬 $A^{i}_{j}$와 $n \times 1$ 열벡터 $v^{j}$가 주어졌을 때 두 곱의 결과를 $u^{i}$라고 했을 때 행렬과 벡터의 곱 : $u^{i} = A^{i}_{j}v^{j}$

이제 위의 내용을 바탕으로 넘파이에서는 어떤 식으로 구현되어 있는 지 알아보도록 하겠습니다. 

 

1. numpy.einsum(subscripts)

여기서 subscripts는 어떤 첨자로 연산을 취할 것인지 적어주시면 됩니다. 이전에 설명한 아인슈타인 표기법에서는 "첨자"를 기준으로 연산을 취한다는 것을 상기하시면 됩니다. 간단한 예제를 보도록 하겠습니다.

a = np.arange(25).reshape(5,5)
# [[ 0  1  2  3  4]
#  [ 5  6  7  8  9]
#  [10 11 12 13 14]
#  [15 16 17 18 19]
#  [20 21 22 23 24]]

np.einsum('ii', a) # 60

위의 예시의 경우에는 einsum 함수의 첨자로 "ii"를 입력하였습니다. a라는 행렬의 ii란 대각선을 의미합니다. 따라서, 이는 대각선 방향으로 더한다는 것이죠!! 즉, 행렬의 대각합(trace)와 동일한 결과입니다. 

a = np.arange(25).reshape(5,5)
# [[ 0  1  2  3  4]
#  [ 5  6  7  8  9]
#  [10 11 12 13 14]
#  [15 16 17 18 19]
#  [20 21 22 23 24]]

np.einsum('ii->i', a) # array([ 0,  6, 12, 18, 24])

위의 경우에는 einsum 함수의 첨자로 "ii -> i"를 입력하였습니다. a라는 행렬의 ii, 즉 대각성분을 취하여 i번째 대각성분을 가지는 벡터를 생성하라는 뜻입니다. 

a = np.arange(25).reshape(5,5)
# [[ 0  1  2  3  4]
#  [ 5  6  7  8  9]
#  [10 11 12 13 14]
#  [15 16 17 18 19]
#  [20 21 22 23 24]]

np.einsum('ij->i', a) # array([ 10,  35,  60,  85, 110])

이번에느 살짝 해석하기 어렵지만 그래도 한번 해보도록 하겠습니다. 위의 경우에는 einsum 함수의 첨자로 "ij -> i"를 입력하였습니다. 먼저, a라는 행렬의 ij번째 원소를 선택합니다. 간단하게 보기 위해서 i = 0이라고 가정하면 0j번째 원소이므로 00, 01, 02, 03, 04 번째 원소들을 선택한다는 것이죠. 해당 인덱스에 해당하는 값을 모두 더하여 새롭게 만드는 벡터의 0번째 요소로 하게 됩니다!! 즉, 10이 되는 것이죠. 이를 i = 0, 1, 2, 3, 4에 대해서 반복하기 때문에 행렬 내에서 axis=1 방향으로의 덧셈, 즉 np.sum(a, axis=1)과 동일한 결과를 가집니다.

c = np.arange(6).reshape(2,3)
# [[0 1 2]
#  [3 4 5]]

np.einsum('ji', c)
# array([[0, 3],
#        [1, 4],
#        [2, 5]])

np.einsum('ji->ij', c)
# array([[0, 3],
#        [1, 4],
#        [2, 5]])

이번에는 2가지를 보도록 하겠습니다. 두 함수 결과 모두 동일한 결과를 가지게 됩니다. 첫번째 함수를 해석해보면 행렬 c의 ji번째 요소를 선택한다는 뜻입니다. 즉, 행렬 c의 전치행렬의 의미하게 되는 것이죠. 다음 함수가 사실 더 직관적일 것입니다. 행렬 c의 ji번째 요소를 선택하여 새로 생성되는 행렬의 ij번째 요소로 추가한다는 뜻이기 때문에 전치행렬임을 쉽게 알 수 있습니다. 

b = np.arange(5)
# array([0, 1, 2, 3, 4])

np.einsum('i,i', b, b) # 30

이번 함수도 그리 어렵지는 않습니다. 그저 벡터 b가 주어졌을 때 i번째 요소를 제곱하여 덧셈을 취하는 것을 의미합니다. 그러므로 내적이죠. 

a = np.arange(25).reshape(5,5)
# [[ 0  1  2  3  4]
#  [ 5  6  7  8  9]
#  [10 11 12 13 14]
#  [15 16 17 18 19]
#  [20 21 22 23 24]]

b = np.arange(5)
# array([0, 1, 2, 3, 4])

np.einsum('ij,j', a, b) # array([ 30,  80, 130, 180, 230])

마지막으로 행렬/벡터 곱을 알아보도록 하겠습니다. a라는 행렬의 ij번째 요소와 b라는 벡터의 j번째 요소를 곱한 뒤 덧셈을 취한다는 뜻입니다. 

 

2. numpy.einsum_path(subscripts, optimize='greedy')

이번에도 아인슈타인 표기법을 이용한 연산을 소개하도록 하겠습니다. 이번 함수에서 다른 점은 이전 포스팅에서 설명한 matrix chain algorithm 처럼 연산횟수가 최소화되는 연산 순서를 찾아주는 것입니다. optimize는 greedy 알고리즘과 optimal 알고리즘 중 선택할 수 있습니다. 바로 아래의 예제를 보도록 하겠습니다.

 

a = np.random.rand(2, 2)
b = np.random.rand(2, 5)
c = np.random.rand(5, 2)
path_info = np.einsum_path('ij,jk,kl->il', a, b, c, optimize='greedy')
print(path_info[0]) # ['einsum_path', (1, 2), (0, 1)]
print(path_info[1])
# ['einsum_path', (1, 2), (0, 1)]
#   Complete contraction:  ij,jk,kl->il
#          Naive scaling:  4
#      Optimized scaling:  3
#       Naive FLOP count:  1.200e+02
#   Optimized FLOP count:  5.700e+01
#    Theoretical speedup:  2.105
#   Largest intermediate:  4.000e+00 elements
# --------------------------------------------------------------------------
# scaling                  current                                remaining
# --------------------------------------------------------------------------
#    3                   kl,jk->jl                                ij,jl->il
#    3                   jl,ij->il                                   il->il

shape이 서로다른 3개의 행렬을 만들었습니다. 각 행렬의 shape은 (2, 2), (2, 5), (5, 2)로 곱셈도 잘 정의되는 것을 볼 수 있습니다. 그 결과는 path_info에 저장되는데 0번 인덱스에는 비용을 최소화한 연산 순서를 알려주고 있습니다. 이 경우에는 a(bc)를 계산하는 것이 더 빠르다는 것을 알려주고 있습니다. 즉, 이 함수는 이전 포스팅의 np.linalg.multi_dot 함수와 유사한 결과입니다.