[개발자를 위한 실전 선형대수학] Chapter 6 연습문제

[chapter6] 행렬 응용 : 데이터 분석에서의 행렬

6-1

공분산 행렬을 상관 행렬로 변환하기

변환하는 과정에서 각 행렬 원소(변수 쌍 간의 공분산)를 두 변수의 분산들의 곱으로 나눈다.

이것은 공분산 행렬에 각 변수의 역 표준 편차(분산의 제곱근)를 가진 대각행렬을 앞과 뒤에서 곱해서 구현한다.

행렬을 곱하면서 분산으로 나누어야 하기 때문에 표준편차를 반전한 것이다.

잘 이해가 안가서 해석하자면..

공분산 행렬을 상관 행렬로 변환하려면 표준편차로 나누어야 한다.

표준편차는 분산(대각선 원소)의 제곱근이다.

공분산 행렬의 원소들 모두를 표준편차로 나눠야하기 때문에, 왼쪽 오른쪽에 각각 역 표준편차를 곱해준다.

= 행렬의 왼쪽에서 곱하면 행 전체에 작용하고, 행렬의 오른쪽에서 곱하면 열 전체에 작용하기 때문.

$$ R = SCS $$

이 문제에서는 식 6-1을 파이썬 코드로 변환해서 공분산 행렬로부터 상관 행렬을 계산하는 것이다.

식 6-1

# read the data into a pandas dataframe
url  = 'https://archive.ics.uci.edu/ml/machine-learning-databases/communities/communities.data'
data = pd.read_csv(url,sep=',',header=None)

# attach column labels (don't worry, I didn't type this all in by hand, lol)
data.columns = [ 'state', 'county', 'community', 'communityname', 'fold', 'population', 'householdsize', 'racepctblack', 'racePctWhite',
'racePctAsian', 'racePctHisp', 'agePct12t21', 'agePct12t29', 'agePct16t24', 'agePct65up', 'numbUrban', 'pctUrban', 'medIncome', 'pctWWage',
'pctWFarmSelf', 'pctWInvInc', 'pctWSocSec', 'pctWPubAsst', 'pctWRetire', 'medFamInc', 'perCapInc', 'whitePerCap', 'blackPerCap', 'indianPerCap',
'AsianPerCap', 'OtherPerCap', 'HispPerCap', 'NumUnderPov', 'PctPopUnderPov', 'PctLess9thGrade', 'PctNotHSGrad', 'PctBSorMore', 'PctUnemployed', 'PctEmploy',
'PctEmplManu', 'PctEmplProfServ', 'PctOccupManu', 'PctOccupMgmtProf', 'MalePctDivorce', 'MalePctNevMarr', 'FemalePctDiv', 'TotalPctDiv', 'PersPerFam', 'PctFam2Par',
'PctKids2Par', 'PctYoungKids2Par', 'PctTeen2Par', 'PctWorkMomYoungKids', 'PctWorkMom', 'NumIlleg', 'PctIlleg', 'NumImmig', 'PctImmigRecent', 'PctImmigRec5',
'PctImmigRec8', 'PctImmigRec10', 'PctRecentImmig', 'PctRecImmig5', 'PctRecImmig8', 'PctRecImmig10', 'PctSpeakEnglOnly', 'PctNotSpeakEnglWell', 'PctLargHouseFam', 'PctLargHouseOccup',
'PersPerOccupHous', 'PersPerOwnOccHous', 'PersPerRentOccHous', 'PctPersOwnOccup', 'PctPersDenseHous', 'PctHousLess3BR', 'MedNumBR', 'HousVacant', 'PctHousOccup', 'PctHousOwnOcc',
'PctVacantBoarded', 'PctVacMore6Mos', 'MedYrHousBuilt', 'PctHousNoPhone', 'PctWOFullPlumb', 'OwnOccLowQuart', 'OwnOccMedVal', 'OwnOccHiQuart', 'RentLowQ', 'RentMedian',
'RentHighQ', 'MedRent', 'MedRentPctHousInc', 'MedOwnCostPctInc', 'MedOwnCostPctIncNoMtg', 'NumInShelters', 'NumStreet', 'PctForeignBorn', 'PctBornSameState', 'PctSameHouse85',
'PctSameCity85', 'PctSameState85', 'LemasSwornFT', 'LemasSwFTPerPop', 'LemasSwFTFieldOps', 'LemasSwFTFieldPerPop', 'LemasTotalReq', 'LemasTotReqPerPop', 'PolicReqPerOffic', 'PolicPerPop',
'RacialMatchCommPol', 'PctPolicWhite', 'PctPolicBlack', 'PctPolicHisp', 'PctPolicAsian', 'PctPolicMinor', 'OfficAssgnDrugUnits', 'NumKindsDrugsSeiz', 'PolicAveOTWorked', 'LandArea',
'PopDens', 'PctUsePubTrans', 'PolicCars', 'PolicOperBudg', 'LemasPctPolicOnPatr', 'LemasGangUnitDeploy', 'LemasPctOfficDrugUn', 'PolicBudgPerPop', 'ViolentCrimesPerPop',
 ]

# have a look at the data
data

# 숫자 데이터만 추출
numberDataset = data._get_numeric_data()

# 필요없는 열 제거 + Numpy 배열로 변환
dataMat = numberDataset.drop(['state','fold'],axis=1).values
dataMat

# 각 변수의 평균 계산
datamean = np.mean(dataMat,axis=0)

# 브로드캐스팅 사용해서 평균 제거
# 각 열에서 자기 평균을 뺌. 즉, 각 변수의 평균이 0인 데이터
dataMatM = dataMat - datamean

# 평균이 진짜 0인지 확인
# 거의 0이면 정
print(np.mean(dataMatM[:,0]))


# 공분산 행렬 계산
covMat = dataMatM.T @ dataMatM  
covMat /= (dataMatM.shape[0]-1) 

# 색상 범위 조절
clim = np.max(np.abs(covMat)) * .2

# 이미지로 표현
plt.figure(figsize=(6,6))
plt.imshow(covMat,vmin=-clim,vmax=clim,cmap='gray')
plt.colorbar()
plt.title('Data covariance matrix')
plt.savefig('Figure_06_01.png',dpi=300)
plt.show()

 

실행하면 아래와 같이 나온다.

답 코드(공분산 행렬 > 상관 행렬)

variances = np.diag(covMat) # 대각원소만 추출 = 분산
standard_devs = np.sqrt( variances ) # 분산에 제곱근 = 표준편차
S = np.diag( 1/standard_devs ) # 역 표준편차 대각행렬

# 상관행렬 계산
# 왼쪽 S = 행 기준 정규화, 오른쪽 S = 열 기준 정규화
corrMat = S @ covMat @ S


# and show the matrices
fig,axs = plt.subplots(1,2,figsize=(13,6))
h1 = axs[0].imshow(covMat,vmin=-clim,vmax=clim,cmap='gray')
axs[0].set_title('Data covariance matrix',fontweight='bold')

h2 = axs[1].imshow(corrMat,vmin=-.5,vmax=.5,cmap='gray')
axs[1].set_title('Data correlation matrix',fontweight='bold')

fig.colorbar(h1,ax=axs[0],fraction=.045)
fig.colorbar(h2,ax=axs[1],fraction=.045)

plt.tight_layout()
plt.savefig('Figure_06_06.png',dpi=300)
plt.show()

실행하면 아래처럼 나온다.

그래프가 좀 더 선명해졌다.

= 상관행렬은 각 변수를 표준편차로 정규화해서 스케일 차이를 제거하니까, 값의 크기(분산)에 묻히던 관계 패턴만 남아서 대비가 또렷해지는 것이다.
공분산은 “크기+관계”를 같이 보여주고, 상관행렬은 “관계만” 보여주기 때문에, 시각적으로는 구조가 더 균일하고 또렷하게 느껴진다.

6-2

NumPy의 np.corrcoef()로 데이터 행렬을 입력받아 상관 행렬을 반환하기

# dataMat의 구조는 (샘플 수, 변수 수)
# 즉 행 = 한사람 / 열 = 변수
# transpose하지 않고 사용하게 되면 '사람들 간의 상관관계'를 구하게 된다.
# (np.corrcoef는 행을 변수라고 생각하기 때문)
covMat = np.corrcoef(dataMat.T)

6-3

아래 변환행렬을 사용해서 그래프 만들기

$$ T = \begin{bmatrix}1 & .5 \\0 & .5 \\\end{bmatrix} $$

정답 코드

T = np.array([
    [1, .5],
    [0, .5]
])

# np.linspace(start, end, num) = start부터 end까지 포함해서 num개의 숫자를 균등 간격으로 생성
# 원을 20드분한 각도를 하나씩 만들되, 시작점(0)과 끝점(2pi)을 중복하지 않게 함
theta = np.linspace(0, 2*np.pi - 2*np.pi/20, 20)
origPoints = np.vstack((np.cos(theta), np.sin(theta))) # 단위 원 위의 좌표들

# 변형
transformedPoints = T @ origPoints

plt.figure(figsize=(6,6))
plt.plot(origPoints[0,:],origPoints[1,:],'ko',label='Original')
plt.plot(transformedPoints[0,:],transformedPoints[1,:],'s',
         color=[.7,.7,.7],label='Transformed')

plt.axis('square')
plt.xlim([-2,2])
plt.ylim([-2,2])
plt.legend()
plt.savefig('Figure_06_08.png',dpi=300)
plt.show()

실행하면 아래처럼 나온다.

6-4

아래 변환행렬을 사용해 동영상 만들기

$$ T = \begin{bmatrix}(1 - \phi /3) & 0 \\0 & \phi\\\end{bmatrix} $$

정답 코드

def aframe(ph):
  T = np.array([
      [1-ph/3, 0],
      [0, ph]
  ])

  # cos, sin 곡선에 같은 선형변환을 적용
  P1 = T@Y1 
  P2 = T@Y2

  plth1.set_xdata(P1[0,:])
  plth1.set_ydata(P1[1,:])

  plth2.set_xdata(P2[0,:])
  plth2.set_ydata(P2[1,:])

  # export the plot handles
  return (plth1,plth2)


# 위상 theta에 따른 cos, sin 파형
th = np.linspace(0,2*np.pi,100) # x좌표 theta
Y1 = np.vstack((th,np.cos(th))) # y좌표1 cos(theta)
Y2 = np.vstack((th,np.sin(th))) # y좌표2 sin(theta)


# setup figure
fig,ax = plt.subplots(1,figsize=(12,6))

plth1, = ax.plot(Y1[0,:],Y1[1,:],'ko') # cos 곡선
plth2, = ax.plot(Y2[0,:],Y2[1,:],'s',color=[.7,.7,.7]) # sin 곡선
ax.set_ylim([-2,2])


# define phases and run animation
phi = 1-np.linspace(-1,1-1/40,40)**2 # ∩ 모양 곡선 (-1부터 1까지 40개, 제곱은 가운데가 최대)
animation.FuncAnimation(fig, aframe, phi, interval=50, repeat=True)

6-5

3차원 욕조 사진을 매끄럽게 만들기.

convolve2d 함수 출력의 데이터 유형은 float64이다. 그러나 plt.imshow는 사진이 제대로 그려지지 않을 것이기 때문에 결과를 unit8로 변환해야 한다.

정답 코드

# (H, W, 3) 형태의 컬러 이미지
# H = 높이, W = 너비, 3 = RGB 채널
bathtub = io.imread('bathtub.jpg')

smooth_bathtub = np.zeros(bathtub.shape)


# 컨볼루션 커널(가우시안 필터) 만들기
# 커널 크기 : 29x29
# 값이 클 수록 더 강한 블러
kernelN = 29 

# 좌표 격자 생성
# X,Y는 각각 (29, 29) 배열. 커널 내부의 좌표계
Y,X     = np.meshgrid(np.linspace(-3,3,kernelN),np.linspace(-3,3,kernelN))

# 가우시안 커널 정의 : 가까운 픽셀은 크게, 먼 픽셀은 작게 가중
kernel  = np.exp( -(X**2+Y**2)/20 )

# 커널 정규화
# 커널 전체 합 = 1, 정규화하지 않으면 이미지가 전체적으로 밝아지거나 어두워짐
kernel  = kernel / np.sum(kernel)

# RGB 각 채널을 따로 처리
# convolve2d는 2D 배열만 처리하기 때문
for i in range(smooth_bathtub.shape[2]):
  smooth_bathtub[:,:,i] = convolve2d(bathtub[:,:,i], kernel, mode="same")

fig = plt.figure(figsize=(10,6))
plt.imshow(smooth_bathtub.astype(np.uint8))
plt.show()

6-6

각 채널의 가우스 폭 파라미터를 변경한다.

R채널 폭 0.5 / 패널 폭 5 / B 채널 폭 50

정답 코드

# (H, W, 3) 형태의 컬러 이미지
# H = 높이, W = 너비, 3 = RGB 채널
bathtub = io.imread('bathtub.jpg')

smooth_bathtub = np.zeros(bathtub.shape)

# 홀수 크기 = 중심 픽셀 명확(15 + 1+ 15)
# 큰 폭(50)을 담기 위해 커널 크기도 커짐
kernelN = 31
kernelWidths = [.5, 5, 50] # 값이 클수록 → 더 넓고 완만한 커널 → 더 강한 블러

_, axs = plt.subplots(1,3,figsize=(12,6))

for i in range(smooth_bathtub.shape[2]):
  # 커널 좌표계
  Y,X     = np.meshgrid(np.linspace(-3,3,kernelN),np.linspace(-3,3,kernelN))
  # 가우시안 커널 생성
  kernel  = np.exp( -(X**2+Y**2) / kernelWidths[i] ) # kernelWidth로 블러 강도 조절
  # 커널 정규화
  kernel  = kernel / np.sum(kernel)

  axs[i].imshow(kernel, cmap='gray')
  axs[i].set_title(f'Width:{kernelWidths[i]} ({"RGB"[i]} channel)')

  smooth_bathtub[:,:,i] = convolve2d(bathtub[:,:,i], kernel, mode='same')

plt.savefig('Figure_06_10.png',dpi=300)
plt.show() # close the kernels figure


# show the smoothed image
fig = plt.figure(figsize=(10,6))
plt.imshow(smooth_bathtub.astype(np.uint8))
plt.show()

출력하면 이렇게 나온다.

6-7

필터 커널을 변경해서 수평선과 수직선을 식별하는 다른 이미지 특징 감지

# 세로 에지 커널
VK = np.array([
    [1,0,-1],
    [1,0,-1],
    [1,0,-1]
])

# 가로 에지 커널
HK = np.array([
    [1,1,1],
    [0,0,0],
    [-1,-1,-1]
])

fig,ax = plt.subplots(2,2,figsize=(16,8))

ax[0,0].imshow(VK,cmap='gray')
ax[0,0].set_title('Vertical kernel')
ax[0,0].set_yticks(range(3))

ax[0,1].imshow(HK,cmap='gray')
ax[0,1].set_title('Horizontal kernel')
ax[0,1].set_yticks(range(3))

bathtub2d = color.rgb2gray(bathtub)

# VK와 내적을 구한다. 좌우 밝기 차이가 크면 값이 크다.
convres = convolve2d(bathtub2d,VK,mode='same')
ax[1,0].imshow(convres,cmap='gray',vmin=0,vmax=.01)
ax[1,0].axis('off')

# HK와 내적을 구한다. 위아래 밝기 차이를 검출한다.
convres = convolve2d(bathtub2d,HK,mode='same')
ax[1,1].imshow(convres,cmap='gray',vmin=0,vmax=.01)
ax[1,1].axis('off')

# 즉 주변 픽셀 벡터와 커널의 내적을 계산해서
# 내적의 절댓값이 크다 = 그 방향으로 밝기 변화가 크다 = 경계다.

plt.savefig('Figure_06_11.png',dpi=300)
plt.show()

출력하면 아래처럼 나온다.