Lecture 7. Ax = 0 풀기 — 파이썬 실습 - 선형대수

L7 서술 파트에서 소거로 영공간을 구하는 절차를 다루었다. 이 노트북에서는 그 절차를 함수로 짜고, sympy 와 scipy 의 결과와 대조해 보도록 하자.

서술 파트의 내용	여기서 확인하는 방법
피벗이 없으면 그 열을 건너뛴다	`row_echelon` 을 직접 짠다
랭크 $=$ 피벗의 개수	손으로 센 값과 `matrix_rank` 를 대조한다
특수해 구하는 절차	`special_solutions` 를 짜고 $A\vv{s} = \vv{0}$ 을 검산한다
$N = \begin{bmatrix} -F \\ I \end{bmatrix}$	`sympy` 의 기약 사다리꼴에서 직접 조립한다
$\dim N(A) = n - r$	랭크를 낮춰 가며 영공간이 자라는 것을 본다

0. 준비¶

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, plane, spin_figure, show_matrix

np.set_printoptions(precision=3, suppress=True)
print("numpy", np.__version__)

numpy 2.5.2

A = np.array([[1, 2, 2,  2],
              [2, 4, 6,  8],
              [3, 6, 8, 10]], dtype=float)

print(show_matrix(A, "A ="))
print("크기 :", A.shape, " -> 식 3개, 미지수 4개")

A =
[   1    2    2    2 ]
[   2    4    6    8 ]
[   3    6    8   10 ]
크기 : (3, 4)  -> 식 3개, 미지수 4개

1. 사다리꼴 만들기¶

L2의 소거 함수와 다른 점은 하나이다. 피벗 자리가 0이고 아래에도 0이 아닌 값이 없으면 그 열을 건너뛰고 오른쪽으로 간다. 어느 열에서 피벗을 찾았는지 함께 기록해 둔다.

def row_echelon(A, tol=1e-10):
    """소거해서 사다리꼴을 만들고 (U, 피벗열 목록) 을 돌려준다."""
    U = np.asarray(A, dtype=float).copy()
    m, n = U.shape
    피벗열, 행 = [], 0

    for 열 in range(n):
        if 행 >= m:
            break

        # 이 열에서 0이 아닌 성분을 가진 행을 찾는다
        후보 = None
        for i in range(행, m):
            if abs(U[i, 열]) > tol:
                후보 = i
                break
        if 후보 is None:
            continue                       # 피벗이 없다 -> 자유 열

        if 후보 != 행:
            U[[행, 후보]] = U[[후보, 행]]   # 행 교환

        for 아래 in range(행 + 1, m):
            곱수 = U[아래, 열] / U[행, 열]
            U[아래] = U[아래] - 곱수 * U[행]

        피벗열.append(열)
        행 += 1

    U[np.abs(U) < tol] = 0.0               # 반올림 찌꺼기를 정리한다
    return U, 피벗열

U, 피벗열 = row_echelon(A)

print(show_matrix(U, "U ="))
print("피벗이 있는 열 (0부터) :", 피벗열)
print("피벗 값 :", [float(U[i, j]) for i, j in enumerate(피벗열)])

U =
[  1   2   2   2 ]
[  0   0   2   4 ]
[  0   0   0   0 ]
피벗이 있는 열 (0부터) : [0, 2]
피벗 값 : [1.0, 2.0]

서술 파트에서 손으로 구한 것과 같다. 두 번째 열에는 피벗이 없어서 건너뛰었고, 마지막 행이 통째로 0이 되었다.

2. 랭크와 자유 열¶

m, n = A.shape
r = len(피벗열)
자유열 = [j for j in range(n) if j not in 피벗열]

print(f"rank r      = {r}")
print(f"자유 열     = {자유열}   ({n} - {r} = {n - r} 개)")
print()
print("np.linalg.matrix_rank 와 같은가 :", r == np.linalg.matrix_rank(A))

rank r      = 2
자유 열     = [1, 3]   (4 - 2 = 2 개)

np.linalg.matrix_rank 와 같은가 : True

랭크는 행의 개수 3보다 작다. 식을 세 개 적었지만 실제로 서로 다른 이야기를 하는 식은 둘뿐이다. 세 번째 행이 어떻게 만들어지는지 확인해 보자.

print("row1 + row2 =", A[0] + A[1])
print("row3        =", A[2])
print("같은가 :", np.allclose(A[0] + A[1], A[2]))

row1 + row2 = [ 3.  6.  8. 10.]
row3        = [ 3.  6.  8. 10.]
같은가 : True

3. 특수해 구하기¶

자유 변수 하나에만 1을 주고 나머지에 0을 준 다음, 피벗 행을 아래에서 위로 훑으며 피벗 변수를 확정한다. 사다리꼴에서 피벗 행은 그 피벗 변수와 오른쪽 변수들만 담고 있으므로 역순으로 올라가면 하나씩 정해진다.

def special_solutions(A, tol=1e-10):
    """자유 변수마다 특수해를 하나씩 만들어 열로 쌓는다."""
    U, 피벗열 = row_echelon(A, tol)
    n = U.shape[1]
    자유열 = [j for j in range(n) if j not in 피벗열]

    해들 = []
    for 자유 in 자유열:
        x = np.zeros(n)
        x[자유] = 1.0
        # 피벗을 오른쪽부터 왼쪽으로 확정한다
        for 행, 열 in reversed(list(enumerate(피벗열))):
            x[열] = -(U[행, 열 + 1:] @ x[열 + 1:]) / U[행, 열]
        해들.append(x + 0.0)               # -0.0 을 0.0 으로 정리한다

    return np.column_stack(해들) if 해들 else np.zeros((n, 0))

S = special_solutions(A)

print(show_matrix(S, "특수해 (열마다 하나)"))
print()
for j in range(S.shape[1]):
    print(f"s{j + 1} = {S[:, j]}")

특수해 (열마다 하나)
[  -2    2 ]
[   1    0 ]
[   0   -2 ]
[   0    1 ]

s1 = [-2.  1.  0.  0.]
s2 = [ 2.  0. -2.  1.]

서술 파트에서 손으로 구한 $(-2, 1, 0, 0)$ 과 $(2, 0, -2, 1)$ 이 그대로 나왔다. 검산해 보자.

print(show_matrix(A @ S, "A @ S ="))
print("모두 영벡터인가 :", np.allclose(A @ S, 0))

A @ S =
[  0   0 ]
[  0   0 ]
[  0   0 ]
모두 영벡터인가 : True

특수해들의 아무 선형결합이나 만들어도 영공간에 있어야 한다.

rng = np.random.default_rng(0)
for _ in range(3):
    c = rng.standard_normal(S.shape[1])
    x = S @ c
    print(f"c = {np.round(c, 3)}  ->  A x = {A @ x}")

c = [ 0.126 -0.132]  ->  A x = [0. 0. 0.]
c = [0.64  0.105]  ->  A x = [0. 0. 0.]
c = [-0.536  0.362]  ->  A x = [0. 0. 0.]

4. `scipy` 와 대조하기¶

scipy.linalg.null_space 는 길이가 1이고 서로 직교하는 기저를 돌려준다. 우리가 구한 것과 벡터 자체는 다르지만 만들어내는 공간은 같아야 한다.

두 기저를 나란히 붙여 놓고 랭크를 재면 확인할 수 있다. 같은 공간을 만든다면 붙여도 랭크가 늘지 않는다.

from scipy.linalg import null_space

N_sp = null_space(A)
print(show_matrix(N_sp, "scipy 의 기저"))
print()
print("우리 기저의 열 개수 :", S.shape[1])
print("scipy 기저의 열 개수 :", N_sp.shape[1])
print("둘을 붙였을 때 랭크  :", np.linalg.matrix_rank(np.column_stack([S, N_sp])))
print("같은 공간을 만드는가 :",
      np.linalg.matrix_rank(np.column_stack([S, N_sp])) == S.shape[1])

scipy 의 기저
[    0.907     0.069 ]
[   -0.409     0.378 ]
[  -0.0893    -0.826 ]
[   0.0446     0.413 ]

우리 기저의 열 개수 : 2
scipy 기저의 열 개수 : 2
둘을 붙였을 때 랭크  : 2
같은 공간을 만드는가 : True

5. 기약 사다리꼴에서 영공간 행렬 읽기¶

sympy 는 분수를 정확히 다루므로 기약 사다리꼴을 깔끔하게 구해 준다.

import sympy as sp

R_sym, 피벗_sym = sp.Matrix(A.astype(int)).rref()
피벗열_sym = list(피벗_sym)

print("R =")
sp.pprint(R_sym)
print()
print("피벗 열 :", 피벗열_sym)

R =
⎡1  2  0  -2⎤
⎢           ⎥
⎢0  0  1  2 ⎥
⎢           ⎥
⎣0  0  0  0 ⎦

피벗 열 : [0, 2]

서술 파트의 (24)에 나온 $R$ 과 같다. 이제 $F$ 를 꺼내 $N = \begin{bmatrix} -F \\ I \end{bmatrix}$ 를 조립한다. 열 순서를 실제로 바꿀 필요는 없고, 피벗 행에 $-F$ 를, 자유 행에 $I$ 를 넣으면 된다.

R_np = np.array(R_sym.tolist(), dtype=float)
자유열_sym = [j for j in range(n) if j not in 피벗열_sym]

F = R_np[:len(피벗열_sym)][:, 자유열_sym]        # 피벗 행 x 자유 열
print(show_matrix(F, "F ="))

N_block = np.zeros((n, len(자유열_sym)))
N_block[피벗열_sym] = -F                          # 피벗 변수 자리에 -F
N_block[자유열_sym] = np.eye(len(자유열_sym))     # 자유 변수 자리에 I

print(show_matrix(N_block, "N = [-F ; I]  (원래 변수 순서)"))
print()
print("우리가 구한 특수해와 같은가 :", np.allclose(N_block, S))
print("A @ N 이 0 인가 :", np.allclose(A @ N_block, 0))

F =
[   2   -2 ]
[   0    2 ]
N = [-F ; I]  (원래 변수 순서)
[  -2    2 ]
[   1    0 ]
[  -0   -2 ]
[   0    1 ]

우리가 구한 특수해와 같은가 : True
A @ N 이 0 인가 : True

$[\,I \;\; F\,]$ 와 $\begin{bmatrix} -F \\ I \end{bmatrix}$ 를 곱하면 $-F + F$ 라서 0이 된다는 것도 직접 확인해 보자. 피벗 열을 앞으로 모은 순서로 보면 된다.

순서 = 피벗열_sym + 자유열_sym
R_정렬 = R_np[:, 순서]
N_정렬 = N_block[순서]

print(show_matrix(R_정렬, "R (피벗 열을 앞으로)"))
print(show_matrix(N_정렬, "N (같은 순서로)"))
print()
print(show_matrix(R_정렬 @ N_정렬, "R N ="))

R (피벗 열을 앞으로)
[   1    0    2   -2 ]
[   0    1    0    2 ]
[   0    0    0    0 ]
N (같은 순서로)
[  -2    2 ]
[  -0   -2 ]
[   1    0 ]
[   0    1 ]

R N =
[  0   0 ]
[  0   0 ]
[  0   0 ]

6. 랭크가 줄면 영공간이 자란다¶

$3 \times 3$ 행렬로 랭크를 하나씩 낮춰 가며 $\dim N(A) = n - r$ 을 확인해 보자.

예제들 = {
    "rank 3": np.array([[2.0, 1, 1], [1, 3, 2], [1, 0, 0]]),
    "rank 2": np.array([[2.0, 1, 1], [1, 3, 2], [3, 4, 3]]),   # row3 = row1 + row2
    "rank 1": np.array([[2.0, 1, 1], [4, 2, 2], [6, 3, 3]]),   # 모든 행이 row1 의 배수
}

print(f"{'':>8}{'rank r':>8}{'n - r':>8}{'dim N':>8}   모양")
for 이름, M in 예제들.items():
    r_i = np.linalg.matrix_rank(M)
    dimN = special_solutions(M).shape[1]
    모양 = {0: "원점", 1: "직선", 2: "평면", 3: "R^3 전체"}[dimN]
    print(f"{이름:>8}{r_i:>8}{3 - r_i:>8}{dimN:>8}   {모양}")

          rank r   n - r   dim N   모양
  rank 3       3       0       0   원점
  rank 2       2       1       1   직선
  rank 1       1       2       2   평면

랭크가 1 줄 때마다 영공간의 차원이 1 늘어난다. 랭크 2인 경우의 영공간은 직선이다.

M2 = 예제들["rank 2"]
d = special_solutions(M2)[:, 0]
d = d / np.linalg.norm(d)

t = np.array([-2.5, 2.5])
선 = t[:, None] * d

traces = [
    go.Scatter3d(x=선[:, 0], y=선[:, 1], z=선[:, 2], mode="lines",
                 line=dict(color=COLORS["nullspace"], width=7), name="N(A)"),
    go.Scatter3d(x=[0], y=[0], z=[0], mode="markers",
                 marker=dict(size=7, color="red"), name="원점"),
]
spin_figure(traces, title="rank 2 : 영공간은 직선", extent=2.6)

랭크 1인 경우에는 평면이 된다.

M1 = 예제들["rank 1"]
S1 = special_solutions(M1)
print(show_matrix(S1, "rank 1 행렬의 특수해 두 개"))
print("A @ S =", np.abs(M1 @ S1).max(), "(0 이어야 한다)")

traces = [
    plane(spans=[S1[:, 0], S1[:, 1]], extent=2.5,
          color=COLORS["nullspace"], opacity=0.4, name="N(A)"),
    go.Scatter3d(x=[0], y=[0], z=[0], mode="markers",
                 marker=dict(size=7, color="red"), name="원점"),
]
spin_figure(traces, title="rank 1 : 영공간은 평면", extent=2.6)

rank 1 행렬의 특수해 두 개
[  -0.5   -0.5 ]
[     1      0 ]
[     0      1 ]
A @ S = 0.0 (0 이어야 한다)

for j in range(S1.shape[1]):
    v = S1[:, j]
    print(f"특수해 {j + 1} = {np.round(v, 3)}   2x + y + z = {2*v[0] + v[1] + v[2]:+.2e}")

특수해 1 = [-0.5  1.   0. ]   2x + y + z = +0.00e+00
특수해 2 = [-0.5  0.   1. ]   2x + y + z = +0.00e+00

마치며...¶

서술 파트의 내용	이 노트북의 코드
피벗이 없으면 열을 건너뛴다	`row_echelon` 의 `continue`
랭크 $=$ 피벗의 개수	`len(피벗열)`
특수해 절차	`special_solutions` 의 역순 `for` 문
기저는 여러 가지, 공간은 하나	붙여서 랭크를 재는 방법
$N = [-F ; I]$	`N_block[피벗열] = -F`, `N_block[자유열] = I`
$\dim N(A) = n - r$	6절의 표

더 해 볼 것¶

row_echelon 이 L2처럼 절댓값이 가장 큰 성분을 피벗으로 고르도록 바꿔 보자. $U$ 가 달라지는가? 특수해는 달라지는가? 영공간은 달라지는가?
special_solutions 에서 자유 변수에 1 대신 2를 넣으면 무엇이 달라지는가? 여전히 영공간의 기저인가?
무작위 $m \times n$ 행렬을 여러 개 만들어 rank + dim N == n 이 늘 성립하는지 확인해 보자. 무작위 행렬의 랭크는 보통 얼마인가? 왜 그런가?
$A$ 와 $A^{\mathsf{T}}$ 의 랭크를 비교해 보자. 여러 행렬로 시험해 보면 무엇이 보이는가?

다음 강의에서는 우변이 0이 아닌 $A\vv{x} = \vv{b}$ 를 푼다. 이번 강의에서 구한 영공간이 그 답의 절반을 이룬다.