Lecture 17. 직교행렬과 그람-슈미트 — 파이썬 실습 - 선형대수

L17 서술 파트의 결론은 두 줄이었다. 직교기저에서는 계수가 내적 한 번으로 나오고, 아무 벡터나 직교로 만드는 절차가 그람-슈미트이며 그 기록이 $A = QR$ 이다.

이 노트북에서는 그람-슈미트를 단계마다 기록하는 함수를 만들어 애니메이션으로 보고, $QR$ 로 푼 최소제곱이 정규방정식보다 실제로 얼마나 안전한지 숫자로 잰다.

서술 파트의 내용	여기서 확인하는 방법
$Q^{\mathsf{T}}Q = I$ 인데 $QQ^{\mathsf{T}} \neq I$	직사각 $Q$ 로 두 곱을 다 계산
길이와 각도 보존	회전행렬을 곱해 보고 재기
계수 $=$ 내적 한 번	연립방정식 풀이와 연산 횟수 비교
그람-슈미트	단계 기록 + 슬라이더 애니메이션
$R$ 이 상삼각	대각 아래가 왜 0인지 성분으로
$QR$ 이 더 안전하다	고차 다항식 맞추기로 오차 비교

0. 준비¶

import numpy as np

from linalg_viz import COLORS, arrow, layout3d, show_matrix, slider_figure

np.set_printoptions(precision=3, suppress=True)
rng = np.random.default_rng(17)
print("numpy", np.__version__)

numpy 2.5.2

1. $Q^{\mathsf{T}}Q$ 와 $QQ^{\mathsf{T}}$ 는 다르다¶

가장 흔한 오해부터 짚는다. 직사각 $Q$ 에서는 한쪽만 단위행렬이다.

Q직사각 = np.array([[1.0, 0.0],
                  [0.0, 1.0],
                  [0.0, 0.0]])            # 3 x 2, 열이 정규직교

print(show_matrix(Q직사각.T @ Q직사각, "Q^T Q  (2 x 2)"))
print(show_matrix(Q직사각 @ Q직사각.T, "Q Q^T  (3 x 3)"))
print("Q^T Q = I 인가 :", np.allclose(Q직사각.T @ Q직사각, np.eye(2)))
print("Q Q^T = I 인가 :", np.allclose(Q직사각 @ Q직사각.T, np.eye(3)))
print("rank(Q Q^T) :", np.linalg.matrix_rank(Q직사각 @ Q직사각.T),
      "  -> 3 이 될 수 없다")

Q^T Q  (2 x 2)
[  1   0 ]
[  0   1 ]
Q Q^T  (3 x 3)
[  1   0   0 ]
[  0   1   0 ]
[  0   0   0 ]
Q^T Q = I 인가 : True
Q Q^T = I 인가 : False
rank(Q Q^T) : 2   -> 3 이 될 수 없다

정사각이면 둘 다 $I$ 가 되고, 그때 $Q^{-1} = Q^{\mathsf{T}}$ 이다. 회전행렬로 확인하자.

각도 = np.pi / 6
Q회전 = np.array([[np.cos(각도), -np.sin(각도)],
               [np.sin(각도),  np.cos(각도)]])

print(show_matrix(Q회전, "회전행렬 (30도)"))
print("Q^T Q = I :", np.allclose(Q회전.T @ Q회전, np.eye(2)))
print("Q Q^T = I :", np.allclose(Q회전 @ Q회전.T, np.eye(2)))
print("Q^-1 이 Q^T 와 같은가 :", np.allclose(np.linalg.inv(Q회전), Q회전.T))

회전행렬 (30도)
[  0.866    -0.5 ]
[    0.5   0.866 ]
Q^T Q = I : True
Q Q^T = I : True
Q^-1 이 Q^T 와 같은가 : True

길이도 각도도 그대로이다.

x = np.array([3.0, 1.0])
y = np.array([-1.0, 2.0])

print("|x|     =", np.linalg.norm(x), "   |Qx|     =", np.linalg.norm(Q회전 @ x))
print("x . y   =", x @ y, "         (Qx).(Qy) =", (Q회전 @ x) @ (Q회전 @ y))
print()
print("길이 보존 :", np.isclose(np.linalg.norm(x), np.linalg.norm(Q회전 @ x)))
print("내적 보존 :", np.isclose(x @ y, (Q회전 @ x) @ (Q회전 @ y)))

|x|     = 3.1622776601683795    |Qx|     = 3.162277660168379
x . y   = -1.0          (Qx).(Qy) = -1.0000000000000004

길이 보존 : True
내적 보존 : True

2. 계수는 내적 한 번이면 나온다¶

같은 벡터를 일반 기저와 정규직교 기저로 각각 표현해 보고, 드는 계산의 양을 비교한다.

n = 600
일반기저 = rng.normal(size=(n, n))                     # 아무 기저
직교기저, _ = np.linalg.qr(rng.normal(size=(n, n)))     # 정규직교 기저
b = rng.normal(size=n)

계수_일반 = np.linalg.solve(일반기저, b)                 # 연립방정식을 푼다
계수_직교 = 직교기저.T @ b                               # 내적 n 번

print("둘 다 b 를 제대로 복원하는가 :",
      np.allclose(일반기저 @ 계수_일반, b),
      np.allclose(직교기저 @ 계수_직교, b))

둘 다 b 를 제대로 복원하는가 : True True

답은 둘 다 맞다. 다른 것은 드는 계산의 양이다. 소거는 곱셈과 덧셈이 대략 $n^3/3$ 번 필요하지만(L4), 내적 $n$ 번은 $n^2$ 번이면 끝난다.

시계로 재면 기계와 라이브러리 설정에 따라 값이 크게 흔들리므로, 여기서는 연산 횟수를 센다.

print(f"{'n':>7}{'소거 n^3/3':>16}{'내적 n^2':>14}{'비율':>10}")
for n in (100, 1000, 10000):
    소거 = n ** 3 / 3
    내적 = n ** 2
    print(f"{n:>7}{소거:>16,.0f}{내적:>14,.0f}{소거 / 내적:>10,.0f}")
print()
print("n 이 열 배가 되면 격차도 열 배가 벌어진다.")

      n        소거 n^3/3        내적 n^2        비율
    100         333,333        10,000        33
   1000     333,333,333     1,000,000       333
  10000 333,333,333,333   100,000,000     3,333

n 이 열 배가 되면 격차도 열 배가 벌어진다.

직교기저에서는 계수 하나만 따로 구할 수도 있다. 다른 계수를 몰라도 된다.

j = 42
print("전체를 구한 뒤 j 번째  :", 계수_직교[j])
print("j 번째만 내적 한 번    :", 직교기저[:, j] @ b)
print("같은가 :", np.isclose(계수_직교[j], 직교기저[:, j] @ b))

전체를 구한 뒤 j 번째  : 0.14676699553575634
j 번째만 내적 한 번    : 0.14676699553575706
같은가 : True

3. 그람-슈미트¶

각 단계를 기록해 두면 나중에 애니메이션 재료로 쓸 수 있다.

def 그람슈미트(A, 기록=None):
    """열을 차례로 직교화한다. 기록 리스트를 주면 단계마다 중간값을 담는다."""
    A = np.asarray(A, dtype=float)
    Q = []
    for j in range(A.shape[1]):
        v = A[:, j].copy()
        뺀것 = np.zeros_like(v)
        for q in Q:
            뺀것 = 뺀것 + (q @ v) * q       # 이미 만든 방향들의 성분
        v = v - 뺀것
        if 기록 is not None:
            기록.append(dict(열=j, 원본=A[:, j].copy(), 뺀양=뺀것, 남은것=v.copy()))
        Q.append(v / np.linalg.norm(v))     # 정규화
    return np.column_stack(Q)

a1 = np.array([1.0, 1.0, 1.0])
a2 = np.array([1.0, 0.0, 2.0])
a3 = np.array([-1.0, 1.0, 3.0])
A = np.column_stack([a1, a2, a3])

기록 = []
Q = 그람슈미트(A, 기록)

print(show_matrix(Q, "Q"))
print(show_matrix(Q.T @ Q, "Q^T Q"))
print("단위행렬인가 :", np.allclose(Q.T @ Q, np.eye(3), atol=1e-12))

Q
[      0.577   -1.57e-16      -0.816 ]
[      0.577      -0.707       0.408 ]
[      0.577       0.707       0.408 ]
Q^T Q
[          1   -2.68e-16   -3.78e-18 ]
[  -2.68e-16           1     2.7e-16 ]
[  -3.78e-18     2.7e-16           1 ]
단위행렬인가 : True

서술 파트에서 손으로 구한 것과 비교해 보자. $\vv{q}_1 = (1,1,1)/\sqrt3$ , $\vv{q}_2 = (0,-1,1)/\sqrt2$ , $\vv{q}_3 = (-2,1,1)/\sqrt6$ 이었다.

손으로 = np.column_stack([
    np.array([1.0, 1, 1]) / np.sqrt(3),
    np.array([0.0, -1, 1]) / np.sqrt(2),
    np.array([-2.0, 1, 1]) / np.sqrt(6),
])
print(show_matrix(손으로, "손으로 구한 Q"))
print("코드와 같은가 :", np.allclose(Q, 손으로))
print()
for 단계 in 기록:
    print(f"열 {단계['열'] + 1} : 원본 {단계['원본']}  -  뺀양 {단계['뺀양']}"
          f"  =  남은것 {단계['남은것']}")

손으로 구한 Q
[   0.577        0   -0.816 ]
[   0.577   -0.707    0.408 ]
[   0.577    0.707    0.408 ]
코드와 같은가 : True

열 1 : 원본 [1. 1. 1.]  -  뺀양 [0. 0. 0.]  =  남은것 [1. 1. 1.]
열 2 : 원본 [1. 0. 2.]  -  뺀양 [1. 1. 1.]  =  남은것 [-0. -1.  1.]
열 3 : 원본 [-1.  1.  3.]  -  뺀양 [1. 0. 2.]  =  남은것 [-2.  1.  1.]

4. 단계별로 보기¶

슬라이더를 움직이며 벡터가 하나씩 직교로 바뀌는 것을 보자. 주황이 원본, 파랑이 빼는 성분, 초록이 남은 것, 청록이 완성된 $\vv{q}$ 이다.

def 프레임(화살표들, 슬롯=5):
    """(시작, 끝, 색, 이름) 목록으로 트레이스를 만든다. 개수는 항상 같게 채운다."""
    채움 = list(화살표들) + [([0, 0, 0], [0, 0, 0], COLORS["basis"], "")] * (
        슬롯 - len(화살표들))
    트레이스 = []
    for 시작, 끝, 색, 이름 in 채움:
        트레이스 += arrow(시작, 끝, 색, 이름, legend=bool(이름))
    return 트레이스

원본색, 뺄색, 남은색, 완성색 = (COLORS["output"], COLORS["input"],
                          COLORS["second"], COLORS["colspace"])
q1, q2, q3 = Q[:, 0], Q[:, 1], Q[:, 2]
뺀양2, 남은2 = 기록[1]["뺀양"], 기록[1]["남은것"]
뺀양3, 남은3 = 기록[2]["뺀양"], 기록[2]["남은것"]
영 = [0, 0, 0]

장면 = [
    ("1. 원본 세 벡터", [(영, a1, 원본색, "a1"), (영, a2, 원본색, "a2"),
                    (영, a3, 원본색, "a3")]),
    ("2. q1 = a1 / |a1|", [(영, q1, 완성색, "q1"), (영, a2, 원본색, "a2"),
                          (영, a3, 원본색, "a3")]),
    ("3. a2 에서 q1 성분을 뺀다", [(영, q1, 완성색, "q1"), (영, a2, 원본색, "a2"),
                            (영, 뺀양2, 뺄색, "뺄 성분"),
                            (뺀양2, a2, 남은색, "남은 것")]),
    ("4. q2 = 남은 것 / 길이", [(영, q1, 완성색, "q1"), (영, q2, 완성색, "q2"),
                          (영, a3, 원본색, "a3")]),
    ("5. a3 에서 평면 성분을 뺀다", [(영, q1, 완성색, "q1"), (영, q2, 완성색, "q2"),
                            (영, a3, 원본색, "a3"), (영, 뺀양3, 뺄색, "뺄 성분"),
                            (뺀양3, a3, 남은색, "남은 것")]),
    ("6. q3 까지 완성", [(영, q1, 완성색, "q1"), (영, q2, 완성색, "q2"),
                    (영, q3, 완성색, "q3")]),
]

slider_figure([프레임(화살표) for _, 화살표 in 장면],
              [이름 for 이름, _ in 장면],
              layout3d("그람-슈미트 : 빼고, 남은 것을 쓴다", extent=3.2),
              initial=0)

각 단계의 "남은 것"이 정말 앞의 $\vv{q}$ 들과 직교하는지 숫자로 확인하자.

print("3단계의 남은 것 :", 남은2)
print("  q1 과의 내적 :", f"{q1 @ 남은2:+.2e}")
print()
print("5단계의 남은 것 :", 남은3)
print("  q1 과의 내적 :", f"{q1 @ 남은3:+.2e}")
print("  q2 와의 내적 :", f"{q2 @ 남은3:+.2e}")

3단계의 남은 것 : [-0. -1.  1.]
  q1 과의 내적 : -3.33e-16

5단계의 남은 것 : [-2.  1.  1.]
  q1 과의 내적 : +0.00e+00
  q2 와의 내적 : +6.66e-16

5. $A = QR$ ¶

R = Q.T @ A

print(show_matrix(R, "R = Q^T A"))
print("Q R 이 A 와 같은가 :", np.allclose(Q @ R, A))
print("대각 아래가 0 인가 :", np.allclose(np.tril(R, -1), 0, atol=1e-12))
print()
print("대각 성분 :", np.diag(R))
print("각 단계의 |남은 것| :", [f"{np.linalg.norm(단계['남은것']):.3f}" for 단계 in 기록])

R = Q^T A
[       1.73        1.73        1.73 ]
[  -4.44e-16        1.41        1.41 ]
[  -5.55e-17    3.33e-16        2.45 ]
Q R 이 A 와 같은가 : True
대각 아래가 0 인가 : True

대각 성분 : [1.732 1.414 2.449]
각 단계의 |남은 것| : ['1.732', '1.414', '2.449']

$R_{ij} = \vv{q}_i \cdot \vv{a}_j$ 이고, $i > j$ 이면 0인 이유는 $\vv{a}_j$ 가 $\vv{q}_1, \dots, \vv{q}_j$ 의 결합이기 때문이었다. 그것도 확인해 보자.

for j in range(3):
    앞부분 = Q[:, :j + 1] @ (Q[:, :j + 1].T @ A[:, j])   # q1..qj 로만 조립
    print(f"a{j + 1} 을 q1..q{j + 1} 로만 조립 : {앞부분}"
          f"   원본과 같은가 {np.allclose(앞부분, A[:, j])}")
print()
for i in range(3):
    for j in range(3):
        if i > j:
            print(f"  R[{i + 1}][{j + 1}] = q{i + 1} . a{j + 1}"
                  f" = {Q[:, i] @ A[:, j]:+.1e}")

a1 을 q1..q1 로만 조립 : [1. 1. 1.]   원본과 같은가 True
a2 을 q1..q2 로만 조립 : [1. 0. 2.]   원본과 같은가 True
a3 을 q1..q3 로만 조립 : [-1.  1.  3.]   원본과 같은가 True

  R[2][1] = q2 . a1 = -4.4e-16
  R[3][1] = q3 . a1 = -5.6e-17
  R[3][2] = q3 . a2 = +3.3e-16

numpy 의 qr 과 비교하면 부호가 다를 수 있다. $\vv{q}$ 의 방향을 반대로 잡아도 직교기저이기 때문이다.

Qn, Rn = np.linalg.qr(A)

print(show_matrix(Qn, "numpy 의 Q"))
print("우리 Q 와 부호만 다른가 :",
      np.allclose(np.abs(Qn), np.abs(Q)))
print("두 분해 모두 A 를 복원하는가 :",
      np.allclose(Qn @ Rn, A), np.allclose(Q @ R, A))
print()
부호 = np.sign(np.diag(Qn.T @ Q))
print("열마다의 부호 :", 부호)
print("부호를 맞추면 같은가 :", np.allclose(Qn * 부호, Q))

numpy 의 Q
[     -0.577   -1.22e-16      -0.816 ]
[     -0.577      -0.707       0.408 ]
[     -0.577       0.707       0.408 ]
우리 Q 와 부호만 다른가 : True
두 분해 모두 A 를 복원하는가 : True True

열마다의 부호 : [-1.  1.  1.]
부호를 맞추면 같은가 : True

6. $QR$ 이 왜 더 안전한가¶

고차 다항식을 맞추면 열들이 서로 닮아 간다. 그때 정규방정식과 $QR$ 의 차이가 드러난다.

print(f"{'차수':>4}{'cond(A)':>12}{'cond(A^T A)':>14}"
      f"{'정규방정식 오차':>16}{'QR 오차':>12}")
for 차수 in (4, 8, 12):
    t = np.linspace(0.0, 1.0, 30)
    Ad = np.column_stack([t ** k for k in range(차수 + 1)])
    참계수 = np.ones(차수 + 1)
    bd = Ad @ 참계수                                   # 잡음 없는 정확한 우변

    정규 = np.linalg.solve(Ad.T @ Ad, Ad.T @ bd)        # A^T A 를 만든다
    Qd, Rd = np.linalg.qr(Ad)
    qr해 = np.linalg.solve(Rd, Qd.T @ bd)              # A^T A 를 만들지 않는다

    print(f"{차수:>4}{np.linalg.cond(Ad):>12.2e}{np.linalg.cond(Ad.T @ Ad):>14.2e}"
          f"{np.abs(정규 - 참계수).max():>16.2e}{np.abs(qr해 - 참계수).max():>12.2e}")

  차수     cond(A)   cond(A^T A)        정규방정식 오차       QR 오차
   4    6.16e+02      3.80e+05        2.88e-12    1.07e-14
   8    6.08e+05      3.70e+11        3.37e-07    3.34e-11
  12    7.48e+08      1.26e+17        2.80e+00    3.14e-09

cond(A^T A) 가 cond(A) 의 제곱인 것이 보인다. 차수 12에서는 정규방정식의 답이 아예 못 쓸 정도가 되는데, $QR$ 은 여전히 쓸 만하다.

차수 = 12
t = np.linspace(0.0, 1.0, 30)
Ad = np.column_stack([t ** k for k in range(차수 + 1)])
참계수 = np.ones(차수 + 1)
bd = Ad @ 참계수

정규 = np.linalg.solve(Ad.T @ Ad, Ad.T @ bd)
Qd, Rd = np.linalg.qr(Ad)
qr해 = np.linalg.solve(Rd, Qd.T @ bd)

print("참 계수는 전부 1 이다.")
print("정규방정식 :", np.round(정규[:6], 3), "...")
print("QR         :", np.round(qr해[:6], 6), "...")
print()
print(f"cond(A)   = {np.linalg.cond(Ad):.2e}")
print(f"cond(A)^2 = {np.linalg.cond(Ad) ** 2:.2e}")
print(f"cond(A^T A) = {np.linalg.cond(Ad.T @ Ad):.2e}   <- 거의 제곱이다")

참 계수는 전부 1 이다.
정규방정식 : [1.    1.    1.    1.005 0.956 1.235] ...
QR         : [1. 1. 1. 1. 1. 1.] ...

cond(A)   = 7.48e+08
cond(A)^2 = 5.60e+17
cond(A^T A) = 1.26e+17   <- 거의 제곱이다

7. 그람-슈미트 자체는 안전하지 않다¶

직교기저가 좋다는 것과 그것을 그람-슈미트로 구하는 것이 안전하다는 것은 다른 이야기이다. 서술 파트에서 경고한 대로 확인해 보자.

print(f"{'eps':>8}{'고전 GS 의 |Q^T Q - I|':>24}{'numpy.qr 의 |Q^T Q - I|':>26}")
for eps in (1e-4, 1e-7, 1e-10):
    L = np.array([[1.0, 1.0, 1.0],
                  [eps, 0.0, 0.0],
                  [0.0, eps, 0.0],
                  [0.0, 0.0, eps]])            # 열이 거의 같다
    Qg = 그람슈미트(L)
    Qn2, _ = np.linalg.qr(L)
    print(f"{eps:>8.0e}"
          f"{np.abs(Qg.T @ Qg - np.eye(3)).max():>24.2e}"
          f"{np.abs(Qn2.T @ Qn2 - np.eye(3)).max():>26.2e}")

     eps     고전 GS 의 |Q^T Q - I|    numpy.qr 의 |Q^T Q - I|
   1e-04                2.26e-09                  2.22e-16
   1e-07                1.33e-02                  4.44e-16
   1e-10                5.00e-01                  3.33e-16

열이 거의 같을수록 고전 그람-슈미트의 직교성이 무너진다. numpy.linalg.qr 은 그람-슈미트가 아니라 하우스홀더 반사를 쓰므로 계속 기계 정밀도를 지킨다. 공식을 아는 것과 그 공식을 그대로 코드로 옮기는 것은 다른 일이다.

마치며...¶

서술 파트의 내용	이 노트북의 코드
$Q^{\mathsf{T}}Q = I$ , $QQ^{\mathsf{T}} \neq I$	직사각 $Q$ 로 두 곱을 계산
길이·각도 보존	회전행렬로 확인
계수 $=$ 내적	소거 $n^3/3$ 대 내적 $n^2$
그람-슈미트	`그람슈미트(A, 기록)` 과 슬라이더
$R$ 이 상삼각	$\vv{a}_j$ 가 $\vv{q}_1..\vv{q}_j$ 로 조립되는 것 확인
$QR$ 이 안전	차수 12에서 정규방정식은 실패, $QR$ 은 성공
고전 GS 의 불안정	라우클리 행렬로 직교성 손실 측정

더 해 볼 것¶

그람슈미트 에 열의 순서를 바꿔 넣어 보자. $Q$ 가 달라지는가? $QR$ 은 여전히 $A$ 인가?
3절의 $A$ 에 종속인 열을 하나 넣으면 어떻게 되는가? 어디에서 무엇이 터지는가?
수정 그람-슈미트(각 단계에서 남은 벡터를 바로 갱신하며 빼는 방식)를 짜서 7절의 표에 한 줄 더해 보자. 고전 방식보다 나아지는가?
6절에서 $\vv{b}$ 에 작은 잡음을 넣어 보자. 두 방법의 차이가 더 벌어지는가?

여기까지가 전반부이다. 다음 강의부터는 행렬 자체를 들여다본다. 이 행렬은 무엇을 하는가, 바꾸지 않는 방향은 무엇인가를 묻기 시작한다.

Lecture 17. 직교행렬과 그람-슈미트 — 파이썬 실습