N6 · SVD·행렬함수·행렬방정식을 실제로 계산하기#
1. 충격 이후의 경로와 누적 크기를 함께 알고 싶다#
두 부문의 기준 상태에서 벗어난 무차원 지수 \(x(t)\)가
를 따른다고 합시다. 시간은 정한 한 기간 단위이고 \(A\)의 계수 단위는 그 역수입니다. 둘째 부문의 변화는 더 빠르게 줄며 첫 부문에 영향을 줍니다. 초기 충격 외에 새 투입은 없습니다. 둘째 식부터 풀면 \(x_2(t)=e^{-2t}x_2(0)\)이고, 첫 식에 \(e^t\)를 곱해 적분하면
경로를 계산하는 데에는 행렬지수가 필요합니다. 누적 충격의 제곱크기
도 알고 싶다면 \(J(x_0)=x_0^TPx_0\)인 행렬 \(P\)를 구하면 됩니다. \(J\)의 단위는 지수 제곱에 시간을 곱한 것입니다. 아래에서 \(P\)를 Lyapunov 방정식으로 구하고 위 적분과 일치함을 보겠습니다.
이런 계산은 값을 정의하는 공식과 실제 알고리즘을 구별해야 합니다. \(e^A=V e^\Lambda V^{-1}\)는 대각화 가능할 때의 식이지만, \(V\)가 거의 특이하면 중간 오차를 키울 수 있습니다. \(AX+XB=C\)를 벡터화하는 것도 맞지만, 작은 행렬방정식을 훨씬 큰 선형계로 바꿀 수 있습니다. 이번 단원은 같은 답을 얻으면서 구조를 보존하는 계산을 다룹니다.
2. SVD를 A의 교차곱으로 만들면 작은 값은 어떻게 되는가#
N4의 \(L_\varepsilon=\begin{pmatrix}1&1\\\varepsilon&0\\0&\varepsilon\end{pmatrix}\)를 다시 봅시다. 그 작은 특이값은 \(\varepsilon\)입니다. \(L_\varepsilon^TL_\varepsilon\)의 작은 고윳값은 \(\varepsilon^2\)인데, \(\varepsilon=10^{-8}\)에서 계산한 교차곱은 이미 계수 1이 되었습니다. 고유분해를 아무리 정확히 해도 잃은 \(\varepsilon^2\)를 복원할 수 없습니다.
성분별 곱셈오차는 N1에서 \(|\widehat G-A^TA|\le\gamma_m|A|^T|A|\)입니다. 따라서 \(\|\widehat G-A^TA\|_2\le\gamma_m\|A\|_F^2\)라는 충분한 상한을 얻습니다. 이 잡음 규모가 \(\sigma_n^2\)와 비슷하면 작은 고윳값의 상대정확도를 보장하기 어렵습니다. 대략 \(\sigma_n/\sigma_1\sim\sqrt u\)가 경계로 등장하지만, 차원과 성분 구조에 따라 달라집니다. 모든 작은 특이값이 반드시 소멸하는 필요충분조건은 아닙니다. \(\operatorname{diag}(1,10^{-12})\)처럼 서로 섞이지 않는 대각 성분은 제곱값도 표현할 수 있습니다.
직접 SVD의 노름 후진오차가 \(O(u)\|A\|\)이면 특이값의 절대오차도 그 규모로 제한됩니다. 그러나 \(\sigma_i\ll\|A\|\)에서는 이것만으로 작은 상대오차를 보장하지 않습니다. 원래 행렬의 형성오차와, 이미 주어진 특수 구조행렬에서 작은 값을 계산하는 상대정확도는 별도로 봐야 합니다.
그림 103 정확한 작은 특잇값은 \(\varepsilon\)이다. 교차곱의 음의 계산 고윳값은 제곱근의 유효 입력이 아니므로 실패로 표시하며, 0은 상대오차 1이다. 직접 SVD의 좋은 결과도 모든 행렬의 성분별 상대정확도 보장은 아니다.#
3. 두 쪽에서 Householder를 적용해 이중대각으로#
\(m\ge n\)인 \(A\)에 왼쪽 Householder를 적용해 첫 열의 첫 행 아래를 없앱니다. 이어 오른쪽 Householder로 첫 행의 둘째 열 뒤를 없앱니다. 다음 행·열에서 같은 일을 반복하면
인 상 이중대각형을 얻습니다. 이전 행·열의 0을 건드리지 않는 영역에서만 다음 변환을 적용합니다. 왼쪽·오른쪽 인자는 모두 직교하므로 특이값을 보존합니다. 고유문제의 닮음변환과 달리 양쪽 인자가 서로 달라도 됩니다.
\(B=\widetilde U\Sigma\widetilde V^T\)를 구하면 원래 인자는 \(U\widetilde U\), \(V\widetilde V\)입니다. \(m\gg n\)이면 먼저 축소 QR의 작은 \(n\times n\) 행렬 \(R\)을 만든 뒤 이중대각화할 수도 있습니다. \(A^TA=R^TR\)라는 대수적 관계를 이용하지만 실제로 이 교차곱을 형성할 필요는 없습니다.
import numpy as np
import scipy.linalg as la
def reflector(x):
x = np.array(x, dtype=float, copy=True)
scale = np.max(np.abs(x))
if scale == 0: return None
v = x/scale
v[0] += np.copysign(np.linalg.norm(v), v[0])
return v/np.linalg.norm(v)
def bidiagonalize(A):
B = np.array(A, dtype=float, copy=True)
m, n = B.shape
if m < n: raise ValueError("이 구현은 m >= n을 사용합니다")
U = np.eye(m); V = np.eye(n)
for k in range(n):
v = reflector(B[k:, k])
if v is not None:
B[k:, k:] -= 2*np.outer(v, v@B[k:, k:])
U[:, k:] -= 2*np.outer(U[:, k:]@v, v)
B[k+1:, k] = 0.
if k+1 < n:
w = reflector(B[k, k+1:])
if w is not None:
B[k:, k+1:] -= 2*np.outer(B[k:, k+1:]@w, w)
V[:, k+1:] -= 2*np.outer(V[:, k+1:]@w, w)
B[k, k+2:] = 0.
return U, B, V
A = np.array([[1., 2., 0.], [0., 1., 3.], [2., -1., 1.], [1., 0., 2.]])
U, B, V = bidiagonalize(A)
assert np.linalg.norm(U@B@V.T-A) < 1e-13
assert np.allclose(la.svdvals(B), la.svdvals(A), rtol=1e-13)
print("양쪽 직교변환과 이중대각 특잇값 보존 확인")
양쪽 직교변환과 이중대각 특잇값 보존 확인
이중대각의 SVD에는 암묵 QR, 분할정복, 높은 상대정확도를 목표로 한 dqds 등의 전문 알고리즘을 사용합니다. Demmel–Kahan의 상대정확도 결과는 이중대각의 성분에 대한 상대 섭동과 해당 알고리즘 조건을 함께 갖는 정리입니다. 일반 밀집행렬의 작은 특이값을 앞선 축약오차와 무관하게 모두 상대 \(u\)로 계산한다는 뜻은 아닙니다. 여기서는 이중대각화와 보존 성질을 직접 구현·증명하고 최종 SVD는 LAPACK의 SVD 알고리즘에 연결합니다.
4. 지수급수는 맞지만 몇 항을 어떻게 더할 것인가#
1절의 \(A\)는 작은 차수의 잘 분리된 고윳값이라 지수의 손계산이 간단합니다. 비대각 2를 매우 큰 \(M\)으로 바꾸어도 이 특정한 나눗셈차 \(M(e^{-1}-e^{-2})\)는 큰 두 수의 거의 완전한 상쇄가 아닙니다. 고유벡터 조건수가 크다는 사실만으로 그 실험의 상대오차가 반드시 \(M\)에 비례해 커진다고 단정하지 않습니다.
반면 H8의 \(A_\varepsilon=\begin{pmatrix}-1&1\\0&-1-\varepsilon\end{pmatrix}\)에서는 \((e^{-1}-e^{-1-\varepsilon})/\varepsilon\)가 나타나 가까운 두 지수의 차가 민감합니다. 이 예는 expm1으로 안정화할 수 있었지만 모든 행렬의 지수를 그런 한 성분 공식으로 풀 수는 없습니다.
또 다른 문제는 지수급수의 큰 중간항입니다. \(K=\begin{pmatrix}-49&24\\-64&31\end{pmatrix}\)의 고윳값은 \(-1,-17\)입니다. \(e^K\)는 유한한 작은 크기지만 \(K^j/j!\)에는 큰 양·음 성분이 나타납니다. 많은 큰 항을 더해 작은 답을 만들면 누적오차가 남을 수 있고, 20항이라는 고정 개수는 충분하다는 보장도 없습니다.
그림 104 기준값은 고정밀 연산으로 구했다. 가로축은 Taylor 절단 차수이고 세로축은 상대 Frobenius 오차의 로그 눈금이다. 두 수평선은 고정된 다른 알고리즘의 결과로 Taylor 항 수에 의존하지 않는다.#
5. 작게 만든 뒤 유리함수로 계산하고 제곱한다#
스케일링과 제곱은 \(X=A/2^s\)를 작게 만든 후 \(e^A=(e^X)^{2^s}\)를 사용하는 방식입니다. 작은 \(X\)에서 다항식보다 적은 차수로 지수급수를 맞추는 Padé 유리함수를 사용합니다. 예를 들어
는 지수급수의 6차까지 맞춥니다. 역행렬 표기는 분모 행렬을 계수로 하는 선형계 풀이로 구현합니다. 정확한 산술에서는 \(e^X-r_3(X)=O(\|X\|^7)\)입니다. 아래 고정 차수 코드는 원리를 확인하는 교육용이며, 라이브러리의 차수·스케일 적응 선택을 복제한 것은 아닙니다.
def exp_pade3(A, threshold=.05):
A = np.asarray(A, dtype=float)
norm = np.linalg.norm(A, 1)
squarings = max(0, int(np.ceil(np.log2(norm/threshold)))) if norm else 0
X = A/(2.**squarings); I = np.eye(len(A))
X2 = X@X; X3 = X2@X
numerator = I+X/2+X2/10+X3/120
denominator = I-X/2+X2/10-X3/120
result = la.solve(denominator, numerator)
for _ in range(squarings): result = result@result
return result
A = np.array([[-1., 2.], [0., -2.]])
exact = np.array([[np.exp(-1), 2*(np.exp(-1)-np.exp(-2))], [0., np.exp(-2)]])
assert np.linalg.norm(exp_pade3(A)-exact)/np.linalg.norm(exact) < 2e-11
assert np.linalg.norm(la.expm(A)-exact)/np.linalg.norm(exact) < 1e-14
print("정의식·교육용 Padé·표준 expm 대조")
정의식·교육용 Padé·표준 expm 대조
표준 expm은 행렬 크기와 오차 추정에 따라 Padé 차수와 스케일을 선택합니다. SciPy의 expm 문서는 사용하는 알고리즘의 출처를 밝힙니다. 정확한 유리근사의 후진오차와, 실제 선형계 풀이·제곱 단계에서 생기는 오차는 구별해야 합니다. 지나친 스케일링은 더 많은 제곱을 요구하며 오차와 비용을 늘릴 수 있습니다. 전체 계산에 무조건 \(\|E\|\le u\|A\|\)인 \(e^{A+E}\) 표현이 있다고 단정하지 않습니다.
6. Schur 형식에서 일반 행렬함수를 계산하기#
\(A=QTQ^*\)이면 \(f(A)=Qf(T)Q^*\)입니다. \(F=f(T)\)는 \(T\)와 가환하므로 \(TF=FT\)입니다. 대각은 \(f(t_{ii})\)이고 \(i<j\)에서
대각에서 가까운 띠부터 계산하면 오른쪽의 이전 값들을 알고 있습니다. 이것이 스칼라 Parlett 점화입니다. \(t_{ii}\)와 \(t_{jj}\)가 같거나 매우 가까우면 분모가 문제가 됩니다. 이런 고윳값들을 같은 블록으로 묶고 블록 내부는 도함수·급수로 처리하며, 떨어진 블록 사이를 Sylvester 방정식으로 풉니다.
예컨대 \(T=\begin{pmatrix}1&1\\0&1\end{pmatrix}\)에서 \(\log T=\begin{pmatrix}0&1\\0&0\end{pmatrix}\)입니다. 스칼라 점화의 \(0/0\)은 함수가 없어서가 아니라 같은 고윳값의 도함수를 사용해야 한다는 뜻입니다. H8의 \(\log(1+z)=z+O(z^2)\)와 \(N^2=0\)이 정확한 답을 줍니다.
상삼각 제곱근 \(S^2=T\)에서는 대각 \(s_{ii}=\sqrt{t_{ii}}\)의 분지를 정한 뒤 \(s_{ij}=(t_{ij}-\sum_{i<k<j}s_{ik}s_{kj})/(s_{ii}+s_{jj})\)를 씁니다. 주제곱근을 쓰려면 스펙트럼이 음의 실수축과 0을 피하는 등의 분지 조건을 확인합니다. 로그를 구했다고 그 결과가 마르코프 생성자인 것도 아닙니다. 실수성, 행합 0, 비대각 비음수 같은 모형 제약은 별도로 검사해야 합니다.
7. 행렬방정식의 각 성분을 풀기#
Sylvester 방정식 \(AX+XB=C\)에서 \(A\)는 \(m\times m\), \(B\)는 \(n\times n\), 미지수 \(X\)는 \(m\times n\)입니다. 벡터화하면 \((I_n\otimes A+B^T\otimes I_m)\operatorname{vec}X=\operatorname{vec}C\)입니다. 맞는 식이지만 \(m=n=N\)이면 계수행렬이 \(N^2\times N^2\)가 되어 밀집 저장은 \(O(N^4)\), 일반 소거는 \(O(N^6)\)입니다.
대신 \(A=UTU^*\), \(B=VSV^*\)를 Schur 분해하고 \(Y=U^*XV\), \(F=U^*CV\)로 놓으면 \(TY+YS=F\)입니다. 성분식은
행은 아래에서 위로, 열은 왼쪽에서 오른쪽으로 진행하면 필요한 값이 이미 계산되어 있습니다. 모든 \(t_{ii}+s_{jj}\)가 비영이어야 유일한 해를 얻습니다. 각 성분의 두 합 길이는 \(m+n\) 이하이므로 총 \(O(mn(m+n))\)이며 정사각에서는 Schur 분해까지 \(O(N^3)\)입니다.
def sylvester_schur(A, B, C):
T, U = la.schur(A, output="complex")
S, V = la.schur(B, output="complex")
F = U.conj().T@C@V
m, n = F.shape; Y = np.zeros_like(F)
for j in range(n):
for i in range(m-1, -1, -1):
denominator = T[i,i]+S[j,j]
if denominator == 0: raise ValueError("유일한 Sylvester 해가 없습니다")
rhs = F[i,j]-T[i,i+1:]@Y[i+1:,j]-Y[i,:j]@S[:j,j]
Y[i,j] = rhs/denominator
return U@Y@V.conj().T
A = np.array([[-1., 2.], [0., -2.]])
P = sylvester_schur(A.T, A, -np.eye(2))
expected = np.array([[.5, 1/3], [1/3, 7/12]])
assert np.linalg.norm(P-expected) < 2e-14
assert np.linalg.norm(A.T@P+P@A+np.eye(2)) < 2e-14
print("누적 충격 행렬 P:", np.real_if_close(P))
누적 충격 행렬 P: [[0.5 0.33333333]
[0.33333333 0.58333333]]
분모가 모두 비영이라는 것은 유일성 조건입니다. 좋은 조건수를 뜻하지는 않습니다. 실제 민감도는
로 측정할 수 있고, 잔차 \(R=C-A\widehat X-\widehat XB\)에 대해 \(\|X-\widehat X\|_F\le\|R\|_F/\operatorname{sep}(A,-B)\)입니다. 비정규행렬에서는 고윳값 합의 최소 절댓값만으로 이 분리를 대신할 수 없습니다.
8. 처음의 누적 충격과 Lyapunov 방정식#
\(P=\begin{pmatrix}p&q\\q&r\end{pmatrix}\)를 1절의 \(A\)에 넣으면
첫 식에서 \(p=1/2\), 둘째에서 \(q=1/3\), 셋째에서 \(r=7/12\)입니다. \(p>0\), \(\det P=13/72>0\)이므로 양의 정부호입니다. 경로를 따라 \(d(x^TPx)/dt=x^T(A^TP+PA)x=-\|x\|^2\)이고, \(x(t)\to0\)이므로 적분하면 \(x_0^TPx_0=J(x_0)\)입니다. 특히 둘째 부문에 단위 충격을 주면 누적 제곱크기는 \(7/12\)입니다.
일반적으로 \(A\)가 연속시간 안정이고 \(Q\succ0\)이면 \(A^TP+PA=-Q\)의 해는 \(P=\int_0^\infty e^{tA^T}Qe^{tA}\,dt\succ0\)입니다. 이산시간 안정 \(\rho(A)<1\)이면 \(P-A^TPA=Q\)의 해는 \(\sum_{k\ge0}(A^T)^kQA^k\)입니다. \(Q\succeq0\)만 있으면 양의 준정부호는 얻지만 양의 정부호를 무조건 얻지는 않습니다. \(Q=0\)은 간단한 반례입니다.
그림 105 시간은 기간 단위이며 왼쪽은 무차원 지수 변화, 오른쪽은 시간까지 적분한 제곱크기이다. 유한 적분은 무한 적분값 아래에서 접근한다.#
9. 일반화 방정식과 반복정련#
일반화 Sylvester \(AXB+CXD=E\)도 두 행렬쌍의 QZ로 줄일 수 있습니다. \(Q_L^*AZ_L=S_A\), \(Q_L^*CZ_L=S_C\), \(Q_R^*BZ_R=T_B\), \(Q_R^*DZ_R=T_D\)를 삼각으로 만들고 \(X=Z_LYQ_R^*\)로 놓습니다. 그러면
\(j\)번째 열에서 \((t_{B,jj}S_A+t_{D,jj}S_C)y_j\)를 미지항으로 모으고, 이전 열의 두 합을 우변으로 넘깁니다. 왼쪽은 상삼각이므로 후진대입합니다. 모든 \(s_{A,ii}t_{B,jj}+s_{C,ii}t_{D,jj}\)가 비영이면 유일합니다. 두 쌍이 각각 정칙이라는 조건만으로 이 교차 합까지 비영이라고 보장하지는 않습니다. 실제 QZ를 이용한 이 변환은 마지막 정리에서 확인합니다.
또한 낮은 정밀도에서 만든 근사 인자 \(A+E\)를 재사용해 높은 정밀도로 잔차를 계산하는 반복정련이 있습니다. 이상화한 보정은 \((A+E)d_k=b-Ax_k\), \(x_{k+1}=x_k+d_k\)입니다. 오차 \(e_k=x_*-x_k\)는
를 따릅니다. 수렴은 이 반복행렬의 스펙트럼 반지름이 1 미만인 조건과 연결됩니다. 단지 “\(\kappa(A)u_{\rm low}<1\)이면 언제나 높은 정밀도의 전진오차”라고 말하면 인자의 성장·상수·잔차 반올림·문제 조건수를 빠뜨립니다. 잔차오차와 보정오차가 매번 \(\xi_k\)로 들어가면 \(\|e_{k+1}\|\le q\|e_k\|+\|\xi_k\|\)이고, \(\|\xi_k\|\le\xi\)이면 오차의 바닥은 \(\xi/(1-q)\)로 제한됩니다.
10. 연습과 전체 풀이#
1. \(L_\varepsilon^TL_\varepsilon\)가 binary64에서 작은 고윳값을 잃는 예가 직접 SVD의 상대오차 정리는 아닌 이유를 설명하세요.
풀이. 교차곱은 \(\varepsilon^2\)를 1에 더해 저장하므로 작은 정보가 사라집니다. 직접 SVD의 노름 후진오차는 \(A\) 크기에 비례하며 특이값 절대오차를 제한합니다. 이를 \(\varepsilon\)로 나누면 상대오차 상한은 커질 수 있습니다. 두 경로의 보장 범위가 다릅니다.
2. \(B=\begin{pmatrix}2&1\\0&1\end{pmatrix}\)의 Jordan–Wielandt 행렬을 쓰세요.
풀이. \(J=\begin{pmatrix}0&B^T\\B&0\end{pmatrix}=\begin{pmatrix}0&0&2&0\\0&0&1&1\\2&1&0&0\\0&1&0&0\end{pmatrix}\)입니다. \(Bv_i=\sigma_i u_i\), \(B^Tu_i=\sigma_i v_i\)에서 \((v_i,u_i)\)와 \((v_i,-u_i)\)는 각각 \(\sigma_i,-\sigma_i\)의 고유벡터입니다. 절댓값을 취하면 각 특이값이 두 번 나타난다는 중복도를 함께 기록해야 합니다.
3. \(A=\operatorname{diag}(1,2)\), \(B=\operatorname{diag}(-1,3)\)인 Sylvester 문제는 모든 \(C\)에서 유일한가요?
풀이. \((1,1)\)성분의 계수는 \(1-1=0\)입니다. \(c_{11}\ne0\)이면 해가 없고, \(c_{11}=0\)이면 \(x_{11}\)을 임의로 택할 수 있습니다. 나머지 성분 계수는 4,1,5이므로 각각 결정됩니다. 한 우변에서 해가 존재하는 것과 모든 우변의 유일성은 다릅니다.
4. 처음의 \(P\)에서 \(x_0=(1,1)\)의 누적 제곱크기를 구하세요.
풀이. \(x_0^TPx_0=1/2+2/3+7/12=7/4\)입니다. 각 부문 단독 충격값의 합 \(13/12\)에 양의 교차항 \(2/3\)이 더해집니다. 두 충격의 경로가 서로 영향을 주므로 단독 값만 합할 수 없습니다.
5. \(T=\begin{pmatrix}4&3\\0&1\end{pmatrix}\)의 주제곱근을 구하세요.
풀이. 대각은 2와 1이고 비대각은 \(3/(2+1)=1\)입니다. \(S=\begin{pmatrix}2&1\\0&1\end{pmatrix}\)를 제곱하면 \(T\)입니다. 반대 분지를 고르면 분모도 달라지므로 부호를 먼저 정해야 합니다.
6. 스칼라 \(A=1\), 근사 인자 \(A+E=1/4\)의 이상화 반복정련을 분석하세요.
풀이. \(E=-3/4\)이고 오차 반복계수는 \((A+E)^{-1}E=-3\)입니다. \(\|A^{-1}E\|=3/4<1\)이어도 반복은 일반적으로 발산합니다. 그 부등식은 역행렬의 존재에 충분하지만 이 정련의 수렴을 보장하는 충분한 작은 상한은 아닙니다.
7. \(N^2=0\)이면 \(r_3(N)\)와 \(e^N\)을 비교하세요.
풀이. 분자와 분모는 \(I\pm N/2\)로 줄고 \((I-N/2)^{-1}=I+N/2\)입니다. 곱은 \(I+N\)이며 지수도 \(I+N\)입니다. 차수가 고정된 근사가 어떤 특수행렬에서는 정확할 수 있다는 예입니다.
11. 지금까지의 내용을 수학의 언어로 정리해 봅시다#
H8의 행렬함수 정의, N3의 직교변환 오차, H7의 일반화 Schur 형를 출발점으로 사용합니다. 지수와 Sylvester의 계산에서 필요한 나머지 연결을 증명합니다.
정리 1. 양쪽 직교축약과 특이값#
직교 \(U,V\)에 대해 \(B=U^TAV\)는 \(A\)와 같은 특이값을 갖습니다. 양쪽 Householder 이중대각화는 정상 범위에서 각 국소오차가 작으면 작은 노름 후진오차를 갖습니다.
증명. \(B^TB=V^TA^TUU^TAV=V^TA^TAV\)이므로 고윳값이 같고 그 비음수 제곱근인 특이값도 같습니다. 두 변환 단계마다 \(\widehat B_j=H_j\widehat B_{j-1}G_j+E_j\)로 쓰고 \(\|E_j\|_F\le\eta\|\widehat B_{j-1}\|_F\)라 합시다. 양쪽 직교곱이 Frobenius 노름을 보존하므로 N3의 망원합 증명을 양쪽에 그대로 적용하여 \(\|\Delta A\|_F\le((1+\eta)^k-1)\|A\|_F\)를 얻습니다. 여기서 \(k\)는 적용한 변환 쌍의 수입니다. 계산된 \(U,V\)와 이 후진식의 정확한 직교인자는 구별합니다. ∎
정리 2. 작은 Padé 근사의 후진 해석#
\(\|X\|=a<2\)인 유도노름에서 \(r_1(X)=(I-X/2)^{-1}(I+X/2)\)는 \(e^{X+E_X}\)로 표현되며
영행렬이 아닌 입력을 정확한 산술로 스케일링하여 \(X=A/2^s\)로 놓고 제곱하면 상대 후진오차가 \(a^2/[12(1-a^2/4)]\) 이하입니다.
증명. \(\|X/2\|<1\)이므로 두 로그급수가 절대수렴합니다. 같은 행렬의 거듭제곱은 가환하므로
이 급수의 지수를 취하면 \(r_1(X)\)입니다. 꼬리 노름에서 \((2j+1)\ge3\)을 쓰면 \(\sum_{j\ge1}a^{2j+1}/(3\cdot2^{2j})=a^3/[12(1-a^2/4)]\)입니다. \(2^s\)제곱은 \(e^{2^s(X+E_X)}\)이고 \(2^sa=\|A\|\)여서 상대 상한이 따릅니다. 이 정리는 정확한 유리함수와 제곱에 대한 것입니다. 반올림된 행렬곱의 오차는 여기에 추가해야 합니다. ∎
\(r_3\)의 6차 일치는 분모 다항식에 \(\sum_{j=0}^6X^j/j!\)를 곱해 확인합니다. 상수부터 3차는 분자가 되고 4·5·6차의 계수는 각각 0입니다. 예를 들어 4차는 \(1/24-1/12+1/20-1/120=0\)입니다. 분모가 0 근처에서 가역이므로 해석적 나머지는 \(O(\|X\|^7)\)입니다. 높은 차수의 최적 스케일 선택과 반올림 분석은 라이브러리 알고리즘의 추가 내용입니다.
정리 3. Sylvester의 유일성과 분리#
선형사상 \(\mathcal L(X)=AX+XB\)가 가역일 필요충분조건은 모든 고윳값 합 \(\lambda_i(A)+\lambda_j(B)\)가 비영인 것입니다. 가역이면 \(\|\mathcal L^{-1}\|_{F\to F}=1/\operatorname{sep}(A,-B)\)입니다.
증명. Schur 좌표로 바꾸는 \(X\mapsto U^*XV\)는 가역입니다. 상삼각 좌표에서 열을 왼쪽부터, 각 열은 아래부터 나열하면 7절 성분식의 계수행렬은 삼각이고 대각은 \(t_{ii}+s_{jj}\)입니다. 따라서 행렬식은 그 곱이며 가역성 조건을 얻습니다. 이는 일부 합이 0일 때 특정 우변에서 해가 없거나 여러 개일 수 있음을 포함합니다.
유한차원 단위구면에서 연속함수 \(\|\mathcal L(X)\|_F\)는 최소를 달성합니다. 가역이면 그 최소는 양수이고 \(\|\mathcal L(X)\|_F\ge\operatorname{sep}\|X\|_F\)입니다. \(Y=\mathcal L(X)\)를 대입하면 역노름 상한이며, 최소를 달성하는 \(X\)를 사용하면 등호입니다. 잔차에 이 역을 적용하면 7절의 오차 인증식이 나옵니다. ∎
정리 4. Lyapunov와 안정성#
실수 \(A\)의 모든 고윳값의 실수부가 음수이고 \(Q=Q^T\succ0\)이면 \(A^TP+PA=-Q\)의 유일한 해는 8절의 적분이며 \(P\succ0\)입니다. 역으로 그러한 \(P,Q\succ0\)가 있으면 \(A\)는 연속시간 안정입니다. 이산시간에서는 \(\rho(A)<1\)과 \(P-A^TPA=Q\)의 양의 정부호 해가 같은 관계를 갖습니다.
증명. Jordan 지수의 각 항은 \(t^j e^{\lambda t}\)입니다. \(\Re\lambda<0\)이면 적분에 필요한 제곱항이 수렴하고 \(e^{tA}\to0\)입니다. \(P\)를 적분으로 정해 \(A^TP+PA\)를 계산하면 적분함수의 도함수 \(d(e^{tA^T}Qe^{tA})/dt\)를 적분한 것이므로 경계값은 \(-Q\)입니다. \(x\ne0\)이면 \(x^Te^{tA^T}Qe^{tA}x>0\)이고 적분도 양수여서 \(P\succ0\)입니다. 유일성은 고윳값 합의 실수부가 음수라는 Sylvester 조건에서 따릅니다.
역으로 복소 고유벡터 \(Av=\lambda v\)에 식을 좌우로 곱하면 \(2\Re\lambda\,v^*Pv=-v^*Qv<0\)이고 \(v^*Pv>0\)이므로 \(\Re\lambda<0\)입니다.
이산 급수는 H8의 거듭제곱 감쇠에 의해 절대수렴합니다. \(A^TPA\)는 급수의 첫 항을 뺀 것이므로 \(P-A^TPA=Q\)입니다. 양의 정부호는 첫 항 \(Q\succ0\)에서 얻습니다. 동차해 \(D=A^TDA\)는 \(D=(A^T)^kDA^k\to0\)이므로 유일합니다. 역으로 고유벡터를 곱하면 \((1-|\lambda|^2)v^*Pv=v^*Qv>0\)여서 \(|\lambda|<1\)입니다. ∎
정리 5. 일반화 Sylvester의 삼각 열 풀이#
9절의 두 QZ 분해가 주어졌다고 하자. 선형사상 \(X\mapsto AXB+CXD\)가 가역일 필요충분조건은 모든 \(s_{A,ii}t_{B,jj}+s_{C,ii}t_{D,jj}\)가 비영인 것입니다. 이 경우 삼각 열 풀이의 비용은 \(O(m^2n+mn^2)\)입니다.
증명. \(A=Q_LS_AZ_L^*\), \(B=Q_RT_BZ_R^*\)와 나머지 두 식에 \(X=Z_LYQ_R^*\)를 대입합니다. 안쪽 유니터리 곱들이 항등으로 사라져 \(AXB+CXD=Q_L(S_AYT_B+S_CYT_D)Z_R^*\)입니다. 양쪽에 \(Q_L^*,Z_R\)를 곱하면 9절의 식입니다.
상삼각 오른쪽 인자 때문에 \(j\)열에서는 \(y_1,\ldots,y_j\)만 나타납니다. 이미 계산한 앞 열을 넘기면
왼쪽 삼각계의 대각이 바로 제시한 합입니다. 모든 합이 비영이면 차례로 유일하게 풀고, 하나라도 0이면 전체 삼각 연산자의 행렬식이 0입니다. 각 열의 두 선형결합은 \(O(mj)\), 두 행렬·벡터곱과 삼각풀이는 \(O(m^2)\)이므로 합하면 \(O(mn^2+m^2n)\)입니다. ∎
정리 6. 이상화 반복정련의 충분조건#
\(A\)와 \(A+E\)가 가역일 때 모든 초기값에서 정확해로 가는 이상화 정련의 필요충분조건은 \(\rho((A+E)^{-1}E)<1\)입니다. \(\eta=\|A^{-1}E\|<1/2\)는 충분조건입니다.
증명. 9절의 오차 재귀에 H8의 거듭제곱 수렴 정리를 적용하면 첫 동치입니다. \(A+E=A(I+A^{-1}E)\)이므로 \(\|(A+E)^{-1}E\|\le\eta/(1-\eta)<1\)입니다. 매번 오차항의 노름이 \(\xi\) 이하이면 재귀 부등식을 반복하여 \(\|e_k\|\le q^k\|e_0\|+\xi(1-q^k)/(1-q)\)를 얻습니다. 이것이 유한정밀도 정련에서 정지 기준과 오차 바닥을 함께 보아야 하는 이유입니다. ∎
행렬함수와 행렬방정식을 계산할 때도 분해 방법과 문제의 조건을 함께 보아야 합니다. 다음 장에서는 순환·Toeplitz·띠 구조가 저장량과 계산량을 어떻게 줄이는지 살펴봅니다. N7로 이어 읽기.