강인한 작업병렬 삼각 Sylvester 방정식 해법

초록

본 논문은 LAPACK의 dtrsyl과 같은 스칼라 순차 솔버가 제공하는 오버플로우 방지 기능을 작업‑병렬 수준으로 확장한다. 2×2 이하 블록을 이용한 삼각 Sylvester 방정식에 대해 동적 스케일링을 적용하고, 레벨‑3 BLAS 연산을 중심으로 타일링하여 높은 병렬 효율을 유지한다. 실험 결과, 기존 비강인 솔버와 거의 동일한 실행 속도를 보이며, 오버플로우가 발생할 수 있는 경우에도 안정적으로 해를 구한다.

상세 분석

논문은 먼저 Sylvester 방정식 AX+XB=C 를 Bartels‑Stewart 알고리즘으로 변환해 삼각 형태 ˜A Y+Y ˜B=˜C 로 축소한다는 전제를 명확히 한다. 삼각 형태에서 역방향 대입(backward substitution) 과정은 행렬 원소가 급격히 커져 부동소수점 오버플로우가 발생할 위험이 있다. LAPACK 3.7.0 의 dtrsyl 은 전역 스케일링 팩터 α∈(0,1] 를 도입해 동적으로 ˜C 를 축소함으로써 오버플로우를 방지한다. 그러나 dtrsyl 은 순차적인 레벨‑2 BLAS 연산에 머물러 다중 코어 환경에서 확장성이 제한된다.

저자들은 두 단계의 개선을 제시한다. 첫 번째는 스칼라 수준에서 기존 dtrsyl 에 보호 연산 ProtectUpdate 와 ProtectDivision 을 삽입해 모든 선형 업데이트와 나눗셈을 안전하게 만든 RobustSyl 알고리즘이다. 이 알고리즘은 각 블록(1×1 혹은 2×2) 에 대해 작은 Sylvester 방정식을 Gaussian elimination 으로 풀고, 업데이트 전후에 최대 ∞‑노름을 기준으로 스케일링 팩터 ζ 를 계산한다. ζ는 ‖C‖∞, ‖A‖∞, ‖Y‖∞ 의 합이 오버플로우 한계 Ω 를 초과하지 않도록 보장한다.

두 번째 단계는 위 스칼라 알고리즘을 타일링하여 레벨‑3 BLAS 연산 중심으로 재구성한 RobustTile 알고리즘이다. 행렬을 M×M, N×N 타일로 분할하고, 각 타일마다 로컬 스케일링 팩터 αₖℓ 를 부여한다. 타일 업데이트는 RobustUpdate 서브루틴을 통해 수행되며, 여기서는 세 스케일링 팩터(α,β,γ) 중 최소값 η 를 선택하고, ProtectUpdate 로 추가 스케일링 ζ 를 계산해 최종 스케일링 δ=η·ζ 를 적용한다. 이렇게 하면 타일 간에 일관된 스케일링이 유지되면서도, 각 타일 연산은 대규모 행렬‑행렬 곱(GEMM) 형태로 실행돼 메모리 대역폭과 캐시 활용이 최적화된다.

병렬 실행 모델은 작업(task) 기반 런타임을 가정한다. RobustSyl 은 특정 타일 Yₖℓ 에 대한 읽기‑쓰기 의존성을 갖고, RobustUpdate 은 해당 타일을 독점적으로 수정한다. 의존성 그래프는 위쪽‑아래, 좌우 블록 순서에 따라 자연스럽게 형성되며, OpenMP 혹은 StarPU 같은 런타임이 이를 스케줄링한다. 또한, 행렬 ‖·‖∞ 계산을 사전에 별도 작업으로 수행해 런타임 오버헤드를 최소화한다.

실험에서는 Broadwell Xeon E5‑2690v4 (28코어, 2 NUMA) 환경에서 순차 및 강인/비강인 솔버를 비교했다. µ=m, ν=n 인 경우 스케일링이 전혀 필요 없으며, drsylv(제안 솔버)는 FLA Sylv와 FLASH Sylv에 비해 약 5‑10% 정도 느리지만, 오버플로우 방지 기능을 갖춘 dtrsyl 대비 2배 이상 빠른 성능을 보였다. 강인성을 요구하는 경우(µ,ν 를 작게 설정해 성장률을 높인 경우) drsylv 은 스케일링이 발생해도 효율이 크게 저하되지 않으며, 강인 비강인 솔버 간 성능 격차가 20‑30% 수준에 머문다. 강인 솔버의 정확도는 상대 잔차 ‖R‖F/(‖˜A‖F+‖˜B‖F)‖Y‖F+‖α˜C‖F 로 측정했을 때, 모두 10⁻¹⁴ 이하의 수치적 안정성을 유지한다.

핵심 기여는 (1) 오버플로우 방지를 위한 동적 스케일링을 레벨‑3 BLAS 기반 타일 작업에 적용한 점, (2) 로컬 스케일링 팩터를 이용해 전역 스케일링에 비해 병렬 효율을 크게 향상시킨 점, (3) 기존 비강인 작업‑병렬 Sylvester 솔버와 동일한 스케일링 없이도 거의 동일한 실행 시간을 달성한 점이다. 이 접근법은 Bartels‑Stewart 전체 파이프라인에 바로 통합될 수 있어, 대규모 제어·시스템 시뮬레이션, 신호 처리 등에서 실시간 혹은 고신뢰성이 요구되는 상황에 유용하다.