이번 포스팅에서는 값 배열에 가중치 배열을 곱해서 합한 가중합(weighted sum)을 구하는 3가지 방법을 소개하겠습니다. 


a 를 가중치, b 를 값 배열이라고 했을 때, 


(1) 내적을 이용한 가중합 계산: np.dot(a, b) or np.matmul(a, b)

(2) 브로드캐스팅(broadcasting)을 이용하여 가중치와 값을 원소끼리 곱한 후 합하는

     np.sum(a.reshape(5, 1) * b, axis=0)

(3) repeat()로 가중치를 값 배열 1축만큼 반복 생성한 후, 가중치와 값의 원소끼리 곱한 후 합하는, 

     np.sum(a.reshape(5, 1).repeat(3, axis=1) * b, axis=0)




먼저, numpy를 import하고, 예제로 사용할 가중치 배열 a와, 값의 행렬 b를 만들어보겠습니다. 



import numpy as np


# weights

a = np.array([0.5, 0.3, 0.1, 0.08, 0.02])


print('a shape:', a.shape)

a shape: (5,)


print(a)

[0.5  0.3  0.1  0.08 0.02]



# values

b = np.arange(15).reshape(5, 3)


print('b shape:', b.shape)

b shape: (5, 3)


print(b)

[[ 0  1  2]

 [ 3  4  5]

 [ 6  7  8]

 [ 9 10 11]

 [12 13 14]]

 




  (1) 내적을 이용한 가중합 계산: np.dot(a, b) 또는 np.matmul(a, b)


가장 편리한 방법은 np.dot() 또는 np.matmul() 메소드를 사용하여 내적(inner prodct, dot product)을 계산하는 것입니다. 이때 가중치 벡터 a 에 대해서는 형태 변환(reshape)을 할 필요가 없이 그대로 사용할 수 있습니다.  



np.dot(a, b)

Out[2]: array([2.46, 3.46, 4.46])


np.matmul(a, b)

Out[3]: array([2.46, 3.46, 4.46])

 




  (2) Broadcasting을 이용하여 가중치와 값을 원소끼리 곱한 후, axis=0으로 합하기


이번에는 위의 (1) 내적을 계산의 각 단계별로 분리해서 순서대로 해보겠습니다. 가중치 a와 값 b의 원소끼리 곱한 후에, axis=0을 기준으로 합할 것입니다. 


먼저, 가중치 a와 값 b를 원소끼리 곱하기 위해 가중치 a의 형태(shape)를 기존의 (5,)에서 a.reshape(5, 1) 을 적용하여 (5, 1) 의 형태로 변환을 해줍니다. 값이 들어있는 배열 b의 형태는 (5, 3) 이므로 가중치 배열 a의 (5, 1) 형태를 값 배열 b에 곱해주면 ==> 서로 형태가 같지 않으므로 numpy 는 가중치 a 배열 (5, 1) 을 (5, 3)으로 자동으로 형태 변환을 시켜서 값 배열 b 의 (5, 3) 형태와 동일하게 맞추어 주어 원소간 곱을 해줍니다. 이러한 기능을 브로드캐스팅(boradcasting) 이라고 합니다. 



# shape of a_rs and b are different

a_rs = a.reshape(5, 1)

print(a_rs.shape)

print(a_rs)

(5, 1)


print(b.shape)

(5, 3)


# multiply using boradcasting of a_rs

a_rs_b_mult = a_rs * b


print(a_rs_b_mult.shape)

(5, 3)


print(a_rs_b_mult)

[[0.   0.5  1.  ]

 [0.9  1.2  1.5 ]

 [0.6  0.7  0.8 ]

 [0.72 0.8  0.88]

 [0.24 0.26 0.28]]



# weighted sum

np.sum(a_rs_b_mult, axis=0)

Out[9]: array([2.46, 3.46, 4.46])



* numpy 배열들의 다른 차원의 배열 간 산술연산 시 Broadcasting 은 아래 포스팅을 참고하세요. 

https://rfriend.tistory.com/287




  (3) repeat()로 가중치를 반복 생성한 후, 가중치와 값을 원소끼리 곱한 후 합하기


위의 (2)번에서는 가중치 배열 a의 형태를 바꾼 후의 a_rs 배열과 값 b 배열을 곱할 때, 사람 눈에는 보이지않게 numpy가 알아서 자동으로 가중치 a_rs 배열 (5, 1) 형태를 브로드캐스팅(broadcasting)을 해주어서 (5, 3) 형태로 만들어서 원소끼리 곱해주었습니다. 




반면에, 이번 (3)번에서는 사람이 repeat(n, axis) 메소드를 사용해서 명시적으로 배열을 n번 만큼 axis 축을 기준으로 반복해주어서 (2)번의 브로드캐스팅의 역할을 수행해주는 것입니다. 


구현 관점에서 보면 브로드케스팅이 편리한 장점이 있고, 반면에 repeat() 메소드로 명시적으로 기입을 해주면 코딩하는 사람이 이해하기 쉬운 장점이 있습니다. 



# match the shape of a and b by repeatition 

a_rs_rp = a.reshape(5, 1).repeat(3, axis=1)


print(a_rs_rp.shape)

(5, 3)


print(a_rs_rp)

[[0.5  0.5  0.5 ]

 [0.3  0.3  0.3 ]

 [0.1  0.1  0.1 ]

 [0.08 0.08 0.08]

 [0.02 0.02 0.02]]



# multiplication of a_rs_rp and b per each elements

a_rs_rp_b_mult = a_rs_rp * b


print(a_rs_rp_b_mult.shape)

(5, 3)


print(a_rs_rp_b_mult)

[[0.   0.5  1.  ]

 [0.9  1.2  1.5 ]

 [0.6  0.7  0.8 ]

 [0.72 0.8  0.88]

 [0.24 0.26 0.28]]



# weighted sum

np.sum(a_rs_rp_b_mult, axis=0)

Out[17]: array([2.46, 3.46, 4.46])

 



많은 도움이 되었기를 바랍니다. 

이번 포스팅이 도움이 되었다면 아래의 '공감~'를 꾹 눌러주세요. 



728x90
반응형
Posted by Rfriend
,

학교 다닐 때 행렬로 연립방정식 풀었던 기억이 날 듯 합니다. 선형대수(Linear Algebra)는 통계, 기계학습, 공학, 영상/이미지 처리 등 여러 분야에서 활용이 됩니다. 선형대수를 전부 다루려면 너무나 방대하므로, 이번 포스팅에서는 Python의 NumPy에 있는 선형대수(Linear Algebra) 함수들 중에서 자주 사용하는 함수에 대해서만 선별적으로 소개하겠습니다. 그리고 선형대수의 이론적인 부분은 별도로 참고할 수 있는 링크를 달도록 하겠습니다. 


  • 단위행렬 (Unit matrix): np.eye(n)
  • 대각행렬 (Diagonal matrix): np.diag(x)
  • 내적 (Dot product, Inner product): np.dot(a, b)
  • 대각합 (Trace): np.trace(x)
  • 행렬식 (Matrix Determinant): np.linalg.det(x)
  • 역행렬 (Inverse of a matrix): np.linalg.inv(x)
  • 고유값 (Eigenvalue), 고유벡터 (Eigenvector): w, v = np.linalg.eig(x)
  • 특이값 분해 (Singular Value Decomposition): u, s, vh = np.linalg.svd(A)
  • 연립방정식 해 풀기 (Solve a linear matrix equation): np.linalg.solve(a, b)
  • 최소자승 해 풀기 (Compute the Least-squares solution): m, c = np.linalg.lstsq(A, y, rcond=None)[0]




 1. 단위행렬 혹은 항등행렬 (Unit matrix, Identity matrix): np.eye(n)


단위행렬은 대각원소가 1이고, 나머지는 모두 0인 n차 정방행렬을 말하며, numpy의 eye() 함수를 사용해서 만들 수 있습니다. 


* 참고 링크 : https://rfriend.tistory.com/141



import numpy as np


unit_mat_4 = np.eye(4)


print(unit_mat_4)

[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]

 




  1. 대각행렬 (Diagonal matrix): np.diag(x)


대각행렬은 대각성분 이외의 모든 성분이 모두 '0'인 n차 정방행렬을 말합니다. 아래 예시의 행렬에서 빨간색으로 표시한 원소를 '0'으로 바꾼 행렬이 대각행렬입니다. 


* 참고 링크 : http://rfriend.tistory.com/141


 

In [1]: import numpy as np


In [2]: x = np.arange(9).reshape(3, 3)


In [3]: print(x)

[[0 1 2]

 [3 4 5]

 [6 7 8]]


In [4]: np.diag(x)

Out[4]: array([0, 4, 8])


In [5]: np.diag(np.diag(x))

Out[5]:

array([[0, 0, 0],

        [0, 4, 0],

        [0, 0, 8]])





  2. 내적 (Dot product, Inner product): np.dot(a, b), a.dot(b)


matrix dot product, inner product, scalar product, projection product

Python에서 '*' 를 사용한 두 행렬 간 곱은 원소 간 곱(element-wise product)을 반환하며, 선형대수에서 말하는 행렬 간 내적 곱을 위해서는 np.dot() 함수를 이용해야 합니다. 


원소 간 곱 (element-wise product)

: a*b

내적 (dot product, inner product)

: np.dot(a, b)

 

In [6]: a = np.arange(4).reshape(2, 2)


In [7]: print(a)

[[0 1]

 [2 3]]


In [8]: a*a

Out[8]:

array([[0, 1],

        [4, 9]])


In [6]: a = np.arange(4).reshape(2, 2)


In [7]: print(a)

[[0 1]

 [2 3]]


In [9]: np.dot(a, a)

Out[9]:

array([[ 2, 3],

        [ 6, 11]])



np.dot(a, b) NumPy 함수와 a.dot(b)의 배열 메소드의 결과는 동일합니다. 



In [10]: a.dot(a)

Out[10]:

array([[ 2, 3],

        [ 6, 11]])

 




  3. 대각합 (Trace): np.trace(x)


정방행렬의 대각에 위치한 원소를 전부 더해줍니다. 

아래의 2차 정방행렬 예의 대각합은 0+5+10+15 = 30 이 됩니다. (파란색으로 표시함)



In [12]: b = np.arange(16).reshape(4, 4)


In [13]: print(b)

[[ 0 1 2 3]

 [ 4 5 6 7]

 [ 8 9 10 11]

 [12 13 14 15]]


In [14]: np.trace(b)

Out[14]: 30

 



3차원 행렬에 대해서도 대각합을 구할 수 있습니다. 2차원은 대각선 부분의 원소 값을 전부 더하면 되지만 3차원 행렬에서는 대각(diagonal)이 어떻게 되나 좀 헷갈릴 수 있겠습니다. 아래의 3차원 행렬의 대각합을 구하는 예를 살펴보면, [0+12+24, 1+13+25, 2+14+26] = [36, 39, 42] 가 됩니다. 



In [15]: c = np.arange(27).reshape(3, 3, 3)


In [16]: print(c)

[[[ 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 25 26]]]


In [17]: np.trace(c)

Out[17]: array([36, 39, 42])

 




  4. 행렬식 (Matrix Determinant): np.linalg.det(x)


역행렬이 존재하는지 여부를 확인하는 방법으로 행렬식(determinant, 줄여서 det)이라는 지표를 사용합니다. 이 행렬식이 '0'이 아니면 역행렬이 존재하고, 이 행렬식이 '0'이면 역행렬이 존재하지 않습니다. 


* 참고 링크 : http://rfriend.tistory.com/142


아래의 예에서 array([[1, 2], [3, 4]]) 의 행렬식이 '-2.0'으로서, '0'이 아니므로 역행렬이 존재한다고 판단할 수 있습니다. 



In [18]: d = np.array([[1, 2], [3, 4]])


In [19]: np.linalg.det(a)

Out[19]: -2.0

 




  5. 역행렬 (Inverse of a matrix): np.linalg.inv(x)


역행렬은 n차정방행렬 Amn과의 곱이 항등행렬 또는 단위행렬 In이 되는 n차정방행렬을 말합니다. A*B 와 B*A 모두 순서에 상관없이 곱했을 때 단위행렬이 나오는 n차정방행렬이 있다면 역행렬이 존재하는 것입니다.

역행렬은 가우스 소거법(Gauss-Jordan elimination method), 혹은 여인수(cofactor method)로 풀 수 있습니다. 




In [20]: a = np.array(range(4)).reshape(2, 2)


In [21]: print(a)

[[0 1]

 [2 3]]


In [22]: a_inv = np.linalg.inv(a)


In [23]: a_inv

Out[23]:

array([[-1.5, 0.5],

        [ 1. , 0. ]])




위의 예제에서 np.linalg.inv() 함수를 사용하여 푼 역행렬이 제대로 푼 것인지 확인을 해보겠습니다. 역행렬의 정의에 따라서 원래의 행렬에 역행렬을 곱하면, 즉, a.dot(a_inv) 또는 np.dot(a, a_inv) 를 하면 단위행렬(unit matrix)가 되는지 확인해보겠습니다. 



In [24]: a.dot(a_inv)

Out[24]:

array([[1., 0.],

         [0., 1.]])

 




  6. 고유값 (Eigenvalue), 고유벡터 (Eigenvector): w, v = np.linalg.eig(x)


정방행렬 A에 대하여 Ax = λx  (상수 λ) 가 성립하는 0이 아닌 벡터 x가 존재할 때 상수 λ 를 행렬 A의 고유값 (eigenvalue), x 를 이에 대응하는 고유벡터 (eigenvector) 라고 합니다. 


np.linalg.eig() 함수는 고유값(eigenvalue) w, 고유벡터(eigenvector) v 의 두 개의 객체를 반환합니다. 


In [25]: e = np.array([[4, 2],[3, 5]])


In [26]: print(e)

[[4 2]

[3 5]]


In [27]: w, v = np.linalg.eig(e)


#  w: the eigenvalues lambda

In [28]: print(w)

[2. 7.]


# v: the corresponding eigenvectors, one eigenvector per column

In [29]: print(v)

[[-0.70710678 -0.5547002 ]

[ 0.70710678 -0.83205029]]

 



고유벡터는 배열 인덱싱하는 방법을 사용해서 각 고유값에 대응하는 고유벡터를 선택할 수 있습니다. 



# eigenvector of eigenvalue lambda 2

In [30]: print(v[:, 0]

[-0.70710678 0.70710678]


# eigenvector of eigenvalue labmda 7

In [31]: print(v[:, 1]

[-0.5547002 -0.83205029]

 




  7. 특이값 분해 (Singular Value Decomposition): u, s, vh = np.linalg.svd(A)


특이값 분해는 고유값 분해(eigen decomposition)처럼 행렬을 대각화하는 한 방법으로서, 정방행렬뿐만 아니라 모든 m x n 행렬에 대해 적용 가능합니다. 특이값 분해는 차원축소, 데이터 압축 등에 사용할 수 있습니다. 이론적인 부분은 설명하자면 너무 길기 때문에 이 포스팅에서는 설명하지 않겠으며, 아래의 링크를 참고하시기 바랍니다. 


* 참고 링크 : http://rfriend.tistory.com/185


아래의 np.linalg.svd(A) 예제는 위의 참고 링크에서 사용했던 예제와 동일한 것을 사용하였습니다. 


In [32]: A = np.array([[3,6], [2,3], [0,0], [0,0]])


In [33]: print(A)

[[3 6]

 [2 3]

 [0 0]

 [0 0]]


In [34]: u, s, vh = np.linalg.svd(A)


In [35]: print(u)

[[-0.8816746 -0.47185793 0. 0. ]

 [-0.47185793 0.8816746 0. 0. ]

 [ 0. 0. 1. 0. ]

 [ 0. 0. 0. 1. ]]


In [36]: print(s)

[7.60555128 0.39444872]


In [37]: print(vh)

[[-0.47185793 -0.8816746 ]

 [ 0.8816746 -0.47185793]]

 




  8. 연립방정식 해 풀기 (Solve a linear matrix equation): np.linalg.solve(a, b)


아래의 두 개 연립방정식의 해(x0, x1)를 np.linalg.solve(a, b) 함수를 사용하여 풀어보겠습니다. 



위의 연립방정식을 어떻게 행렬로 입력하고, np.linalg.solve(a, b)에 입력하는지 유심히 살펴보시기 바랍니다. 



In [38]: a = np.array([[4, 3], [3, 2]])


In [39]: b = np.array([23, 16])


In [40]: x = np.linalg.solve(a, b)


In [41]: print(x)

[2. 5.]

 



NumPy 가 제대로 x0, x1의 해를 풀었는지 확인해보겠습니다. x0=2, x1=5 가 해 맞네요!



In [42]: np.allclose(np.dot(a, x), b)

Out[42]: True

 




  9. 최소자승 해 풀기 (Compute the Least-squares solution)

     : m, c = np.linalg.lstsq(A, y, rcond=None)[0]


회귀모형 적합할 때 최소자승법(Least-squares method)으로 잔차 제곱합을 최소화하는 회귀계수를 추정합니다. 


* 참고 링크 : https://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.lstsq.html



아래의 예에서는 회귀계수 m, y절편 c를 최소자승법을 사용해서 구해보겠습니다. 



In [43]: x = np.array([0, 1, 2, 3])


In [44]: y = np.array([-1, 0.2, 0.9, 2.1])


In [45]: A = np.vstack([x, np.ones(len(x))]).T


In [46]: A

Out[46]:

array([[0., 1.],

        [1., 1.],

        [2., 1.],

        [3., 1.]])


In [47]: m, c = np.linalg.lstsq(A, y, rcond=None)[0]


In [48]: print(m, c)

0.9999999999999999 -0.9499999999999997

 



아래의 그래프에서 점은 원래의 데이터이며, 빨간색 선은 최소자승법으로 추정한 회귀식의 적합선이 되겠습니다. 



In [49]: import matplotlib.pyplot as plt

    ...: plt.plot(x, y, 'o', label='Original data', markersize=10)

    ...: plt.plot(x, m*x + c, 'r', label='Fitted line')

    ...: plt.legend()

    ...: plt.show()




많은 도움이 되었기를 바랍니다. 


728x90
반응형
Posted by Rfriend
,