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

[chapter5] 행렬, 파트 2 : 행렬의 확장 개념

5-1

행렬 노름은 행렬이 가진 숫자의 크기와 관련이 있다. 이를 증명하기 위한 실험하기.

실험을 10번 반복하는데, 매번 10x10 난수 행렬을 만들고 프로베니우스 노름 계산하기.

그리고 이 실험을 40번 반복하고, 매번 행렬에 0과 50 사이의 다른 스칼라를 곱하기.

정답 코드

scalingVals = np.linspace(0,50,40) # 0부터 50까지를 40개 값으로 균등 분할한 배열
nExperiments = 10 # 같은 스케일에서 랜덤 행렬을 10번 생성

matrixNorms = np.zeros((len(scalingVals),nExperiments)) # (40, 10), 각 칸에 행렬의 노름 값 하나를 저장

# 스케일 값 순회
for si in range(len(scalingVals)):
  for expi in range(nExperiments):
    R = np.random.randn(10,10) * scalingVals[si] # 10×10 난수 행렬에 스칼라를 곱함
    matrixNorms[si,expi] = np.linalg.norm(R,'fro') # R의 프로베니우스 노름을 구해서 결과값에 넣음

plt.plot(scalingVals,np.mean(matrixNorms,axis=1),'ko-')
plt.xlabel('Matrix scalar')
plt.ylabel('Matrix Frobenius norm')
plt.savefig('Figure_05_07.png',dpi=300)
plt.show()

# check that norm=0 for zeros matrix
print(matrixNorms[0,:])

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

이 문제는 아래 두 가지를 보여준다.

• 스케일이 커지면 노름이 비례해서 커진다.
• 노름은 행렬의 구체적 형태와 상관없이 ‘전체 크기’만 측정한다

5-2

두 행렬 사이의 프로베니우스 거리를 1로 만드는 스칼라를 찾는 알고리즘 작성하기.

먼저 두 개의 행렬(같은 크기)를 입력으로 받아 두 행렬 사이의 프로베니우스 거리를 반환하는 파이썬 함수를 작성

그리고 두개의 NxN 난수 행렬을 만들고, 두 행렬에 스칼라곱셈을 하는 변수 s =1 을 만든다.

크기를 조정한 행렬들 사이에 프로베니우스 거리를 계산한다.

거리가 1 이상인 경우 스칼라를 0.9배로 설정하고 크기가 조정된 행렬들 사이의 거리를 다시 계산한다.(while문 안에서 수행)

프로베니우스 거리가 1 미만이 되면 while 루프는 종료하고 반복 횟수와 스칼라값을 기록한다.

정답 코드

# 프로베니우스 노름 구하기
def EuclideanDistance(M1,M2):
  D = M1-M2
  return np.sqrt(np.sum(D**2))

# 7X7 난수행렬 만들기
N = 7
A = np.random.randn(N,N)
B = np.random.randn(N,N)

numIters = 0 # 몇번 반복했는지
s = 1
while EuclideanDistance(s*A,s*B)>1:
  s *= .9
  numIters += 1 

print(f'Number of iterations: {numIters-1}')
print(f'Final value of scalar: {s/.9:.3f}')
print(f'Final Euclidean distance: {EuclideanDistance(s/.9*A,s/.9*B):.3f}')

행렬 거리(Frobenius 거리)는 스케일에 대해 정확히 비례해서 변한다.

5-3

대각합 연산과 유클리드 공식의 결과(프로베니우스 노름)가 동일하다는 것을 구현

M = 50
A = np.random.randn(M,M)

# 대각합
norm1 = np.sqrt(np.sum(np.diag(A.T@A)))

# 프로베니우스 노름
norm2 = np.sqrt(np.sum(A**2))

print(norm1-norm2) #0.0

5-4

행렬 이동이 행렬 노름에 미치는 영향 알아보기.

• 10x10 난수 행렬 만들고 프로베니우스 노름 계산
• for문 내부에 다음 단계 코딩
  • 노름의 일부만큼 행렬 이동
  • 원래 행렬에서 노름의 변화율 계산
  • 이동된 행렬과 원래 행렬 사이의 프로베니우스 거리 계산
  • 행렬과 원소사이의 상관관계 계산(np.flatten 사용해 벡터화된 행렬의 상관관계 계산
  • (이동시킬 값인 노름의 비율은 0~1 사이를 30개의 단계로 나누어 진행)

정답 코드

def EuclideanDistance(M1,M2):
  
  # matrix difference
  D = M1-M2

  # matrix distance
  return np.sqrt(np.sum(D**2))

shifting = np.linspace(0,1,30)

A = np.random.randn(10,10)
normA = np.linalg.norm(A, 'fro')

shiftingResults = np.zeros( (len(shifting),3) )
resultsNames = [ 'Change in norm (%)','Corr. with original','Frobenius distance' ]

for si in range(len(shifting)):
  # 노름의 일부만큼 행렬 이동
  As = A + shifting[si]*normA*np.eye(10)
  
  # 이동된 행렬과 원래 행렬 사이의 프로베니우스 거리 계산
  normShift = np.linalg.norm(As, 'fro')

  # 노름 변화
  # - 행렬 전체 크기가 얼마나 커졌는지 (%)
  shiftingResults[si,0] = 100 * (normShift - normA) / normA
  # 원래 행렬과의 상관계수
  # - 두 행렬의 방향이 얼마나 비슷한지(corrcoef = 피어슨 상관계수)
  shiftingResults[si,1] = np.corrcoef(A.flatten(), As.flatten())[0,1]
  # 프로베니우스 거리
  # - 얼마나 이동했는지
  shiftingResults[si,2] = EuclideanDistance(A, As)

## plotting!
_,axs = plt.subplots(1,3,figsize=(12,4))

for i in range(3):

  # plot the results
  axs[i].plot(shifting,shiftingResults[:,i],'ks-')
  axs[i].set_xlabel('Shifting (prop. of norm)')
  axs[i].set_ylabel(resultsNames[i])

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

실행하면 아래와 같은 그래프 3개가 나온다.

노름 변화, 상관계수는 이동량과 곡선으로 비례하고, 프로베니우스 거리는 직선으로 비례한다.

즉, 거리는 '차이'만 보기 때문에 선형적이고, 노름과 상관은 '합과 각도'를 보기 때문에 비선형이다.

행렬을 항등행렬 방향으로 조금 이동했을 때, 노름도 커지고 거리도 증가하지만, 방향도 달라질 수 있다.

5-5

임의의 계수를 가진 난수 행렬을 만든다.

계수가 r인 MxN 행렬을 만들려면 난수 Mxr 행렬에 rxN 행렬을 곱한다.

파이썬에서 구현하고 계수가 실제로 r인지 확인한다.

정답 코드

M = 5
N = 8
r = 3

A = np.random.randn(M,r) @ np.random.randn(r,N)

print(A.shape)
print(np.linalg.matrix_rank(A)) # 3
# Mxr 행렬도 계수가 3이하이고, rxN 행렬도 계수가 3이하임.

5-6

더한 결과가 계수-0, 계수-1, 계수-2 가 되는 세 쌍의 계수-1 행렬을 만들어 행렬 계수의 덧셈 법칙을 구현하기

A = np.diag([ 1,0,0,0,0])
B = np.diag([-1,0,0,0,0])
C = A+B

np.linalg.matrix_rank(A),np.linalg.matrix_rank(B),np.linalg.matrix_rank(C)

A = np.diag([1,0,0,0,0])
B = np.zeros(A.shape)
B[0,1] = 10
C = A+B

np.linalg.matrix_rank(A),np.linalg.matrix_rank(B),np.linalg.matrix_rank(C)

A = np.diag([1,0,0,0,0])
B = np.diag([0,1,0,0,0])
C = A+B

np.linalg.matrix_rank(A),np.linalg.matrix_rank(B),np.linalg.matrix_rank(C)

5-7

랭크가 정해진 두 행렬을 만들었을 때 합의 랭크/ 곱의 랭크가 어떻게 되는지 확인하는 문제

정답 코드

def makeAmatrix(M,r):
  return np.random.randn(M,r) @ np.random.randn(r,M)

matSize = 20 # matrix size (square)
rs = range(2,16) # range of ranks

Ranks = np.zeros((len(rs),len(rs),2))

for i in range(len(rs)):
  for j in range(len(rs)):
    S = makeAmatrix(matSize,rs[i]) + makeAmatrix(matSize,rs[j]) # A+B
    M = makeAmatrix(matSize,rs[i]) @ makeAmatrix(matSize,rs[j]) # A@B
    
    Ranks[i,j,0] = np.linalg.matrix_rank(S) # A+B의 랭크
    Ranks[i,j,1] = np.linalg.matrix_rank(M) # A@B의 랭크

    # 이론적으로 아래와 같다.
    # rank(S) ≤ min(20, rs[i] + rs[j])
    # rank(M) ≤ min(rank(A), rank(B))

    # 결과
    # - 덧셈 : rank(A)와 rank(B)가 둘다 커질수록 결과 랭크도 커짐.
    # 두 행렬을 더하면, 각자 가진 독립 정보가 거의 그대로 합쳐지되, 전체 차원(20)을 넘지는 못한다.
    # - 곱셈 : min(rank(A), rank(B)) 
    # 곱셈에서는 x축이 작거나 or y축이 작으면 결과도 작음.(작은 쪽이 병목)
    # 행렬을 곱하면, 정보량은 두 행렬 중 랭크가 작은 쪽으로 제한된다.

5-8

행렬 $A$, $A^T$, $A^TA$, $AA^T$의 계수가 모두 동일한 것 보이기.

내 코드

A = np.random.randn(6,6)
AT = A.T
ATA = AT @ A
AAT = A @ AT
print(np.linalg.matrix_rank(A)) # 6
print(np.linalg.matrix_rank(AT)) # 6
print(np.linalg.matrix_rank(ATA)) # 6
print(np.linalg.matrix_rank(AAT)) # 6

5-9

$v \in C(A)$의 답을 구하는 문제

• 정규분포에서 무작위로 추출한 숫자를 사용해 계수-3 행렬 $A \in R^{4*3}$과 벡터 $v \in R^4$ 만들기
• 벡터가 행렬의 열공간의 있는지 여부를 확인하기
• 코드를 여러번 실행해서 일관된 패턴을 보이는지 확인하기
• 다음으로는 $A \in R^{4*3}$의 계수-4 행렬을 사용
• A가 4x4 난수 행렬일 때 항상 $v \in C(A)$를 찾을 수 있다.

def is_V_inColA(A,v):

  if A.shape[0]!=v.shape[0]:
    raise Exception('Size mismatch! A and v must have the same column dimensionality!.')

  rankA  = np.linalg.matrix_rank(A)
  rankAv = np.linalg.matrix_rank( np.hstack((A,v)) ) # A 오른쪽에 v를 붙임

  # v가 A 열들의 조합이면 → rank 동일
  return rankA==rankAv

A = np.random.randn(4,3) # 4,3일때는 false, 4x4 이상일때는 true
v = np.random.randn(4,1)

is_V_inColA(A,v)

# 행렬 A가 full rank(랭크가 최댓값)일 때,
# 그 행렬의 열공간은 전체 공간이고,
# 따라서 그 공간의 모든 벡터가 A의 열공간에 포함된다.

5-10

• 정방 난수행렬 만들기
• 행렬의 계수를 줄이기(한 열을 다른 열의 배수로 설정)
• 행렬식을 계산하고 절댓값을 저장

이 세 단계를 이중 for루프에서 실행하기. 첫 번째 루프는 3x3에서 30x30까지 수행/두 번째 루프는 세 단계를 100회 반복

마지막으로 행렬의 원소 수에 대한 함수로 100회에 걸쳐 평균화된 행렬식을 선 차트로 그리기

ns = np.arange(3,31)
dets = np.zeros((len(ns),100))

for i in range(len(ns)):
  for j in range(100):
    A = np.random.randn(ns[i],ns[i])
    A[1] = A[0]*2 # 1번째와 2번째가 종속 > 따라서 이론 행렬식은 0이 되어야 한다.
    dets[i,j] = np.abs(np.linalg.det(A))

# 컴퓨터에서는 0이 아님. 그리고 행렬의 사이즈가 커질수록 그 크기가 점점 커진다.
# 큰 n일수록 연산 수 증가, 반올림 오차 폭증
# 컴퓨터에서는 n이 클수록 수치 계산에서 det로 rank / 가역성 판단하면 위험하다

# plotting
plt.figure(figsize=(6,4))
plt.plot(ns**2,np.log(np.mean(dets,axis=1)),'ks-',linewidth=3)
plt.xlabel('Number of elements in the matrix')
plt.ylabel('Log determinant')
plt.title('Empirical determinants of singular matrices')
plt.savefig('Figure_05_10.png',dpi=300)
plt.show()

원소수가 많아질수록 오차가 커진다.