수치 알고리즘으로 근과 선형계를 안정적으로 푸는 방법

문제의 구조가 해법을 고른다

비선형 방정식의 근을 찾는 일과 선형계 Ax=b를 푸는 일은 서로 다른 수치 알고리즘을 요구한다. Newton-Raphson은 f(x)=0의 근을 반복적으로 근사하고, Gaussian Elimination은 선형계를 상삼각 형태로 바꾼 뒤 해를 구한다. LU와 Cholesky 분해는 선형계 풀이를 반복하거나 행렬의 구조를 활용해야 할 때 특히 유용하다.

선택의 기준은 행렬이 대칭 양의 정부호(SPD)인지, 희소한지, 같은 계수 행렬에 여러 우변을 적용하는지, 그리고 Newton 반복에서 초기값과 도함수를 얼마나 신뢰할 수 있는지에 달려 있다.

Newton-Raphson은 근방에서 빠르지만 경로를 관리해야 한다

Newton-Raphson은 다음 반복식으로 비선형 방정식의 해를 갱신한다.

x_{k+1} = x_k − f(x_k)/f′(x_k)

적절한 초기값에서 근에 가까워지면 2차 수렴 특성을 보인다. 반면 도함수가 0에 가까운 지점, 좋지 않은 초기값, 감쇠 없이 지나치게 큰 스텝을 쓰는 경우에는 발산할 수 있다. 다변수 문제에서는 야코비안 기반으로 선형계 JΔx = −f를 풀어야 하므로, 반복 1회의 비용은 주로 O(n^3) 선형계 풀이가 지배한다.

Gaussian Elimination과 LU는 같은 계산 경로를 공유한다

Gaussian Elimination은 전진 소거로 행렬을 상삼각화하고, 후진 대입으로 해를 얻는다. 실제 구현에서는 부분 피벗팅을 포함한 LU 분해와 같은 절차로 이해할 수 있다.

부분 피벗팅은 성장인자를 억제해 수치 안정성을 높인다. 완전 피벗팅도 가능하지만, 추가 비용과 복잡도를 감수해야 한다. 해의 존재와 유일성은 행렬의 조건수와 순위에 좌우되며, 분해가 끝난 뒤의 해 계산은 삼각계 풀이로 결정적으로 진행된다.

SPD 행렬에는 Cholesky가 맞는다

LU 분해는 일반 정방행렬에 대해 A=LU 형태로, 피벗 행렬을 포함하면 A=PLU 형태로 적용한다. 같은 계수 행렬에 여러 우변 b를 적용할 때는 분해를 한 번만 수행하고 이후 풀이에 재사용할 수 있다.

Cholesky 분해는 SPD 행렬에 대해 A=LLᵀ로 표현한다. Dense 행렬 기준으로 LU는 약 2/3 n^3 플롭, Cholesky는 약 1/3 n^3 플롭이 필요하다. 대칭 저장도 가능해 메모리 사용량을 줄일 수 있으며, SPD 조건 아래에서는 수치적으로 매우 안정적이다. 다만 SPD 조건이 깨지면 분해는 실패한다.

수렴과 오차를 함께 감시하는 반복 과정

Newton-Raphson에서는 |f(x)|와 |Δx|를 함께 정지 조건으로 쓰고, 최대 반복 횟수를 둔다. 도함수 소멸과 NaN도 감지해야 하며, 필요한 경우 이분법이나 선 탐색으로 전환할 수 있어야 한다. 라인서치, 댐핑, 트러스트 리전은 큰 스텝으로 인한 발산을 줄이는 데 사용한다. 스케일링과 정규화도 수렴성 개선에 기여한다.

⤢✕

선형계에서는 잔차 ||Ax−b||/||b||와 조건수 추정을 확인해야 한다. 분해 실패나 경고를 처리하는 루틴도 운영 코드에 포함하는 편이 안전하다.

비용과 적용 조건의 차이

알고리즘	성능(연산비용)	확장성(대규모/희소)	일관성(결정성/수렴)	안정성(수치 오차)	운영 편의(재사용/라이브러리)
Newton-Raphson	1회 반복 저비용(1D), 다변수는 O(n^3) 선형계 포함	희소 야코비안 활용 시 유리	초기값·함수 특성에 좌우, 수렴 불확실성	감쇠/트러스트 리전 필요	구현 중간 난이도, 도함수 필요
Gaussian Elimination	Dense O(n^3)	희소 시 fill-in 관리 필수	결정적 결과	부분 피벗으로 양호	표준 구현 풍부, 교육·디버깅 용이
LU 분해	Dense ~2/3 n^3, 다중 b에 유리	희소 LU로 대규모 가능	결정적 결과	피벗팅으로 안정	분해 재사용 매우 편리
Cholesky 분해	Dense ~1/3 n^3, 메모리 절감	SPD 희소 행렬에 최적	결정적 결과	매우 안정(SPD 필요)	운영 간편, 캐시 친화

삼각계 풀이는 O(n^2)이므로, 다중 우변을 처리할 때는 LU 또는 Cholesky 분해를 한 번 수행한 후 반복 solve하는 방식이 총 비용을 낮춘다. 희소 행렬에서는 밴드 구조와 fill-in을 고려해야 하며, 노드 순서화와 AMD, METIS 등을 통해 대규모 시스템의 확장성을 확보할 수 있다.

해석, 최적화, 제어에서의 사용 방식

구조해석, 열전달, 회로 시뮬레이션에는 대규모 SPD 선형계가 등장하며, Cholesky 또는 LDLᵀ 기반 희소 솔버를 적용할 수 있다. 비선형 재료와 접촉 문제에서는 Newton-Raphson과 스파스 LU·Cholesky 선형계 풀이를 결합해 운용한다.

로지스틱 회귀, Newton-CG, Gauss-Newton 기반 파라미터 추정에서는 해스시안 근사와 선형계 풀이가 병목이 된다. 사전조건자를 활용할 수 있으며, 정규화 문제에서 AᵀA가 SPD이면 Cholesky로 빠르게 풀 수 있다. 다만 조건수가 극단적인 경우에는 QR 또는 SVD를 고려한다.

금융 모델의 파라미터 추정과 금리 곡선 부트스트랩에는 Newton-Raphson을 적용할 수 있고, 라인서치로 안정성을 확보한다. 배치 가격 산출처럼 같은 계수 행렬을 반복해서 다룰 때는 LU 분해 재사용이 처리량 향상으로 이어진다.

실시간 제어와 MPC에서는 KKT 시스템이 대칭(준)양의 정부호 형태로 나타날 수 있다. 이 경우 Cholesky 또는 LDLᵀ 분해가 효율적이며, 고정 시간 예산 안에서 해를 구해야 한다.

분해 재사용과 희소화가 만드는 차이

LU 대비 Cholesky의 연산량은 약 50% 수준이며, 1/3 n^3 대 2/3 n^3의 차이와 대칭 저장에 따른 메모리 절감이 있다. 다중 우변 k개를 처리할 때도 분해 1회와 k회 삼각해 풀이로 총 소요시간을 줄일 수 있다.

Newton-Raphson은 근방에 도달한 뒤 2차 수렴 특성으로 반복 횟수가 급감할 수 있다. 부분 피벗팅, 잔차와 조건수 모니터링은 수치 오차를 완화하고 실패를 조기에 감지하는 수단이다. 감쇠와 트러스트 리전은 Newton의 발산률을 줄이고 성공률을 높인다.

희소 분해와 순서화를 적용하면 수십·수백만 변수 규모의 문제까지 다룰 수 있다. 온라인과 배치가 섞인 워크로드에서는 분해 캐시와 재사용이 처리량 증대에 도움이 된다.

구현에서 확인할 코드 경로

감쇠를 포함한 Newton-Raphson

전제: Python 3.10+, NumPy 선택사항(미사용 가능)

def newton(f, df, x0, tol=1e-8, max_iter=50, damping=True):
    x = float(x0)
    fx = f(x)
    for k in range(max_iter):
        dfx = df(x)
        if not (abs(dfx) > 1e-14):
            # 도함수 불량: 미세 스텝 혹은 종료
            raise RuntimeError("Derivative near zero")
        step = -fx / dfx
        # 단순 백트래킹 감쇠
        if damping:
            t = 1.0
            for _ in range(20):
                xn = x + t * step
                fn = f(xn)
                if abs(fn) <= (1 - 1e-4 * t) * abs(fx):
                    x, fx = xn, fn
                    break
                t *= 0.5
            else:
                raise RuntimeError("Line search failed")
        else:
            x = x + step
            fx = f(x)
        if abs(fx) < tol and abs(step) < tol:
            return x
    raise RuntimeError("Not converged")

# 예: x^3 - 2 = 0
root = newton(lambda x: x**3 - 2, lambda x: 3*x**2,  x0=1.0)
print(root)  # 2^(1/3)

LU와 Cholesky로 Ax=b 풀기

전제: Python 3.10+, SciPy >= 1.9

import numpy as np
from scipy.linalg import lu_factor, lu_solve, cho_factor, cho_solve

A = np.array([[4., 1., 0.],
              [1., 3., 1.],
              [0., 1., 2.]])

b = np.array([1., 2., 3.])

# LU (일반 행렬)
lu, piv = lu_factor(A)
x_lu = lu_solve((lu, piv), b)

# Cholesky (SPD 가정)
c, lower = cho_factor(A, lower=True, check_finite=True)
x_chol = cho_solve((c, lower), b)

print(x_lu, x_chol)

Cholesky를 적용하기 전에는 A의 대칭성과 정부호성을 확인해야 하며, SPD가 아니면 예외가 발생할 수 있다. 조건수가 큰 경우에는 잔차와 후행 검증을 수행하고, 필요하면 정규화·스케일링 또는 QR·SVD 같은 대체 분해를 고려한다.

운영에서 놓치기 쉬운 선택 기준

LU에는 부분 피벗팅을 적용하고, 행렬 스케일링으로 성장인자와 라운딩오차를 낮춘다. Newton-Raphson에는 라인서치나 트러스트 리전을 두며, 브래킷팅과 다초점 초기화로 초기값을 준비할 수 있다. 야코비안은 자동미분과 색칠법으로 효율적으로 계산한다.

희소 행렬에서는 AMD·METIS 순서화와 블록 구조 인식으로 fill-in을 줄인다. 정밀도는 배정밀도를 기본으로 두고, 정확도 요구에 따라 분해 FP32와 보정 FP64를 조합하는 혼합정밀 방식을 검토할 수 있다. 극단적인 악조건수나 랭크 결핍에는 QR·SVD가 적합하며, 정규방정식 AᵀA는 빠르지만 수치적으로 취약하다.