클렌쇼 점화식

1. 개요

클렌쇼 점화식 Clenshaw recurrence
하는 일	세 항 점화식을 따르는 함수족의 급수 $\sum c_k\phi_k(x)$ 값 계산
연산량	체비쇼프: 항당 곱셈 1~2번 · 덧셈 2번, 저장은 스칼라 2개
대표 쓰임	체비쇼프·르장드르 급수 평가, 괴르첼 알고리즘, 특수함수 라이브러리
약점	$x\to\pm1$ 에서 오차가 차수와 함께 자란다 → 라인슈 수정
이름	클렌쇼 (1955, 체비쇼프 급수)

체비쇼프 다항식을 하나씩 만들어 더하지 마라. 점화식을 뒤에서부터 접어라.

클렌쇼 점화식(Clenshaw recurrence, Clenshaw’s algorithm)은 세 항 점화식 $\phi_{k+1}(x) = \alpha_k(x)\phi_k(x) + \beta_k(x)\phi_{k-1}(x)$ 를 만족하는 함수족으로 쓴 유한 급수

$$S(x) = \sum_{k=0}^{n} c_k\,\phi_k(x)$$

의 값을, 개별 $\phi_k(x)$ 를 하나도 만들지 않고 계수 쪽에서 거꾸로 점화를 돌려 계산하는 방법이다. 영국 국립물리연구소(NPL)의 찰스 클렌쇼가 1955년에 체비쇼프 급수용으로 내놓았고, 같은 구조가 르장드르 다항식·에르미트 다항식 같은 모든 직교다항식, 그리고 삼각급수에도 그대로 통한다.

구조적으로는 호너 방법의 일반화다. 호너는 「단항식 $x^k$ 의 점화식 $x^{k+1} = x\cdot x^k$」를 거꾸로 접은 것이고, 클렌쇼는 그 자리에 항이 하나 더 붙은 세 항 점화식을 넣은 것이다. 그래서 체비쇼프 다항식 기반 근사(Chebfun, 스펙트럴 코드, 수학 함수 라이브러리의 다항식 근사)에서 급수를 값으로 바꾸는 표준 부품이 이것이다.¹

2. 알고리즘

뒤에서부터 보조 수열 $b_k$ 를 만든다($b_{n+1} = b_{n+2} = 0$).

$$b_k = c_k + \alpha_k(x)\,b_{k+1} + \beta_{k+1}(x)\,b_{k+2}, \qquad k = n, n-1, \dots, 1$$

그러면 급수 값은

$$S(x) = c_0\,\phi_0(x) + b_1\,\phi_1(x) + \beta_1(x)\,\phi_0(x)\,b_2$$

다. 증명은 한 줄짜리 망원(telescoping)이다: $c_k = b_k - \alpha_k b_{k+1} - \beta_{k+1}b_{k+2}$ 를 원래 합에 대입하고 $\phi_k$ 끼리 묶으면, 점화식 덕분에 $k \ge 2$ 의 계수가 전부 0 이 되고 처음 두 항만 남는다.

체비쇼프 판. $T_{k+1} = 2xT_k - T_{k-1}$ 이므로 $\alpha_k = 2x$, $\beta_k = -1$ 이고, $T_0 = 1$, $T_1 = x$ 를 넣으면

$$b_k = c_k + 2x\,b_{k+1} - b_{k+2}, \qquad S(x) = c_0 + x\,b_1 - b_2$$

가 된다. 코드로는

b1 = b2 = 0
for k = n … 1:  b1, b2 = 2*x*b1 - b2 + c[k], b1
S = x*b1 - b2 + c[0]

이다. $2x$ 를 미리 한 번 계산해 두면 항당 곱셈 1번, 덧셈 2번이다. $T_k(x)$ 를 점화식으로 하나씩 만들어 $c_kT_k$ 를 누적하는 순진한 방식은 곱셈이 항당 2번이고 변수도 하나 더 든다.

르장드르 판. $(k+1)P_{k+1} = (2k+1)xP_k - kP_{k-1}$ 이므로 $\alpha_k = \frac{(2k+1)x}{k+1}$, $\beta_k = -\frac{k}{k+1}$ 다. 무작위 계수 30개, $x = 0.37$ 로 위 공식을 돌리면 NumPy legval 과 비트 단위로 같은 값(0.01760985060560727)이 나왔다. 주의할 것은 $\beta$ 의 첨자다 — 점화식의 $b_{k+2}$ 에 붙는 것은 $\beta_k$ 가 아니라 $\beta_{k+1}$ 이다. 체비쇼프처럼 $\beta$ 가 상수인 경우엔 티가 안 나다가 르장드르로 옮기는 순간 틀린 값이 나오는 단골 버그다.

도함수. $S'(x)$ 도 같은 틀로 나온다. 점화식을 $x$ 로 미분하면 $b_k' = 2b_{k+1} + 2x\,b_{k+1}' - b_{k+2}'$ (체비쇼프)가 되므로 $b$ 와 $b'$ 를 한 루프에서 같이 돌리면 된다. 호너의 「조립제법 두 번으로 $p$ 와 $p'$」와 같은 발상이고, 체비쇼프 근사 위에서 뉴턴-랩슨법을 돌릴 때 쓴다.

3. 왜 그냥 더하지 않나 — 그리고 삼각급수

순진한 방식이 꼭 틀리는 것은 아니다. 체비쇼프 점화식은 $[-1, 1]$ 위에서 안정이라 $T_k$ 를 앞으로 만들어도 오차가 폭주하지 않는다. 그래도 클렌쇼를 쓰는 이유는 셋이다.

• 싸다. 항당 곱셈이 반이고, 저장이 스칼라 둘이다.
• $\phi_k$ 를 몰라도 된다. 점화 계수 $\alpha_k$, $\beta_k$ 만 있으면 되므로 정규화가 제각각인 직교다항식족을 같은 코드 하나로 처리한다.
• 삼각함수를 안 부른다. $T_k(\cos\theta) = \cos k\theta$ 이니 $\sum c_k\cos(k\theta)$ 를 그대로 계산해도 되지만, $\cos$ 를 $n$ 번 부르는 것은 비싸고 $k\theta$ 의 인자 오차가 $k$ 에 비례해 커진다. 무작위 계수 1000개, $x = 0.5$ 에서 재 보면 코사인 직접 합의 오차(계수 절댓값 합 대비)는 $3.6\times10^{-15}$, 클렌쇼는 $3.2\times10^{-17}$ 이었다.

마지막 점을 거꾸로 쓰면 삼각급수 $\sum c_k\cos k\theta$, $\sum c_k\sin k\theta$ 를 삼각함수 호출 두 번($\cos\theta$, $\sin\theta$)으로 계산하는 법이 된다. $\alpha = 2\cos\theta$, $\beta = -1$ 인 클렌쇼가 곧 신호처리의 괴르첼 알고리즘(Goertzel, 1958)이다 — 이산 푸리에 변환(고속 푸리에 변환 참고)의 주파수 몇 개만 필요할 때(DTMF 전화 다이얼 톤 검출이 고전적 예) FFT 대신 쓴다.

4. 오차 — 끝점 근처의 함정

클렌쇼의 반올림 오차는 「보조 수열 $b_k$ 의 크기」에 비례한다. $x$ 가 내부에 있으면 $b_k$ 가 계수 크기 수준에 머물러 문제가 없다. 그런데 $x \to \pm1$ 이면 점화식 $b_k = 2xb_{k+1} - b_{k+2}$ 의 특성근 $e^{\pm i\theta}$ 가 $\theta \to 0$ 에서 겹근 1 로 합쳐지고, 반올림 오차가 선형으로 쌓인다(겹근이 있는 점화식의 해가 $k$ 에 비례해 자라는 것과 같은 이유). 젠틀먼(1969)이 괴르첼 알고리즘에 대해 이 오차를 분석했고, 같은 결론이 체비쇼프 클렌쇼에 적용된다.

실측(무작위 정규 계수, 참값은 50자리 정밀도 mpmath, 오차는 계수 절댓값 합 $\sum|c_k|$ 로 나눈 값):

$n$	$x$	클렌쇼	순진한 합	코사인 직접 합	라인슈 수정
100	0.5	5.3e-20	4.1e-17	1.4e-16	5.3e-20
100	$1-10^{-8}$	4.3e-15	2.1e-15	4.9e-17	1.9e-17
1000	0.5	3.2e-17	8.8e-17	3.6e-15	3.2e-17
1000	$1-10^{-4}$	1.3e-14	9.2e-16	4.3e-17	3.3e-18
1000	$1-10^{-8}$	6.2e-14	2.3e-14	2.3e-17	4.8e-17

내부에서는 클렌쇼가 가장 좋지만 끝점 근처에선 기계 엡실론의 수백 배까지 나빠지고, $n$ 이 클수록 더하다. 대부분의 응용(유효숫자 13~14자리면 충분한 근사)에서는 무해하지만, 끝점 값이 0 근처인 함수(예: 경계조건을 만족하는 해)를 상대오차로 볼 때는 바로 문제가 된다.

라인슈 수정. 차분 $d_k = b_k - b_{k+1}$ 을 주 변수로 삼아 점화식을 다시 쓰면($x$ 가 1 에 가까울 때)

$$d_k = 2(x-1)\,b_{k+1} + d_{k+1} + c_k, \qquad b_k = d_k + b_{k+1}$$

이 되고, 급수 값은 $S = c_0 + (x-1)b_1 + d_1$ 이다. 작은 수 $x - 1$ 을 정확히(입력에서 직접, 또는 $\theta$ 로부터 $-2\sin^2(\theta/2)$ 로) 넣는 것이 핵심이다 — $2x$ 를 먼저 만들고 나서 2 를 빼면 정보가 이미 반올림으로 날아간 뒤다. $x$ 가 $-1$ 쪽이면 $d_k = b_k + b_{k+1}$ 과 $x+1$ 로 대칭인 판을 쓴다. 표의 마지막 열이 이것으로, 끝점에서도 $10^{-17}$ 수준을 유지한다. 이 판은 슈퇴어·불리르슈의 교과서를 통해 퍼졌고 젠틀먼의 분석과 짝으로 인용된다.²

5. 실무 메모

• 계수 순서. 루프는 높은 차수부터 돈다. 체비쇼프 계수가 빠르게 감쇠하는 매끄러운 함수라면 꼬리 $c_k$ 가 먼저 작은 값으로 들어가고 큰 값이 나중에 더해져, 합산 순서 면에서도 유리하다.
• $c_0$ 의 반값 관례. 체비쇼프 급수를 $\tfrac12 c_0 + \sum_{k\ge1} c_kT_k$ 로 쓰는 책(뉴메리컬 레시피 등)과 $\sum_{k\ge0}$ 로 쓰는 책이 섞여 있다. 마지막 줄의 $c_0$ 를 그대로 쓸지 반으로 쓸지는 출처의 관례를 확인할 것 — 수식은 맞는데 상수항이 정확히 $c_0/2$ 만큼 틀린 버그가 여기서 나온다.
• 구간 변환. $[a, b]$ 위의 근사라면 $x \mapsto (2x - a - b)/(b - a)$ 로 먼저 옮긴다. 이 변환의 반올림이 위 표의 「$x - 1$ 을 정확히」 조건을 깨는 흔한 경로다.
• 행렬 인자. $x$ 자리에 행렬 $A$ 를 넣으면 $\sum c_kT_k(A)$ 의 행렬-벡터 곱 $\sum c_kT_k(A)v$ 를 $A$ 와 벡터 곱만으로 계산한다. 커널 다항식법(KPM)의 상태밀도 계산, 체비쇼프 필터 부분공간 반복, 다항식 전처리기가 모두 이 형태다.
• 클렌쇼-커티스 구적과 혼동 금지. 같은 클렌쇼(그리고 커티스, 1960)의 이름이 붙었지만 이쪽은 적분 규칙이다. 「클렌쇼 알고리즘」이라 하면 점화식, 「클렌쇼-커티스」라 하면 구적이다.³

6. 관련 문서

• 체비쇼프 다항식 · 직교다항식 · 르장드르 다항식 · 에르미트 다항식
• 호너 방법 · 뉴턴 보간 · 바리센트릭 보간 · 보간과 근사
• 클렌쇼-커티스 구적 · 스펙트럴 방법 · 고속 푸리에 변환
• 후진 오차 해석 · 부동소수점 연산 · 오차 없는 변환

7. Footnotes

1. 옥스퍼드의 Chebfun 은 「함수를 체비쇼프 계수 배열로 들고 다니면 부동소수점 수처럼 다룰 수 있다」는 철학의 MATLAB 패키지다. 그 안에서 함수값이 필요할 때마다 돌아가는 게 결국 이 다섯 줄짜리 루프다. 거대한 추상화의 밑바닥엔 늘 이런 루프가 하나 깔려 있다. ↩

2. 그래서 이 수정은 책마다 「라인슈 수정」, 「젠틀먼-라인슈」, 「수정 클렌쇼」 등 이름이 제각각이다. 원 출처를 둘러싼 이름 붙이기 싸움은 수치해석의 국룰이니 너무 신경 쓰지 말 것. ↩

3. 한 사람이 1955년엔 점화식, 1960년엔 구적법으로 이름을 두 번 남겼다. 체비쇼프 다항식 하나로 평생 갈아넣은 결과다. 당시 NPL 에선 체비쇼프 급수로 수학 함수표를 만드는 것이 현업이었다. ↩