대각합과 고윳값¶

대각합은 \(\operatorname{tr}(\mathbf{A}) = 3 + 3 = 6\)이다.

고윳값은 \(\det(\mathbf{A} - \lambda\mathbf{I}) = (3-\lambda)^2 - 1 = \lambda^2 - 6\lambda + 8 = 0\)을 만족하므로 \(\lambda_1 = 4\), \(\lambda_2 = 2\)이다.

고윳값의 합: \(4 + 2 = 6 = \operatorname{tr}(\mathbf{A})\). 또한 \(\det(\mathbf{A}) = 9 - 1 = 8 = 4 \times 2 = \lambda_1 \lambda_2\)이다.

\(\mathbf{A}\mathbf{B}\)의 \((i,i)\) 성분은 \(\sum_{k=1}^n a_{ik} b_{ki}\)이므로

이다. \(\mathbf{B}\mathbf{A}\)의 \((k,k)\) 성분은 \(\sum_{i=1}^m b_{ki} a_{ik}\)이므로

이다. 두 이중합 모두 \(i = 1, \dots, m\)과 \(k = 1, \dots, n\)에 대한 같은 항 \(a_{ik} b_{ki}\)를 더한 것이다. 덧셈의 교환법칙에 의해

이다. 참고로 \(\mathbf{A}\mathbf{B}\)는 \(m \times m\)이고 \(\mathbf{B}\mathbf{A}\)는 \(n \times n\)이다. 크기는 다를 수 있지만 대각합은 언제나 같다. \(\square\)

인자를 묶어 순환 성질을 적용한다.

순환 이동이 \(\mathbf{X}^\top\)를 맨 뒤에서 맨 앞으로 옮겨 \(p \times p\) 단위행렬을 만들어낸다. 이 결과는 \(n\)이나 \(\mathbf{X}\)의 구체적인 성분과 무관하게 성립한다.

완전 열계수 가정이 어디에 쓰였는지 짚어 두자. 그 가정 없이는 \(\mathbf{X}^\top\mathbf{X}\)가 가역이 아니어서 \(\mathbf{H}\) 자체가 정의되지 않는다. 계수가 \(r < p\)인 경우에는 일반화 역행렬을 써서 \(\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{+}\mathbf{X}^\top\)로 정의하며, 그때는 \(\operatorname{tr}(\mathbf{H}) = r\)이 된다.

\(\mathbf{z} \sim N(\mathbf{0}, \mathbf{I}_n)\)이므로 \(\boldsymbol{\mu} = \mathbf{0}\)이고, 항등식에 의해

이다. \(\mathbf{A}\)가 대칭 멱등이므로 고윳값은 모두 0 또는 1이고, 1인 고윳값의 개수가 계수 \(r\)이다. 따라서

이다. 이 결과가 \(\text{SSE}/\sigma^2 \sim \chi^2_{n-p}\)인 이유를 설명해 준다. 잔차생성행렬 \(\mathbf{M} = \mathbf{I} - \mathbf{H}\)의 계수가 \(n - p\)이므로 \(E[\text{SSE}/\sigma^2] = \operatorname{tr}(\mathbf{M}) = n - p\)이고, 이는 \(\chi^2_{n-p}\) 분포의 평균과 일치한다.

선형성. 대각 성분끼리의 덧셈이므로

이다.

곱은 그렇지 않다. 가장 간단한 반례는 \(\mathbf{A} = \mathbf{B} = \mathbf{I}_2\)다. \(\operatorname{tr}(\mathbf{A}\mathbf{B}) = \operatorname{tr}(\mathbf{I}_2) = 2\)이지만 \(\operatorname{tr}(\mathbf{A})\operatorname{tr}(\mathbf{B}) = 2 \times 2 = 4\)다. 아래 코드는 또 다른 반례를 보여 준다(\(5\) 대 \(0\)).

import numpy as np

A = np.array([[1., 2.], [3., 4.]])
B = np.array([[0., 1.], [1., 0.]])

print("tr(2A + 3B) =", np.trace(2*A + 3*B),
      " 2tr(A)+3tr(B) =", 2*np.trace(A) + 3*np.trace(B))
print("tr(AB) =", np.trace(A @ B),
      " tr(A)tr(B) =", np.trace(A) * np.trace(B))

출력:

tr(2A + 3B) = 10.0  2tr(A)+3tr(B) = 10.0
tr(AB) = 5.0  tr(A)tr(B) = 0.0

대각합은 덧셈에 대해서는 잘 행동하지만 곱에 대해서는 그렇지 않다. 곱에서 성립하는 것은 순환 성질 \(\operatorname{tr}(\mathbf{A}\mathbf{B}) = \operatorname{tr}(\mathbf{B}\mathbf{A})\)뿐이다. 참고로 행렬식은 정반대다. \(\det(\mathbf{A}\mathbf{B}) = \det(\mathbf{A})\det(\mathbf{B})\)는 성립하지만 \(\det(\mathbf{A}+\mathbf{B})\)에는 간단한 공식이 없다. \(\square\)

\(\mathbf{A}^\top\mathbf{A}\)의 \((j,j)\) 성분은 \(\sum_i a_{ij}^2\), 곧 \(j\)번째 열의 제곱합이다. 대각합은 이를 모든 열에 대해 더한 것이므로

이다. 한편 \(\mathbf{A}^\top\mathbf{A}\)의 고윳값은 특이값의 제곱 \(\sigma_i^2\)이고 대각합은 고윳값의 합이므로

이다.

import numpy as np

rng = np.random.default_rng(0)
A = rng.normal(size=(4, 4))
sv = np.linalg.svd(A, compute_uv=False)

print("tr(A^T A)   =", round(np.trace(A.T @ A), 6))
print("성분 제곱합 =", round((A ** 2).sum(), 6))
print("특이값 제곱합 =", round((sv ** 2).sum(), 6))

출력:

tr(A^T A)   = 13.498337
성분 제곱합 = 13.498337
특이값 제곱합 = 13.498337

통계적 의미. \(\mathbf{A}\)가 중심화된 자료행렬이면 \(\operatorname{tr}(\mathbf{A}^\top\mathbf{A})\)는 총제곱합이고, 특이값 제곱은 각 주성분이 설명하는 몫이다. "첫 \(k\)개 성분이 설명하는 비율"이 \(\sum_{i \le k}\sigma_i^2 / \sum_i \sigma_i^2\)인 것이 이 등식에서 나온다. \(\square\)

스칼라는 자기 자신의 대각합과 같다는 점(\(\mathbf{z}^\top\mathbf{A}\mathbf{z} = \operatorname{tr}(\mathbf{z}^\top\mathbf{A}\mathbf{z})\))에서 출발해 순환 성질을 쓴다.

대각합과 기댓값은 모두 선형이므로 순서를 바꿀 수 있다.

여기서 \(E[\mathbf{z}\mathbf{z}^\top] = \operatorname{Var}(\mathbf{z}) + E[\mathbf{z}]E[\mathbf{z}]^\top = \boldsymbol{\Sigma} + \boldsymbol{\mu}\boldsymbol{\mu}^\top\)이므로

이다(마지막에서 다시 순환 성질을 썼다).

\(\boldsymbol{\Sigma} = \mathbf{I}\)로 두면 본문의 공식이 된다. 이 유도에서 정규성은 전혀 쓰이지 않았다. 평균과 공분산만 있으면 성립한다. \(\square\)

\(\operatorname{Var}(\mathbf{Y}) = \mathbf{A}\boldsymbol{\Sigma}\mathbf{A}^\top\)이므로 순환 성질에 의해

이다. \(\mathbf{A}\)가 직교행렬이면 \(\mathbf{A}^\top\mathbf{A} = \mathbf{I}\)이므로

로 총분산이 보존된다.

import numpy as np

rng = np.random.default_rng(0)
Sigma = np.array([[4., 1., 0.], [1., 3., 1.], [0., 1., 2.]])

Q, _ = np.linalg.qr(rng.normal(size=(3, 3)))     # 직교행렬
S = np.diag([2., 1., 0.5])                        # 직교가 아닌 대각 척도변환

print("tr(Sigma)        =", round(np.trace(Sigma), 6))
print("직교 Q 로 변환   =", round(np.trace(Q @ Sigma @ Q.T), 6))
print("척도 S 로 변환   =", round(np.trace(S @ Sigma @ S.T), 6))

출력:

tr(Sigma)        = 9.0
직교 Q 로 변환   = 9.0
척도 S 로 변환   = 19.5

회전은 총분산을 보존하지만 척도변환은 그렇지 않다. PCA가 회전만 하는 이유가 여기에 있다. 분산의 총량은 그대로 두고 축 사이의 배분만 바꾼다. 반대로 변수를 표준화하는 것은 척도변환이므로 총분산이 바뀐다(표준화하면 \(\operatorname{tr} = p\)가 된다). 공분산행렬로 PCA를 하는 것과 상관행렬로 하는 것이 다른 결과를 주는 이유다. \(\square\)

실계수 특성다항식의 복소근은 켤레쌍으로 나타난다. \(\lambda = a + bi\)가 근이면 \(\bar{\lambda} = a - bi\)도 근이다. 켤레쌍끼리 더하고 곱하면

로 허수부가 상쇄된다. 따라서 전체 합과 곱이 실수다.

import numpy as np

R = np.array([[0., -1.], [1., 0.]])       # 90도 회전
ev = np.linalg.eigvals(R)

print("고윳값:", ev)
print("합    :", ev.sum().real, " tr(R) =", np.trace(R))
print("곱    :", np.prod(ev).real, " det(R) =", round(np.linalg.det(R), 6))

출력:

고윳값: [0.+1.j 0.-1.j]
합    : 0.0  tr(R) = 0.0
곱    : 1.0  det(R) = 1.0

\(90^\circ\) 회전행렬은 고윳값이 \(\pm i\)로 순허수다. 실벡터 중 방향이 보존되는 것이 하나도 없으므로 당연하다. 그럼에도 합은 \(0 = \operatorname{tr}(\mathbf{R})\), 곱은 \(1 = \det(\mathbf{R})\)로 실수다.

대칭행렬에서는 이런 일이 없다. 스펙트럼 정리가 고윳값이 모두 실수임을 보장한다. 통계에서 다루는 공분산행렬과 사영행렬이 모두 대칭이므로, 실무에서 복소 고윳값을 만날 일은 드물다. \(\square\)

특이값분해 \(\mathbf{X} = \mathbf{U}\mathbf{D}\mathbf{V}^\top\)를 넣으면 \(\mathbf{X}^\top\mathbf{X} = \mathbf{V}\mathbf{D}^2\mathbf{V}^\top\)이고 \(\mathbf{V}^\top\mathbf{V} = \mathbf{I}\)이므로

이다. 따라서

이고, 순환 성질과 \(\mathbf{U}^\top\mathbf{U} = \mathbf{I}\)에서

를 얻는다.

import numpy as np

rng = np.random.default_rng(0)
X = rng.normal(size=(50, 5))
d = np.linalg.svd(X, compute_uv=False)

for lam in (0., 1., 10., 100.):
    H = X @ np.linalg.inv(X.T @ X + lam * np.eye(5)) @ X.T
    print(f"lambda={lam:>6}:  tr(H) = {np.trace(H):7.4f}"
          f"   sum d^2/(d^2+lambda) = {np.sum(d**2/(d**2+lam)):7.4f}")

출력:

lambda=   0.0:  tr(H) =  5.0000   sum d^2/(d^2+lambda) =  5.0000
lambda=   1.0:  tr(H) =  4.8942   sum d^2/(d^2+lambda) =  4.8942
lambda=  10.0:  tr(H) =  4.1267   sum d^2/(d^2+lambda) =  4.1267
lambda= 100.0:  tr(H) =  1.6690   sum d^2/(d^2+lambda) =  1.6690

\(\lambda = 0\)이면 \(\operatorname{tr} = p = 5\)로 보통최소제곱과 같고, \(\lambda\)가 커질수록 값이 줄어 \(\lambda \to \infty\)에서 0으로 간다.

이 값을 유효자유도라 부른다. 능형회귀는 모수를 \(5\)개 그대로 두지만 벌점이 각 방향을 \(d_j^2/(d_j^2+\lambda)\)만큼 축소하므로, 모형이 실제로 쓰는 자유도는 그보다 작다. 특이값이 작은 방향(공선성이 심한 방향)일수록 더 강하게 축소된다는 점도 식에서 바로 읽힌다. 18장에서 이 양이 모형 선택 기준에 쓰인다. \(\square\)