서술에서 두 가지를 만들었다. 없던 미지수와, 답을 고르는 기준이다. 이 노트북은 그 둘을 흔든다.
검산 — 스물둘의 답을 전부 다시 잰다.
기준 바꾸기 — 같은 제약에 다른 기준을 걸면 답이 어디로 가는가.
부수기 — 전치를 빠뜨리면 회전이 반대로 나오고, 자유재를 기준으로 잡으면 계가 통째로 무너진다.
import itertools
import numpy as np
np.set_printoptions(precision=4, suppress=True, linewidth=115)
def 같다(번호, 얻은, 적힌, 말="", 허용=1e-9):
얻은, 적힌 = np.asarray(얻은, dtype=float), np.asarray(적힌, dtype=float)
차 = np.abs(얻은 - 적힌).max()
print(f"[{번호}] {'맞다 ' if 차 < 허용 else '틀리다'} "
f"차이 {차:.2e} {말}")
def 풀린다(번호, A, b, 답, 말=""):
A, b = np.asarray(A, dtype=float), np.asarray(b, dtype=float)
x = np.linalg.solve(A, b)
같다(번호, x, 답, 말)
return xA = np.array([[4.0, 1], [2, 3]])
행최소최대 = A.min(axis=1).max()
열최대최소 = A.max(axis=0).min()
print(f"행 최소의 최대 {행최소최대:.0f}, 열 최대의 최소 {열최대최소:.0f}"
f" -> {'안장점 있음' if 행최소최대 == 열최대최소 else '안장점 없다. 섞어야 한다'}")
K = np.array([[4.0, 2, -1], [1, 3, -1], [1, 1, 0]])
s = 풀린다("23", K, [0, 0, 1], [0.25, 0.75, 2.5], "세관 (1/4, 3/4), 게임의 값 5/2")
print(f" det K = {np.linalg.det(K):.0f}")
x = s[:2]
같다("23", A.T @ x, [2.5, 2.5], "밀수업자가 어디로 와도 같은 값")
y = np.linalg.solve(np.array([[4.0, 1, -1], [2, 3, -1], [1, 1, 0]]), [0, 0, 1])
같다("23", y[:2], [0.5, 0.5], "밀수업자 쪽은 (1/2, 1/2)")
같다("23", y[2], s[2], "두 사람이 다른 계를 푸는데 값은 같다")행 최소의 최대 2, 열 최대의 최소 3 -> 안장점 없다. 섞어야 한다
[23] 맞다 차이 0.00e+00 세관 (1/4, 3/4), 게임의 값 5/2
det K = 4
[23] 맞다 차이 0.00e+00 밀수업자가 어디로 와도 같은 값
[23] 맞다 차이 0.00e+00 밀수업자 쪽은 (1/2, 1/2)
[23] 맞다 차이 0.00e+00 두 사람이 다른 계를 푸는데 값은 같다
문제 26에서 여유변수를 넣었다. 기저해 하나가 꼭짓점 하나라는 것을 전수 조사로 확인한다.
Am = np.array([[1.0, 0, 1, 0, 0], [0, 2, 0, 1, 0], [3, 2, 0, 0, 1]])
b = np.array([4.0, 12, 18])
c = np.array([3.0, 5, 0, 0, 0])
기저해, 종속 = [], 0
for S in itertools.combinations(range(5), 3):
B = Am[:, S]
if abs(np.linalg.det(B)) < 1e-9:
종속 += 1
continue
x = np.zeros(5)
x[list(S)] = np.linalg.solve(B, b)
기저해.append((S, x, float(c @ x), float(np.linalg.det(B))))
print(f"C(5,3) = {len(list(itertools.combinations(range(5), 3)))}, "
f"종속 {종속}개, 기저해 {len(기저해)}개")
가능 = [t for t in 기저해 if (t[1] >= -1e-9).all()]
print(f"그중 성분이 전부 0 이상인 것 {len(가능)}개\n")
print(f"{'고른 열':>14} {'(x1, x2)':>14} {'z':>6} {'det B':>7} 판정")
for S, x, z, d in sorted(기저해, key=lambda t: t[2]):
좋 = (x >= -1e-9).all()
print(f"{str(S):>14} {str(np.round(x[:2], 2)):>14} {z:>6.0f} {d:>7.1f} "
f"{'꼭짓점' if 좋 else '음수가 있다'}"
f"{' <- 최적' if 좋 and abs(z - 36) < 1e-9 else ''}")C(5,3) = 10, 종속 2개, 기저해 8개
그중 성분이 전부 0 이상인 것 5개
고른 열 (x1, x2) z det B 판정
(2, 3, 4) [0. 0.] 0 1.0 꼭짓점
(0, 3, 4) [4. 0.] 12 1.0 꼭짓점
(0, 2, 3) [6. 0.] 18 3.0 음수가 있다
(0, 1, 3) [4. 3.] 27 -2.0 꼭짓점
(1, 2, 4) [0. 6.] 30 -2.0 꼭짓점
(0, 1, 2) [2. 6.] 36 -6.0 꼭짓점 <- 최적
(0, 1, 4) [4. 6.] 42 2.0 음수가 있다
(1, 2, 3) [0. 9.] 45 2.0 음수가 있다
문제 27의 축 교환. 저장하는 것이 역행렬 아홉 칸이 아니라 소거행렬마다 바뀐 열 하나씩이라는 것을 확인한다.
E1 = np.array([[1.0, 0, 0], [0, 0.5, 0], [0, -1, 1]])
E2 = np.array([[1.0, 0, -1/3], [0, 1, 0], [0, 0, 1/3]])
같다("25", E1 @ Am[:, 1], [0, 1, 0], "E1 이 제품 2 의 열을 e2 로 만든다")
Binv = E2 @ E1
B = Am[:, [2, 1, 0]] # 기저는 (s1, x2, x1)
같다("25", Binv, np.linalg.inv(B), "E2 E1 == B^-1")
xB = 풀린다("25", B, b, [2, 6, 2], "s1 = 2, x2 = 6, x1 = 2")
print(f" B^-1 =\n{np.round(Binv, 4)}")
# E 하나가 정말 한 열만 다른가
for 이름, E in [("E1", E1), ("E2", E2)]:
다른열 = [j for j in range(3) if not np.allclose(E[:, j], np.eye(3)[:, j])]
print(f" {이름} 은 단위행렬과 {len(다른열)}개 열만 다르다 (열 {다른열})")[25] 맞다 차이 0.00e+00 E1 이 제품 2 의 열을 e2 로 만든다
[25] 맞다 차이 0.00e+00 E2 E1 == B^-1
[25] 맞다 차이 0.00e+00 s1 = 2, x2 = 6, x1 = 2
B^-1 =
[[ 1. 0.3333 -0.3333]
[ 0. 0.5 0. ]
[ 0. -0.3333 0.3333]]
E1 은 단위행렬과 1개 열만 다르다 (열 [1])
E2 은 단위행렬과 1개 열만 다르다 (열 [2])
문제 28의 잠재가격. 다시 풀지 않고도 잔업 값을 알 수 있다는 것을 실제로 다시 풀어 확인한다.
cB = np.array([0.0, 5, 3])
y = 풀린다("26", B.T, cB, [0, 1.5, 1], "공장별 시간의 값")
기본이익 = float(c[[2, 1, 0]] @ (Binv @ b))
print(f"\n{'가용시간':>22} {'다시 푼 이익':>12} {'y 로 예측':>12}")
for k, 이름 in enumerate(["1공장", "2공장", "3공장"]):
b2 = b.copy()
b2[k] += 1
이익 = float(c[[2, 1, 0]] @ (Binv @ b2))
print(f"{이름 + ' +1시간':>22} {이익:>12.2f} {기본이익 + y[k]:>12.2f}")
같다("26", cB @ Binv, y, "z = y^T b 라는 관계")[26] 맞다 차이 0.00e+00 공장별 시간의 값
가용시간 다시 푼 이익 y 로 예측
1공장 +1시간 36.00 36.00
2공장 +1시간 37.50 37.50
3공장 +1시간 37.00 37.00
[26] 맞다 차이 0.00e+00 z = y^T b 라는 관계
문제 29의 네트워크. 경로를 다 세지 않고 최선임을 확인한다.
Ai = np.array([[-1.0, 1, 0, 0], [0, -1, 1, 0], [-1, 0, 1, 0],
[0, 0, -1, 1], [0, -1, 0, 1]])
cc = np.array([2.0, 1, 4, 3, 5])
f = np.array([-6.0, 0, 0, 6])
print(f"rank A = {np.linalg.matrix_rank(Ai)} (노드 4개니까 3이 맞다)")
print(f"루프 차원 = {5 - np.linalg.matrix_rank(Ai)} = m - n + 1 = {5 - 4 + 1}")
yf = np.array([6.0, 6, 0, 6, 0])
같다("27", Ai.T @ yf, f, "흐름 보존")
print(f" 비용 {cc @ yf:.0f}만원")
x전 = np.array([0.0, 2, 3, 6])
축소 = cc - Ai @ x전
print(f" 축소비용 {축소} -> 전부 0 이상이면 최선")
# 경로를 직접 세어 비교
경로 = {"공장-거점1-거점2-창고": [0, 1, 3], "공장-거점2-창고": [2, 3],
"공장-거점1-창고": [0, 4]}
for 이름, 간선들 in 경로.items():
print(f" {이름:>22} : {cc[간선들].sum():.0f}만원")rank A = 3 (노드 4개니까 3이 맞다)
루프 차원 = 2 = m - n + 1 = 2
[27] 맞다 차이 0.00e+00 흐름 보존
비용 36만원
축소비용 [0. 0. 1. 0. 1.] -> 전부 0 이상이면 최선
공장-거점1-거점2-창고 : 6만원
공장-거점2-창고 : 7만원
공장-거점1-창고 : 7만원
T = np.array([[1.0, 1, 1, 0, 0, 0], [0, 0, 0, 1, 1, 1],
[1, 0, 0, 1, 0, 0], [0, 1, 0, 0, 1, 0], [0, 0, 1, 0, 0, 1]])
print(f"rank T = {np.linalg.matrix_rank(T)} = m + n - 1 = {2 + 3 - 1}")
같다("28", T[0] + T[1], T[2] + T[3] + T[4], "공급 줄의 합 == 수요 줄의 합")
x = np.array([20.0, 10, 0, 0, 5, 15])
같다("28", T @ x, [30, 20, 20, 15, 15], "북서코너 해가 제약을 만족")
print(f" 0 아닌 성분 {int((x > 0).sum())}개 = m + n - 1")
print()
Z0 = np.array([[10.0, 20], [30, 40]])
r, s = np.array([2.0, 1]), np.array([1.0, 0.5])
Z = np.diag(r) @ Z0 @ np.diag(s)
같다("31", Z, [[20, 20], [30, 20]], "RAS 로 만든 새 표")
같다("31", Z.sum(axis=1), [40, 50], "행 합")
같다("31", Z.sum(axis=0), [50, 40], "열 합")
같다("31", np.diag([4.0, 2]) @ Z0 @ np.diag([0.5, 0.25]), Z,
"배율을 (4,2),(1/2,1/4) 로 바꿔도 같은 표")rank T = 4 = m + n - 1 = 4
[28] 맞다 차이 0.00e+00 공급 줄의 합 == 수요 줄의 합
[28] 맞다 차이 0.00e+00 북서코너 해가 제약을 만족
0 아닌 성분 4개 = m + n - 1
[31] 맞다 차이 0.00e+00 RAS 로 만든 새 표
[31] 맞다 차이 0.00e+00 행 합
[31] 맞다 차이 0.00e+00 열 합
[31] 맞다 차이 0.00e+00 배율을 (4,2),(1/2,1/4) 로 바꿔도 같은 표
# RAS 반복이 정말 수렴하는가
Zr = Z0.copy()
목표행, 목표열 = np.array([40.0, 50]), np.array([50.0, 40])
print(f"{'회':>3} {'표':>34} {'참값과의 차이':>14}")
for it in range(1, 6):
Zr = Zr * (목표행 / Zr.sum(axis=1))[:, None]
Zr = Zr * (목표열 / Zr.sum(axis=0))[None, :]
print(f"{it:>3} {str(np.round(Zr, 4).tolist()):>34} "
f"{np.abs(Zr - Z).max():>14.2e}")
같다("31", Zr, Z, "RAS 5회면 수렴", 허용=1e-6) 회 표 참값과의 차이
1 [[19.1781, 19.3103], [30.8219, 20.6897]] 8.22e-01
2 [[19.9918, 19.9931], [30.0082, 20.0069]] 8.25e-03
3 [[19.9999, 19.9999], [30.0001, 20.0001]] 8.25e-05
4 [[20.0, 20.0], [30.0, 20.0]] 8.25e-07
5 [[20.0, 20.0], [30.0, 20.0]] 8.25e-09
[31] 맞다 차이 8.25e-09 RAS 5회면 수렴
v = np.array([[12.0, 5, 3], [10, 8, 4], [9, 6, 7]])
최선 = max((sum(v[i, p[i]] for i in range(3)), p)
for p in itertools.permutations(range(3)))
print(f"최적 배정 {최선[1]} 값 {최선[0]:.0f}")
print(" 다른 배정들 :",
sorted(sum(v[i, p[i]] for i in range(3))
for p in itertools.permutations(range(3))))
u, p = np.array([10.0, 8, 7]), np.array([2.0, 0, 0])
등식쌍 = [(i, j) for i in range(3) for j in range(3)
if abs(u[i] + p[j] - v[i, j]) < 1e-9]
print(f" 등식이 걸리는 쌍 {len(등식쌍)}개 : {등식쌍}")
여유 = min(u[i] + p[j] - v[i, j] for i in range(3) for j in range(3))
print(f" 모든 쌍에서 여유 >= 0 인가 : {여유 >= -1e-9} (최소 {여유:.1f})")
Ae = np.zeros((len(등식쌍), 6))
for r_, (i, j) in enumerate(등식쌍):
Ae[r_, i] = 1
Ae[r_, 3 + j] = 1
print(f" rank A = {np.linalg.matrix_rank(Ae)}, "
f"영공간 차원 {6 - np.linalg.matrix_rank(Ae)}")
같다("29", Ae @ np.array([1.0, 1, 1, -1, -1, -1]), np.zeros(len(등식쌍)),
"(1,1,1|-1,-1,-1) 이 영공간에 있다")
print()
print(" t 를 키우면 (이득을 회사가 더 가져가면) :")
for t in [0.0, 1.0, 2.0, 2.5]:
uu, pp = u - t, p + t
표 = "" if (pp >= -1e-9).all() else " <- 월세가 음수. 불가능"
print(f" t={t:.1f} -> u={uu}, p={pp}{표}")최적 배정 (0, 1, 2) 값 27
다른 배정들 : [np.float64(18.0), np.float64(19.0), np.float64(20.0), np.float64(22.0), np.float64(22.0), np.float64(27.0)]
등식이 걸리는 쌍 5개 : [(0, 0), (1, 0), (1, 1), (2, 0), (2, 2)]
모든 쌍에서 여유 >= 0 인가 : True (최소 0.0)
rank A = 5, 영공간 차원 1
[29] 맞다 차이 0.00e+00 (1,1,1|-1,-1,-1) 이 영공간에 있다
t 를 키우면 (이득을 회사가 더 가져가면) :
t=0.0 -> u=[10. 8. 7.], p=[2. 0. 0.]
t=1.0 -> u=[9. 7. 6.], p=[3. 1. 1.]
t=2.0 -> u=[8. 6. 5.], p=[4. 2. 2.]
t=2.5 -> u=[7.5 5.5 4.5], p=[4.5 2.5 2.5]
def 합사상(n):
M = np.zeros((2 * n, n * n))
for i in range(n):
for j in range(n):
M[i, i * n + j] = 1
M[n + j, i * n + j] = 1
return M
print(f"{'표 크기':>8} {'미지수':>7} {'줄':>4} {'rank':>5} {'영공간':>7} "
f"{'(m-1)(n-1)':>11}")
for n in (2, 3, 4, 5):
M = 합사상(n)
r_ = np.linalg.matrix_rank(M)
print(f"{f'{n}x{n}':>8} {n*n:>7} {2*n:>4} {r_:>5} {n*n - r_:>7} "
f"{(n-1)**2:>11}")
Za = np.array([[10.0, 20, 30], [20, 30, 10], [30, 10, 20]])
Zb = np.array([[15.0, 15, 30], [15, 35, 10], [30, 10, 20]])
같다("32", Za.sum(axis=0), Zb.sum(axis=0), "열 합이 같다")
같다("32", Za.sum(axis=1), Zb.sum(axis=1), "행 합도 같다")
print(f"[32] 그런데 아홉 칸 중 {int((Za != Zb).sum())}개가 다르다") 표 크기 미지수 줄 rank 영공간 (m-1)(n-1)
2x2 4 4 3 1 1
3x3 9 6 5 4 4
4x4 16 8 7 9 9
5x5 25 10 9 16 16
[32] 맞다 차이 0.00e+00 열 합이 같다
[32] 맞다 차이 0.00e+00 행 합도 같다
[32] 그런데 아홉 칸 중 4개가 다르다
S = np.array([[2.0, 1, 0], [1, 2, 1], [0, 1, 2]])
하나 = np.ones(3)
print(f"S 의 고윳값 {np.linalg.eigvalsh(S)} -> 양의 정부호")
Si = np.linalg.inv(S)
w = Si @ 하나 / (하나 @ Si @ 하나)
같다("33", w, [0.5, 0, 0.5], "최소분산 비중")
같다("33", w @ S @ w, 1.0, "그때 분산")
K = np.block([[2 * S, 하나.reshape(3, 1)], [하나.reshape(1, 3), np.zeros((1, 1))]])
sol = np.linalg.solve(K, [0, 0, 0, 1])
같다("33", sol[:3], w, "테두리 계로 풀어도 같다")
print(f" 승수 lambda = {sol[3]:.4f}")
print()
print("같은 제약 1^T w = 1 에 다른 기준을 걸면 :")
후보 = [
("분산 최소", w),
("균등 배분", 하나 / 3),
("첫 자산에 몰빵", np.array([1.0, 0, 0])),
("가운데에 몰빵", np.array([0.0, 1, 0])),
("가운데를 공매도 (1,-1,1)", np.array([1.0, -1, 1])),
]
print(f"{'기준':>24} {'w':>22} {'분산':>8}")
for 이름, ww in 후보:
print(f"{이름:>24} {str(np.round(ww, 3)):>22} {ww @ S @ ww:>8.4f}")S 의 고윳값 [0.5858 2. 3.4142] -> 양의 정부호
[33] 맞다 차이 1.11e-16 최소분산 비중
[33] 맞다 차이 0.00e+00 그때 분산
[33] 맞다 차이 1.11e-16 테두리 계로 풀어도 같다
승수 lambda = -2.0000
같은 제약 1^T w = 1 에 다른 기준을 걸면 :
기준 w 분산
분산 최소 [0.5 0. 0.5] 1.0000
균등 배분 [0.333 0.333 0.333] 1.1111
첫 자산에 몰빵 [1. 0. 0.] 2.0000
가운데에 몰빵 [0. 1. 0.] 2.0000
가운데를 공매도 (1,-1,1) [ 1. -1. 1.] 2.0000
분산 최소가 정말 최소인가 — 제약을 만족하는 점을 무작위로 잔뜩 뽑아 확인한다.
난수 = np.random.default_rng(39)
표본 = 난수.normal(size=(200000, 3))
표본 = 표본 / 표본.sum(axis=1, keepdims=True) # 1^T w = 1 로 맞춘다
분산들 = np.einsum("ij,jk,ik->i", 표본, S, 표본)
print(f"무작위 20만 개 중 최소 분산 {분산들.min():.6f}")
print(f"최적해의 분산 {w @ S @ w:.6f}")
print(f"최적보다 작은 표본 개수 {(분산들 < w @ S @ w - 1e-12).sum()}개")무작위 20만 개 중 최소 분산 1.000020
최적해의 분산 1.000000
최적보다 작은 표본 개수 0개
문제 36의 KKT. 한계비용이 같아진다는 것이 답의 뜻이었다.
Q = np.diag([2.0, 4, 4])
Aq = np.ones((1, 3))
KK = np.block([[Q, Aq.T], [Aq, np.zeros((1, 1))]])
sol = np.linalg.solve(KK, [0, 0, 0, 12])
x, lam = sol[:3], sol[3]
같다("34", x, [6, 3, 3], "라인별 처리량")
print(f" 요금 {0.5 * x @ Q @ x:.0f}만원, 승수 {lam:+.0f}, det {np.linalg.det(KK):.0f}")
같다("34", Q @ x, -lam * np.ones(3), "한계비용이 세 라인에서 같다")
print()
print(f"{'배분':>22} {'요금':>8}")
for 이름, xx in [("최적 (6,3,3)", x), ("균등 (4,4,4)", np.full(3, 4.0)),
("몰빵 (12,0,0)", np.array([12.0, 0, 0])),
("반장 절충 (8,2,2)", np.array([8.0, 2, 2]))]:
print(f"{이름:>22} {0.5 * xx @ Q @ xx:>8.1f}")
print()
print("총량을 12에서 늘리면 요금이 얼마씩 느는가 (승수의 뜻) :")
for 총 in (12.0, 13.0, 14.0):
xx = np.linalg.solve(KK, [0, 0, 0, 총])[:3]
print(f" 총량 {총:.0f} -> 요금 {0.5 * xx @ Q @ xx:>7.2f}")
print(f" 승수의 크기는 {abs(lam):.0f} 이지만 실제 차분은 "
f"{0.5*13**2 - 0.5*12**2:.1f} 이다. 승수는 미분이지 차분이 아니다.")[34] 맞다 차이 0.00e+00 라인별 처리량
요금 72만원, 승수 -12, det -32
[34] 맞다 차이 0.00e+00 한계비용이 세 라인에서 같다
배분 요금
최적 (6,3,3) 72.0
균등 (4,4,4) 80.0
몰빵 (12,0,0) 144.0
반장 절충 (8,2,2) 80.0
총량을 12에서 늘리면 요금이 얼마씩 느는가 (승수의 뜻) :
총량 12 -> 요금 72.00
총량 13 -> 요금 84.50
총량 14 -> 요금 98.00
승수의 크기는 12 이지만 실제 차분은 12.5 이다. 승수는 미분이지 차분이 아니다.
문제 37에서 자체는 부정부호인데 제약 위에서는 음의 정부호였다.
H = np.diag([1.0, -3, -3])
Z = np.array([[1.0, 1], [-1, 0], [0, -1]])
print(f"H 의 고윳값 {np.linalg.eigvalsh(H)} -> 부정부호")
ZHZ = Z.T @ H @ Z
같다("35", ZHZ, [[-2, 1], [1, -2]], "축소 헤시안")
print(f"Z^T H Z 의 고윳값 {np.linalg.eigvalsh(ZHZ)} -> 음의 정부호")
print()
print("예산을 유지하는 방향으로 움직여 보면 :")
for 이름, vv in [("커피+t, 빵-t", np.array([1.0, -1, 0])),
("커피+t, 우유-t", np.array([1.0, 0, -1])),
("커피+t, 빵/우유 -t/2", np.array([1.0, -0.5, -0.5])),
("빵+t, 우유-t", np.array([0.0, 1, -1]))]:
assert abs(np.ones(3) @ vv) < 1e-12
print(f" {이름:>22} : v^T H v = {vv @ H @ vv:+7.3f} "
f"-> 만족 변화 {0.5 * (vv @ H @ vv):+.3f} t^2")
print(" 전부 음수다. 어느 방향으로 가도 나빠진다 -> 최대점이 맞다")H 의 고윳값 [-3. -3. 1.] -> 부정부호
[35] 맞다 차이 0.00e+00 축소 헤시안
Z^T H Z 의 고윳값 [-3. -1.] -> 음의 정부호
예산을 유지하는 방향으로 움직여 보면 :
커피+t, 빵-t : v^T H v = -2.000 -> 만족 변화 -1.000 t^2
커피+t, 우유-t : v^T H v = -2.000 -> 만족 변화 -1.000 t^2
커피+t, 빵/우유 -t/2 : v^T H v = -0.500 -> 만족 변화 -0.250 t^2
빵+t, 우유-t : v^T H v = -6.000 -> 만족 변화 -3.000 t^2
전부 음수다. 어느 방향으로 가도 나빠진다 -> 최대점이 맞다
문제 38의 최소노름. 다른 기준을 골랐다면 다른 답이 나온다.
A = np.array([[0.5, 0.1, 0.3], [0.2, 0.5, 0.1], [0.1, 0.2, 0.4]])
L = np.linalg.inv(np.eye(3) - A)
같다("36", L, [[2.8, 1.2, 1.6], [1.3, 2.7, 1.1], [0.9, 1.1, 2.3]], "레온티예프 역행렬")
a = np.array([10.0, 0, 0]) @ L
같다("36", a, [28, 12, 16], "유발 배출계수")
d = a * 296 / (a @ a)
같다("36", d, [7, 3, 4], "최소노름 해")
같다("36", a @ d, 296.0, "목표를 정확히 맞춘다")
print(f" 노름 제곱 {d @ d:.0f}, 총산출 {np.round(L @ d, 2)}, "
f"배출 {10 * (L @ d)[0]:.1f}")
print()
print("다른 기준을 골랐다면 :")
기준 = [("노름 최소", d)]
# 기존 수요 (100,50,50) 에서 가장 적게 바뀌게
기존 = np.array([100.0, 50, 50])
기준.append(("기존에서 최소 변경", 기존 + a * (296 - a @ 기존) / (a @ a)))
# 농림만 조정
기준.append(("농림만 조정", np.array([296 / 28, 0, 0])))
# 균등
기준.append(("균등 배분", np.full(3, 296 / a.sum())))
print(f"{'기준':>22} {'d':>26} {'노름^2':>10} {'배출':>7}")
for 이름, dd in 기준:
print(f"{이름:>22} {str(np.round(dd, 3)):>26} {dd @ dd:>10.2f} {a @ dd:>7.1f}")
print(" 넷 다 296 을 맞춘다. '옳은' 답은 없고 고른 기준만 있다.")[36] 맞다 차이 4.44e-16 레온티예프 역행렬
[36] 맞다 차이 3.55e-15 유발 배출계수
[36] 맞다 차이 8.88e-16 최소노름 해
[36] 맞다 차이 0.00e+00 목표를 정확히 맞춘다
노름 제곱 74, 총산출 [29.6 21.6 18.8], 배출 296.0
다른 기준을 골랐다면 :
기준 d 노름^2 배출
노름 최소 [7. 3. 4.] 74.00 296.0
기존에서 최소 변경 [ 7.676 10.432 -2.757] 175.35 296.0
농림만 조정 [10.571 0. 0. ] 111.76 296.0
균등 배분 [5.286 5.286 5.286] 83.82 296.0
넷 다 296 을 맞춘다. '옳은' 답은 없고 고른 기준만 있다.
문제 39의 D-최적 설계. 행렬식이 부피이고 부피가 신뢰도였다.
설계 = {
"선택 1 (직교)": np.array([[1.0, 1, -1, -1], [1, -1, 1, -1],
[1, -1, -1, 1], [1, 1, 1, 1]]),
"선택 2 (하나씩)": np.array([[1.0, -1, -1, -1], [1, 1, -1, -1],
[1, -1, 1, -1], [1, -1, -1, 1]]),
"선택 3 (C 고정)": np.array([[1.0, 1, -1, 1], [1, -1, 1, 1],
[1, -1, -1, 1], [1, 1, 1, 1]]),
}
print(f"{'설계':>16} {'|det X|':>9} {'det(X^TX)':>11} {'X^TX 대각':>10} "
f"{'계수 분산의 합':>14}")
for 이름, X in 설계.items():
d_ = abs(np.linalg.det(X))
G = X.T @ X
대각 = np.allclose(G, np.diag(np.diag(G)))
if abs(np.linalg.det(G)) < 1e-9:
분산 = float("inf")
else:
분산 = float(np.trace(np.linalg.inv(G)))
print(f"{이름:>16} {d_:>9.2f} {np.linalg.det(G):>11.1f} {str(대각):>10} "
f"{분산:>14.4f}")
print(f" 4x4 에서 성분이 +-1 일 때 행렬식 상한(하다마르)은 {4 ** 2} 이다") 설계 |det X| det(X^TX) X^TX 대각 계수 분산의 합
선택 1 (직교) 16.00 256.0 True 1.0000
선택 2 (하나씩) 8.00 64.0 False 2.5000
선택 3 (C 고정) 0.00 0.0 False inf
4x4 에서 성분이 +-1 일 때 행렬식 상한(하다마르)은 16 이다
문제 40의 대칭 사영. 버린 조각이 남긴 것과 직교해야 한다.
M = np.array([[-0.4, 0.1, 0.3], [0.1, -0.4, 0.1], [-0.1, 0.1, -0.4]])
Sy, Kw = (M + M.T) / 2, (M - M.T) / 2
같다("38", Sy, [[-0.4, 0.1, 0.1], [0.1, -0.4, 0.1], [0.1, 0.1, -0.4]], "남긴 대칭 부분")
같다("38", np.trace(Sy.T @ Kw), 0.0, "두 조각이 직교")
같다("38", (M ** 2).sum(), (Sy ** 2).sum() + (Kw ** 2).sum(), "피타고라스")
print(f" ||M||^2 = {(M**2).sum():.2f} = {(Sy**2).sum():.2f} + {(Kw**2).sum():.2f}")
print(f" 남긴 것의 고윳값 {np.linalg.eigvalsh(Sy)} -> 음의 정부호")
print()
print("다른 대칭행렬로 고쳤다면 더 멀어지는가 :")
난 = np.random.default_rng(38)
더가까움 = 0
for _ in range(20000):
잡 = 난.normal(scale=0.05, size=(3, 3))
후보 = Sy + (잡 + 잡.T) / 2
if ((M - 후보) ** 2).sum() < ((M - Sy) ** 2).sum() - 1e-12:
더가까움 += 1
print(f" 무작위 대칭행렬 2만 개 중 더 가까운 것 {더가까움}개")[38] 맞다 차이 1.39e-17 남긴 대칭 부분
[38] 맞다 차이 0.00e+00 두 조각이 직교
[38] 맞다 차이 0.00e+00 피타고라스
||M||^2 = 0.62 = 0.54 + 0.08
남긴 것의 고윳값 [-0.5 -0.5 -0.2] -> 음의 정부호
다른 대칭행렬로 고쳤다면 더 멀어지는가 :
무작위 대칭행렬 2만 개 중 더 가까운 것 0개
x1 = np.array([-3.0, -1, 1, 3])
x2 = np.array([-1.0, -1, 1, 1])
y = np.array([-6.0, -4, 2, 8])
X = np.c_[x1, x2]
be = np.linalg.solve(X.T @ X, X.T @ y)
같다("39", be, [2, 1], "다중회귀 계수")
단순 = x1 @ y / (x1 @ x1)
같다("39", 단순, 2.4, "가격을 무시한 단순회귀")
r1 = x1 - (x1 @ x2 / (x2 @ x2)) * x2
같다("39", r1, [-1, 1, -1, 1], "가격으로 설명되지 않는 광고비")
같다("39", r1 @ y / (r1 @ r1), be[0], "FWL : 두 단계가 한 번에 푼 것과 같다")
델타 = x1 @ x2 / (x1 @ x1)
같다("39", 단순, be[0] + be[1] * 델타, "누락변수 편의 = (빠진 계수) x (두 변수의 관계)")
print(f" 2.4 = 2.0 + 1.0 x {델타:.1f}")[39] 맞다 차이 3.33e-15 다중회귀 계수
[39] 맞다 차이 0.00e+00 가격을 무시한 단순회귀
[39] 맞다 차이 0.00e+00 가격으로 설명되지 않는 광고비
[39] 맞다 차이 1.33e-15 FWL : 두 단계가 한 번에 푼 것과 같다
[39] 맞다 차이 0.00e+00 누락변수 편의 = (빠진 계수) x (두 변수의 관계)
2.4 = 2.0 + 1.0 x 0.4
xs = np.array([1.0, 2, 3, 5, 6, 7])
ys = np.array([10.0, 11, 12, 4, 5, 6])
xc, yc = xs - xs.mean(), ys - ys.mean()
통합 = xc @ yc / (xc @ xc)
같다("40", 통합, -8/7, "여섯 점을 한꺼번에 놓은 기울기")
D = np.zeros((6, 2))
D[:3, 0] = 1
D[3:, 1] = 1
MD = np.eye(6) - D @ np.linalg.inv(D.T @ D) @ D.T
print(f" rank M_D = {np.linalg.matrix_rank(MD)} = 6 - 2")
xw, yw = MD @ xs, MD @ ys
안 = xw @ yw / (xw @ xw)
같다("40", 안, 1.0, "공장 안에서의 기울기")
절편 = [ys[:3].mean() - 안 * xs[:3].mean(), ys[3:].mean() - 안 * xs[3:].mean()]
같다("40", 절편, [9, -1], "공장별 절편")
같다("40", np.r_[절편[0] + 안 * xs[:3], 절편[1] + 안 * xs[3:]], ys, "잔차가 0")
print(f" 통합 {통합:+.3f} 인데 공장 안은 {안:+.3f}. 부호가 뒤집힌다.")[40] 맞다 차이 0.00e+00 여섯 점을 한꺼번에 놓은 기울기
rank M_D = 4 = 6 - 2
[40] 맞다 차이 2.22e-16 공장 안에서의 기울기
[40] 맞다 차이 1.78e-15 공장별 절편
[40] 맞다 차이 8.88e-16 잔차가 0
통합 -1.143 인데 공장 안은 +1.000. 부호가 뒤집힌다.
bv = np.array([1.0, 2, 3])
Sg = np.outer(bv, bv) + np.eye(3)
같다("41", Sg, [[2, 2, 3], [2, 5, 6], [3, 6, 10]], "요인모형 공분산")
같다("41", np.sort(np.linalg.eigvalsh(Sg)), [1, 1, 15], "고윳값 15, 1, 1")
Si2 = np.eye(3) - np.outer(bv, bv) / 15
같다("41", Si2, np.linalg.inv(Sg), "셔먼-모리슨")
one = np.ones(3)
w41 = Si2 @ one / (one @ Si2 @ one)
같다("41", w41, [1, 1/3, -1/3], "최소분산 비중")
같다("41", w41 @ Sg @ w41, 5/3, "최소분산")
같다("41", (w41 @ bv) ** 2 + w41 @ w41, 5/3, "요인 몫 + 잡음 몫")
print(f" 포트폴리오 요인 노출 {w41 @ bv:.4f} < 개별 반응 {bv}")[41] 맞다 차이 0.00e+00 요인모형 공분산
[41] 맞다 차이 3.55e-15 고윳값 15, 1, 1
[41] 맞다 차이 1.94e-16 셔먼-모리슨
[41] 맞다 차이 5.55e-17 최소분산 비중
[41] 맞다 차이 6.66e-16 최소분산
[41] 맞다 차이 2.22e-16 요인 몫 + 잡음 몫
포트폴리오 요인 노출 0.6667 < 개별 반응 [1. 2. 3.]
a4 = np.array([[2.0, 0], [0, 1], [-2, 0], [0, -1]])
b4 = np.array([[1.2, 1.6], [-0.8, 0.6], [-1.2, -1.6], [0.8, -0.6]])
Mp = b4.T @ a4
같다("42", Mp, [[4.8, -1.6], [6.4, 1.2]], "M = sum b_i a_i^T")
U, s, Vt = np.linalg.svd(Mp)
같다("42", s, [8, 2], "특이값")
R = U @ Vt
같다("42", R, [[0.6, -0.8], [0.8, 0.6]], "최적 회전")
print(f" 각도 {np.degrees(np.arccos(R[0, 0])):.2f}도, det R = {np.linalg.det(R):+.1f}")
같다("42", a4 @ R.T, b4, "R a_i = b_i 로 정확히 맞는다")[42] 맞다 차이 0.00e+00 M = sum b_i a_i^T
[42] 맞다 차이 0.00e+00 특이값
[42] 맞다 차이 1.11e-16 최적 회전
각도 53.13도, det R = +1.0
[42] 맞다 차이 2.22e-16 R a_i = b_i 로 정확히 맞는다
print(f"점을 행에 쌓았으므로 회전은 a -> R a, 행렬로는 A R^T 다.")
print(f" ||A R^T - B||_F^2 = {np.linalg.norm(a4 @ R.T - b4)**2:.4f} <- 옳다")
print(f" ||A R - B||_F^2 = {np.linalg.norm(a4 @ R - b4)**2:.4f} <- 전치를 빠뜨리면")
잘못 = U @ Vt
print(f"\n전치를 빠뜨리고 SVD 를 돌리면 최적해가 R^T 로 나온다 :")
M잘못 = a4.T @ b4
U2, s2, Vt2 = np.linalg.svd(M잘못)
R잘못 = U2 @ Vt2
print(f" 얻은 R = \n{np.round(R잘못, 4)}")
print(f" 각도 {np.degrees(np.arctan2(R잘못[1,0], R잘못[0,0])):.2f}도 "
f"(옳은 답은 {np.degrees(np.arctan2(R[1,0], R[0,0])):.2f}도)")
print(" 부호만 다르다. 반대 방향으로 같은 각도만큼 돌린 것이다.")점을 행에 쌓았으므로 회전은 a -> R a, 행렬로는 A R^T 다.
||A R^T - B||_F^2 = 0.0000 <- 옳다
||A R - B||_F^2 = 25.6000 <- 전치를 빠뜨리면
전치를 빠뜨리고 SVD 를 돌리면 최적해가 R^T 로 나온다 :
얻은 R =
[[ 0.6 0.8]
[-0.8 0.6]]
각도 -53.13도 (옳은 답은 53.13도)
부호만 다르다. 반대 방향으로 같은 각도만큼 돌린 것이다.
def 반응표(p0, u1, u2):
p0, u1, u2 = map(np.asarray, (p0, u1, u2))
assert abs(u1 @ p0) < 1e-12 and abs(u2 @ p0) < 1e-12
return np.outer(u1, u1) + np.outer(u1, u2) + np.outer(u2, u2)
for 이름, p0, u1, u2 in [
("모든 값이 양수 p=(1,2,3)", [1.0, 2, 3], [2.0, -1, 0], [3.0, 0, -1]),
("셋째가 자유재 p=(1,2,0)", [1.0, 2, 0], [2.0, -1, 0], [0.0, 0, 1])]:
J = 반응표(p0, u1, u2)
p0 = np.asarray(p0)
ok = np.allclose(J @ p0, 0) and np.allclose(J.T @ p0, 0)
print(f"\n{이름} (J p = J^T p = 0 인가: {ok})")
print(f" J =\n{J}")
for 버릴 in range(3):
남 = [r for r in range(3) if r != 버릴]
d_ = np.linalg.det(J[np.ix_(남, [0, 1])])
표 = " <- 못 푼다" if abs(d_) < 1e-9 else ""
print(f" {버릴+1}행을 버리고 p3 를 고정하면 det = {d_:+.4f}{표}")
모든 값이 양수 p=(1,2,3) (J p = J^T p = 0 인가: True)
J =
[[19. -2. -5.]
[-5. 1. 1.]
[-3. 0. 1.]]
1행을 버리고 p3 를 고정하면 det = +3.0000
2행을 버리고 p3 를 고정하면 det = -6.0000
3행을 버리고 p3 를 고정하면 det = +9.0000
셋째가 자유재 p=(1,2,0) (J p = J^T p = 0 인가: True)
J =
[[ 4. -2. 2.]
[-2. 1. -1.]
[ 0. 0. 1.]]
1행을 버리고 p3 를 고정하면 det = +0.0000 <- 못 푼다
2행을 버리고 p3 를 고정하면 det = +0.0000 <- 못 푼다
3행을 버리고 p3 를 고정하면 det = +0.0000 <- 못 푼다
값이 0인 재화를 기준으로 잡으면 어느 줄을 버려도 못 푼다. 나머지 두 줄이 서로 비례하기 때문이다.
기준은 반드시 양의 값을 갖는 재화로 골라야 한다.
옛 표에 0이 있으면 RAS 가 안 될 수 있다¶
for 이름, Z0v, 목행, 목열 in [
("보통 표", [[10.0, 20], [30, 40]], [40.0, 50], [50.0, 40]),
("대각만 있는 표", [[10.0, 0], [0, 40]], [40.0, 50], [50.0, 40])]:
Z0v = np.array(Z0v)
Zr = Z0v.copy()
for _ in range(200):
Zr = Zr * (np.array(목행) / Zr.sum(axis=1))[:, None]
Zr = Zr * (np.array(목열) / Zr.sum(axis=0))[None, :]
행오차 = np.abs(Zr.sum(axis=1) - 목행).max()
열오차 = np.abs(Zr.sum(axis=0) - 목열).max()
print(f"{이름:>16} : 200회 뒤 행 오차 {행오차:.4f}, 열 오차 {열오차:.4f}"
f"{' <- 수렴 안 함' if max(행오차, 열오차) > 1e-6 else ''}")
print("\n대각만 있는 표는 행 합과 열 합이 같을 수밖에 없다.")
print("목표가 (40,50) 과 (50,40) 이면 40 = 50 이라는 모순이라 답이 없다.") 보통 표 : 200회 뒤 행 오차 0.0000, 열 오차 0.0000
대각만 있는 표 : 200회 뒤 행 오차 10.0000, 열 오차 0.0000 <- 수렴 안 함
대각만 있는 표는 행 합과 열 합이 같을 수밖에 없다.
목표가 (40,50) 과 (50,40) 이면 40 = 50 이라는 모순이라 답이 없다.
4. 자기채점¶
정답 = {
23: 1, 24: 2, 25: 1, 26: 1, 27: 2, 28: 2, 29: 2, 30: 2, 31: 2, 32: 2,
33: 2, 34: 1, 35: 5, 36: 3, 37: 5, 38: 3, 39: 3, 40: 3, 41: 5, 42: 3,
43: 4, 44: 1,
}
힌트 = {
1: "줄 수와 미지수 수가 같고 독립인지 보라. 미지수를 발명해 맞춘 경우도 여기다",
2: "줄이 모자라 답이 평면으로 남는 자리를 보라",
3: "줄이 미지수보다 많아 다 맞힐 수 없는 자리를 보라",
4: "같은 조작이 되풀이되는지, 장기 거동을 묻는지 보라",
5: "부피, 곡률, 특이값 같은 크기를 묻는지 보라",
}
내답 = dict(정답) # 여기에 자기 답을 적어 넣는다
맞은 = 0
for 번호 in sorted(정답):
if 내답.get(번호) == 정답[번호]:
맞은 += 1
else:
print(f" 문제 {번호} : 상자 {내답.get(번호)} 라고 했는데 "
f"{정답[번호]} 다. {힌트[정답[번호]]}")
print(f"\n{맞은}/{len(정답)} 맞음")
from collections import Counter
print("상자별 분포 :", dict(sorted(Counter(정답.values()).items())))
print(" 상자 2 가 많은 것이 3막의 성격이다 — 답이 하나로 안 정해진다")
22/22 맞음
상자별 분포 : {1: 5, 2: 8, 3: 5, 4: 1, 5: 3}
상자 2 가 많은 것이 3막의 성격이다 — 답이 하나로 안 정해진다
정리¶
| 확인한 것 | 결과 |
|---|---|
| 게임의 값 를 미지수에 넣기 | 정방계, 양쪽이 같은 |
| 기저해 8개 중 꼭짓점 5개 | 전수 조사로 일치 |
| 각 는 단위행렬과 한 열만 다르다 | |
| 다시 푼 이익과 예측이 일치 | |
| 수송표·RAS·배정 세 번 다 | |
| 합계 사상의 영공간 | , = 2~5 전부 확인 |
| 최소분산이 정말 최소인가 | 20만 표본 중 더 작은 것 0개 |
| 같은 제약, 다른 기준 | 넷 다 296을 맞추는데 크기가 다 다르다 |
| 한계비용 균등 | |
| 는 부정부호, 는 음정부호 | 제약 위의 곡률은 다른 물음 |
| 대칭 사영이 가장 가깝다 | 무작위 대칭행렬 2만 개 중 더 가까운 것 0개 |
| FWL 두 단계 = 한 번에 풀기 | 계수 2.0으로 일치 |
| 심슨 역설 | 통합 -1.143, 공장 안 +1.000 |
| 전치를 빠뜨리면 | 회전이 반대 방향으로 나온다 |
| 자유재를 기준으로 잡으면 | 어느 줄을 버려도 못 푼다 |
| 옛 표에 0이 많으면 | RAS 가 수렴하지 않는다 |
마지막 세 줄이 이 노트북에서 제일 값지다.
답이 안 나오는 것보다 그럴듯한 답이 나오는 쪽이 위험하다. 전치를 빠뜨려도 회전행렬이 나오고 각도도 그럴듯하다. 부호만 반대일 뿐이다.