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

L8 서술 파트에서 완전해 $\vv{x} = \vv{x}_p + \vv{x}_n$ 과 랭크에 따른 네 경우를 다루었다. 이 노트북에서는 완전해를 구하는 함수를 짜고, 네 경우를 각각 만들어 해의 개수를 확인해 보도록 하자.

서술 파트의 내용	여기서 확인하는 방법
소거 후 $0 = (\text{0이 아닌 수})$ 이면 해가 없다	증강행렬의 마지막 열이 피벗 열이 되는지 본다
$\vv{y}^{\mathsf{T}}A = \vv{0}$ 이면 $\vv{y}^{\mathsf{T}}\vv{b} = 0$ 이어야 한다	왼쪽 영벡터를 구해 직접 확인한다
$\vv{x} = \vv{x}_p + \vv{x}_n$	`solve_general` 을 짜서 `lstsq` 와 대조한다
특수해를 바꿔도 해집합은 같다	다른 $\vv{x}_p$ 로 같은 집합이 나오는지 본다
랭크에 따른 네 경우	네 예제의 해의 개수를 세어 표로 만든다

0. 준비¶

L7에서 만든 사다리꼴 함수를 그대로 쓴다. 이번에는 증강행렬에도 적용할 것이다.

import numpy as np
import plotly.graph_objects as go

from linalg_viz import COLORS, arrow, plane, spin_figure, show_matrix

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

numpy 2.5.2

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

    for 열 in range(n):
        if 행 >= m:
            break
        후보 = 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, 피벗열

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

b_있음 = np.array([1.0, 5.0, 6.0])
b_없음 = np.array([1.0, 5.0, 7.0])

print(show_matrix(A, "A ="))
print("rank =", np.linalg.matrix_rank(A))

A =
[   1    2    2    2 ]
[   2    4    6    8 ]
[   3    6    8   10 ]
rank = 2

1. 해가 있는지 판정하기¶

증강행렬 $[A \mid \vv{b}]$ 를 소거해서 마지막 행을 보면 된다. 코드로는 더 깔끔한 판정법이 있다. 우변 열이 피벗 열이 되면 해가 없다. 그 열에서 피벗을 잡았다는 것은 곧 $0 = (\text{0이 아닌 수})$ 인 행이 생겼다는 뜻이기 때문이다.

n = A.shape[1]

for 이름, b in [("b = (1, 5, 6)", b_있음), ("b = (1, 5, 7)", b_없음)]:
    확장 = np.hstack([A, b.reshape(-1, 1)])
    U, 피벗열 = row_echelon(확장)
    print(f"{이름}")
    print(show_matrix(U, "  소거 결과 [U | c]"))
    print(f"  피벗 열 : {피벗열}")
    print(f"  우변 열({n})이 피벗 열인가 : {n in 피벗열}"
          f"  ->  {'해가 없다' if n in 피벗열 else '해가 있다'}")
    print()

b = (1, 5, 6)
  소거 결과 [U | c]
[  1   2   2   2   1 ]
[  0   0   2   4   3 ]
[  0   0   0   0   0 ]
  피벗 열 : [0, 2]
  우변 열(4)이 피벗 열인가 : False  ->  해가 있다

b = (1, 5, 7)
  소거 결과 [U | c]
[  1   2   2   2   1 ]
[  0   0   2   4   3 ]
[  0   0   0   0   1 ]
  피벗 열 : [0, 2, 4]
  우변 열(4)이 피벗 열인가 : True  ->  해가 없다

서술 파트에서 손으로 구한 것과 같다. $\vv{b} = (1,5,6)$ 이면 마지막 행이 $0 = 0$ 이고, $(1,5,7)$ 이면 $0 = 1$ 이다.

랭크로도 같은 판정을 할 수 있다. L6에서 쓴 방법이다.

for 이름, b in [("b = (1, 5, 6)", b_있음), ("b = (1, 5, 7)", b_없음)]:
    확장 = np.hstack([A, b.reshape(-1, 1)])
    print(f"{이름} :  rank(A) = {np.linalg.matrix_rank(A)},"
          f"  rank([A|b]) = {np.linalg.matrix_rank(확장)}")

b = (1, 5, 6) :  rank(A) = 2,  rank([A|b]) = 2
b = (1, 5, 7) :  rank(A) = 2,  rank([A|b]) = 3

2. 조건 $\vv{y}^{\mathsf{T}}\vv{b} = 0$ 확인하기¶

서술 파트에서 세 행 사이의 관계 $\vv{y} = (-1, -1, 1)$ 을 찾았다. 이것은 $\vv{y}^{\mathsf{T}}A = \vv{0}$ 을 만족하는 벡터이고, $A^{\mathsf{T}}$ 의 영공간에 들어 있다.

from scipy.linalg import null_space

y = null_space(A.T)[:, 0]
y = y / y[0] * -1          # 첫 성분이 -1 이 되도록 맞춘다

print("y =", y)
print("y^T A =", y @ A, " (영벡터)")
print()
for 이름, b in [("b = (1, 5, 6)", b_있음), ("b = (1, 5, 7)", b_없음)]:
    print(f"{이름} :  y . b = {y @ b:+.3f}"
          f"  ->  {'해가 있다' if abs(y @ b) < 1e-10 else '해가 없다'}")

y = [-1. -1.  1.]
y^T A = [0. 0. 0. 0.]  (영벡터)

b = (1, 5, 6) :  y . b = +0.000  ->  해가 있다
b = (1, 5, 7) :  y . b = +1.000  ->  해가 없다

소거로 판정한 결과와 같다. $\vv{y}$ 는 행들 사이의 관계를 담고 있고, 그 관계가 우변에도 그대로 성립해야 해가 존재한다.

풀어 쓰면 $-b_1 - b_2 + b_3 = 0$ , 즉 $b_3 = b_1 + b_2$ 이다.

for b in (b_있음, b_없음):
    print(f"b = {b} :  b1 + b2 = {b[0] + b[1]:g},  b3 = {b[2]:g}"
          f"  ->  {'조건 만족' if np.isclose(b[0] + b[1], b[2]) else '조건 위반'}")

b = [1. 5. 6.] :  b1 + b2 = 6,  b3 = 6  ->  조건 만족
b = [1. 5. 7.] :  b1 + b2 = 6,  b3 = 7  ->  조건 위반

3. 완전해 구하기¶

자유 변수를 전부 0으로 두고 후진 대입하면 특수해가 나온다. 영공간 기저는 L7의 방법으로 구한다. 두 가지를 합쳐 완전해를 돌려주는 함수를 만들자.

def special_solutions(A, tol=1e-10):
    """A x = 0 의 특수해들을 열로 쌓아 돌려준다. (L7과 같다)"""
    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)
    return np.column_stack(해들) if 해들 else np.zeros((n, 0))

def solve_general(A, b, tol=1e-10):
    """A x = b 의 완전해를 구한다.

    (해가 있는가, 특수해 x_p, 영공간 기저 N) 을 돌려준다.
    해가 없으면 (False, None, None).
    """
    A = np.asarray(A, dtype=float)
    b = np.asarray(b, dtype=float).reshape(-1)
    n = A.shape[1]

    U, 피벗열 = row_echelon(np.hstack([A, b.reshape(-1, 1)]), tol)
    if n in 피벗열:                       # 우변 열이 피벗 -> 0 = (0 아닌 수)
        return False, None, None

    # 자유 변수를 0으로 두고 피벗 변수를 뒤에서부터 확정한다
    x_p = np.zeros(n)
    for 행, 열 in reversed(list(enumerate(피벗열))):
        x_p[열] = (U[행, n] - U[행, 열 + 1:n] @ x_p[열 + 1:]) / U[행, 열]

    return True, x_p + 0.0, special_solutions(A, tol)

있음, x_p, N = solve_general(A, b_있음)

print("해가 있는가 :", 있음)
print("x_p =", x_p)
print(show_matrix(N, "영공간 기저 (열마다 하나)"))
print()
print("A @ x_p =", A @ x_p, "   b =", b_있음)
print("맞는가 :", np.allclose(A @ x_p, b_있음))

해가 있는가 : True
x_p = [-2.   0.   1.5  0. ]
영공간 기저 (열마다 하나)
[  -2    2 ]
[   1    0 ]
[   0   -2 ]
[   0    1 ]

A @ x_p = [1. 5. 6.]    b = [1. 5. 6.]
맞는가 : True

서술 파트에서 손으로 구한 $(-2, 0, 3/2, 0)$ 과 같다. 이제 영공간을 아무렇게나 더해 보자.

for _ in range(4):
    c = rng.standard_normal(N.shape[1])
    x = x_p + N @ c
    print(f"c = {np.round(c, 3)}  ->  A x = {A @ x}")

c = [ 0.126 -0.132]  ->  A x = [1. 5. 6.]
c = [0.64  0.105]  ->  A x = [1. 5. 6.]
c = [-0.536  0.362]  ->  A x = [1. 5. 6.]
c = [1.304 0.947]  ->  A x = [1. 5. 6.]

어떤 $c$ 를 넣어도 $A\vv{x}$ 가 $\vv{b}$ 그대로이다. 해가 없는 경우도 확인해 보자.

print("b = (1, 5, 7) 의 결과 :", solve_general(A, b_없음))

b = (1, 5, 7) 의 결과 : (False, None, None)

4. 특수해를 바꿔도 해집합은 같다¶

$\vv{x}_p$ 에 영공간의 아무 원소나 더하면 또 다른 특수해가 된다. 그것으로 만든 해집합이 원래와 같은지 확인해 보자. 두 집합이 같은지는 $\vv{x}_p' - \vv{x}_p$ 가 영공간에 있는지로 판정하면 된다.

x_p2 = x_p + 3.0 * N[:, 0] - 1.5 * N[:, 1]      # 다른 특수해

print("x_p  =", x_p)
print("x_p2 =", x_p2)
print()
print("A @ x_p2 =", A @ x_p2, " (b 와 같아야 한다)")
print("두 특수해의 차이 :", x_p2 - x_p)
print("그 차이가 영공간에 있는가 :", np.allclose(A @ (x_p2 - x_p), 0))
print("  -> 두 표현은 같은 해집합을 나타낸다")

x_p  = [-2.   0.   1.5  0. ]
x_p2 = [-11.    3.    4.5  -1.5]

A @ x_p2 = [1. 5. 6.]  (b 와 같아야 한다)
두 특수해의 차이 : [-9.   3.   3.  -1.5]
그 차이가 영공간에 있는가 : True
  -> 두 표현은 같은 해집합을 나타낸다

자유 변수를 0이 아닌 값으로 두고 구해도 마찬가지이다. 직접 만들어 보자.

# 자유 변수(2열, 4열)에 원하는 값을 주고 나머지를 맞춘다
자유값 = np.array([1.0, 2.0])
x_p3 = x_p + N @ 자유값

print("x_p3 =", x_p3)
print("자유 변수 자리 (x2, x4) :", x_p3[[1, 3]], " <- 우리가 준 값")
print("A @ x_p3 =", A @ x_p3)

x_p3 = [ 0.   1.  -2.5  2. ]
자유 변수 자리 (x2, x4) : [1. 2.]  <- 우리가 준 값
A @ x_p3 = [1. 5. 6.]

5. 랭크에 따른 네 경우¶

서술 파트의 표를 코드로 확인해 보자. 각 경우마다 예제를 하나씩 두고, 해가 있는 우변과 없는 우변을 시험한다.

def 보고(이름, A, b):
    """해의 존재와 개수를 한 줄로 요약한다."""
    m, n = A.shape
    r = np.linalg.matrix_rank(A)
    있음, x_p, N = solve_general(A, b)
    if not 있음:
        상태 = "해 없음"
    elif N.shape[1] == 0:
        상태 = "해 1개"
    else:
        상태 = f"해 무수히 많음 (자유도 {N.shape[1]})"
    print(f"  {이름:<22} m={m} n={n} r={r}   b={b}  ->  {상태}")

print("r = m = n  (정방, 가역)")
보고("[[1,2],[3,4]]", np.array([[1., 2], [3, 4]]), np.array([1.0, 2.0]))
보고("[[1,2],[3,4]]", np.array([[1., 2], [3, 4]]), np.array([7.0, -3.0]))
print()

print("r = n < m  (식이 미지수보다 많다)")
A2 = np.array([[1., 0], [0, 1], [1, 1]])
보고("[[1,0],[0,1],[1,1]]", A2, np.array([1.0, 2.0, 3.0]))
보고("[[1,0],[0,1],[1,1]]", A2, np.array([1.0, 2.0, 4.0]))
print()

print("r = m < n  (미지수가 식보다 많다)")
A3 = np.array([[1., 0, 1], [0, 1, 1]])
보고("[[1,0,1],[0,1,1]]", A3, np.array([1.0, 2.0]))
보고("[[1,0,1],[0,1,1]]", A3, np.array([-5.0, 9.0]))
print()

print("r < m,  r < n")
A4 = np.array([[1., 2, 3], [2, 4, 6]])
보고("[[1,2,3],[2,4,6]]", A4, np.array([1.0, 2.0]))
보고("[[1,2,3],[2,4,6]]", A4, np.array([1.0, 3.0]))

r = m = n  (정방, 가역)
  [[1,2],[3,4]]          m=2 n=2 r=2   b=[1. 2.]  ->  해 1개
  [[1,2],[3,4]]          m=2 n=2 r=2   b=[ 7. -3.]  ->  해 1개

r = n < m  (식이 미지수보다 많다)
  [[1,0],[0,1],[1,1]]    m=3 n=2 r=2   b=[1. 2. 3.]  ->  해 1개
  [[1,0],[0,1],[1,1]]    m=3 n=2 r=2   b=[1. 2. 4.]  ->  해 없음

r = m < n  (미지수가 식보다 많다)
  [[1,0,1],[0,1,1]]      m=2 n=3 r=2   b=[1. 2.]  ->  해 무수히 많음 (자유도 1)
  [[1,0,1],[0,1,1]]      m=2 n=3 r=2   b=[-5.  9.]  ->  해 무수히 많음 (자유도 1)

r < m,  r < n
  [[1,2,3],[2,4,6]]      m=2 n=3 r=1   b=[1. 2.]  ->  해 무수히 많음 (자유도 2)
  [[1,2,3],[2,4,6]]      m=2 n=3 r=1   b=[1. 3.]  ->  해 없음

정리하면 이렇다.

$r = m$ 인 경우(3번째)에는 어떤 우변을 넣어도 해가 있었다
$r = n$ 인 경우(1, 2번째)에는 해가 있을 때 딱 하나였다
둘 다 아닌 경우(4번째)에는 우변에 따라 없거나 무수히 많았다

식의 개수와 미지수의 개수만 보고는 알 수 없다. 2번째와 3번째를 비교해 보면 식이 더 많은 쪽이 해가 없을 수 있고, 미지수가 더 많은 쪽은 언제나 해가 있다.

6. 해집합은 영공간을 평행이동한 것¶

$3 \times 3$ 행렬로 그려 보자. 랭크가 1이면 영공간이 평면이고, 해집합도 나란한 평면이 된다.

A_r1 = np.array([[1., 2, 3],
                 [2., 4, 6],
                 [3., 6, 9]])          # 모든 행이 (1,2,3) 의 배수
b_r1 = np.array([1.0, 2.0, 3.0])       # = 첫 번째 열

있음, xp1, N1 = solve_general(A_r1, b_r1)
print("rank =", np.linalg.matrix_rank(A_r1), "  해가 있는가 :", 있음)
print("x_p =", xp1)
print(show_matrix(N1, "영공간 기저"))
print("A @ x_p =", A_r1 @ xp1)

rank = 1   해가 있는가 : True
x_p = [1. 0. 0.]
영공간 기저
[  -2   -3 ]
[   1    0 ]
[   0    1 ]
A @ x_p = [1. 2. 3.]

traces = [
    plane(spans=[N1[:, 0], N1[:, 1]], extent=1.1,
          color=COLORS["nullspace"], opacity=0.4, name="N(A)  (원점을 지난다)"),
    plane(spans=[N1[:, 0], N1[:, 1]], through=xp1, extent=1.1,
          color=COLORS["output"], opacity=0.35, name="해집합"),
    go.Scatter3d(x=[0], y=[0], z=[0], mode="markers",
                 marker=dict(size=7, color="red"), name="원점"),
]
traces += arrow([0, 0, 0], xp1, color=COLORS["input"], name="x_p")

spin_figure(traces, title="해집합은 영공간 평면을 x_p 만큼 옮긴 것", extent=1.4)

해집합 위의 점들이 모두 $x + 2y + 3z = 1$ 을 만족하는지도 확인해 보자.

for _ in range(4):
    c = rng.standard_normal(2)
    x = xp1 + N1 @ c
    print(f"x = {np.round(x, 3)}   x + 2y + 3z = {x[0] + 2*x[1] + 3*x[2]:.3f}"
          f"   A x = {A_r1 @ x}")

x = [ 6.204 -0.704 -1.265]   x + 2y + 3z = 1.000   A x = [1. 2. 3.]
x = [ 2.123 -0.623  0.041]   x + 2y + 3z = 1.000   A x = [1. 2. 3.]
x = [ 6.306 -2.325 -0.219]   x + 2y + 3z = 1.000   A x = [1. 2. 3.]
x = [ 5.689 -1.246 -0.732]   x + 2y + 3z = 1.000   A x = [1. 2. 3.]

마치며...¶

서술 파트의 내용	이 노트북의 코드
해의 존재 판정	우변 열이 피벗 열이 되는지 (`n in 피벗열`)
$\vv{y}^{\mathsf{T}}A = \vv{0}$ 조건	`null_space(A.T)` 로 $\vv{y}$ 를 구해 `y @ b`
특수해	`solve_general` 의 후진 대입
완전해	`x_p + N @ c`
특수해를 바꿔도 해집합은 같다	차이가 영공간에 있는지 확인
네 경우	5절의 `보고` 함수

더 해 볼 것¶

solve_general 이 np.linalg.lstsq 와 어떻게 다른지 확인해 보자. 해가 무수히 많을 때 lstsq 는 어떤 해를 고르는가? 그 해의 길이를 다른 해들과 비교해 보자.
해가 없는 경우에 lstsq 는 무엇을 돌려주는가? 그 결과에 $A$ 를 곱하면 $\vv{b}$ 와 얼마나 다른가?
5절의 네 경우에서 우변을 무작위로 만들어 여러 번 시험해 보자. $r = m$ 인 경우에는 정말 한 번도 실패하지 않는가?
$\vv{y}^{\mathsf{T}}A = \vv{0}$ 인 $\vv{y}$ 들이 이루는 집합도 부분공간이다. 그 차원은 $m - r$ 인가? 여러 행렬로 확인해 보자.

다음 강의에서는 지금까지 정의 없이 써 온 독립, 기저, 차원을 정확히 정의한다.