Clenshaw Algorithm

Gypsophila
  2026-01-28 2026-01-28 Created   2026-01-29 14:03:02 2026-01-29 14:03:02 Updated    

Clenshaw算法是一种可以高效、稳定地计算在某族正交多项式基底下展开的函数逼近值的递归方法。具体来说，对于一个形式为

的阶多项式展开，其中是一族满足以下三项递归关系的定义在上的正交多项式:

其中,,和是已知函数，我们希望计算某一点处的的值，即以下内积

Clenshaw的关键思想是注意到递推式可以被写为如下的矩阵形式:

上式可被简记为

其中是一个下三对角矩阵，是第一列为，其余元素为的单位向量。使用上述关系，我们可以将写成如下形式:

注意到是一个下三对角矩阵，因此其转置是一个上三对角矩阵，计算只需使用一次回代法，最终的值即为的第一个分量。当都已经被预计算得到之后，一共需要次乘法与次加法即可完成回代得到，总的计算复杂度为flops。

实际上，Clenshaw算法也可以被看作是单项式基下的秦九韶算法（Horner方法）在正交多项式基底下的推广：对于单项式基底，递推关系为

此时上面介绍的Clenshaw算法退化为秦九韶算法：

Benchmark

下图展示了Clenshaw算法在Chebyshev基底下计算阶多项式展开在个点处的值时的耗时曲线。因为Clenshaw算法在计算单点值时的复杂度为，所以理论上在计算个点处的值时总的复杂度为，下图中虚线表示的复杂度曲线，可以看到Clenshaw算法的实际耗时与理论复杂度吻合良好。

Something is Wrong
Benchmark of Clenshaw Algorithm
Code: MATLAB code for benchmarking Clenshaw's algorithm. 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
nn = round(logspace(1,5,50));
j = 1; 
for n = nn 
    c = rand(n,1); 
    x = linspace(-1,1,n)'; 
    s = tic;
    vals = chebtech.clenshaw( x, c ); 
    t_clen(j) = toc(s); 
    j = j + 1; 
    n
end
%%
loglog(nn, t_clen, 'linewidth',2), hold on, 
loglog(nn(30:end), nn(30:end).^2/1e10, 'linewidth',2, 'LineStyle','--', 'Color','black')
xlabel('N')
ylabel('Time')
title("Clenshaw's Algorithm")
text(2*1e4, 0.01, 'O(N^2)', 'FontSize',12)
grid on

Remarks

当仅需要计算展开式在少量点处的值时，Clenshaw算法是最常用且几乎总是高效的数值方法。但在一些特殊情况下，可能存在更快的算法，例如

  1. 在Chebyshev基底下要计算阶逼近在全部的个Chebyshev点处的值时，应当使用FFT-based的快速Chebyshev变换算法，其复杂度为，与之对比Clenshaw算法的复杂度为。
  2. 在Legendre基底下要计算阶逼近在全部的个Gauss-Legendre点处的值时，应当使用专门的快速Legendre变换算法([1][2][3])，其复杂度为，与之对比Clenshaw算法的复杂度为。
  3. 对于更一般的Jacobi多项式基底下的展开，在全部的个Gauss-Jacobi点处计算值时，也存在类似的快速变换算法[4]，但是通常具有较大的常数系数，这使得在通常规模问题下Clenshaw算法仍然更为高效，仅在超大规模问题下才会被快速算法超越，并且此类快速算法的误差一般更大。
  4. 当求值点并非经典的正交多项式节点时，可以尝试NUFFT等方法加速求值过程[5]，例如Faster chebtech restrict中讨论了如何使用NUFFT加速计算Chebyshev展开在等距点上的取值，当阶数时，快速算法才会优于Clenshaw算法，并且附带误差增长为。

See Clenshaw Algorithm on Wikipedia for more details.

Code: MATLAB code for Chebyshev based Clenshaw's algorithm in Chebfun.
clenshaw.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
function y = clenshaw(x, c)
%CLENSHAW   Clenshaw's algorithm for evaluating a Chebyshev expansion.
%   If C is a column vector, Y = CLENSHAW(X, C) evaluates the Chebyshev
%   expansion
%
%     Y = P_N(X) = C(1)*T_0(X) + ... + C(N)*T_{N-1}(X) + C(N+1)*T_N(X)
%
%   using Clenshaw's algorithm.
%
%   If C is an (N+1) x M matrix, then CLENSHAW interprets each of the columns
%   of C as coefficients of a degree N polynomial and evaluates the M Chebyshev
%   expansions
%
%     Y_m = P_N(X) = C(1,m)*T_0(X) + ... + C(N,m)*T_{N-1}(X) + C(N+1,m)*T_N(X)
%
%   for 1 <= m <= M, returning the results as columns of a matrix Y =
%   [Y_1 ... Y_M].
%
%   In both cases, X must be a column vector.
%
% See also CHEBTECH.FEVAL, CHEBTECH.BARY.

% Copyright 2017 by The University of Oxford and The Chebfun Developers. 
% See http://www.chebfun.org/ for Chebfun information.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Developer note: Clenshaw is not typically called directly, but by FEVAL().
%
% Developer note: The algorithm is implemented both for scalar and for vector
% inputs. Of course, the vector implementation could also be used for the scalar
% case, but the additional overheads make it a factor of 2-4 slower. Since the
% code is short, we live with the minor duplication.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% X should be a column vector.
if ( size(x, 2) > 1 )
    warning('CHEBFUN:CHEBTECH:clenshaw:xDim', ...
        'Evaluation points should be a column vector.');
    x = x(:);
end

% Scalar or array-valued function?
if ( size(c, 2) == 1 )
    y = clenshaw_scl(x, c);
else
    y = clenshaw_vec(x, c);
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% function y = clenshaw_scl(x, c)
% % Clenshaw scheme for scalar-valued functions.
% bk1 = 0*x; 
% bk2 = bk1;
% x = 2*x;
% for k = (size(c,1)-1):-1:1
%     bk = c(k) + x.*bk1 - bk2;
%     bk2 = bk1; 
%     bk1 = bk;
% end
% y = c(1) + .5*x.*bk1 - bk2;
% end

% function y = clenshaw_vec(x, c)
% % Clenshaw scheme for array-valued functions.
% x = repmat(x(:), 1, size(c, 2));
% bk1 = zeros(size(x, 1), size(c, 2)); 
% bk2 = bk1;
% e = ones(size(x, 1), 1);
% x = 2*x;
% for k = (size(c,1)-1):-1:1
%     bk = e*c(k,:) + x.*bk1 - bk2;
%     bk2 = bk1; 
%     bk1 = bk;
% end
% y = e*c(1,:) + .5*x.*bk1 - bk2;
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% These code are the same as above with the exception that it does two steps of
% the algorithm at once. This means that the intermediate variables do not need
% to be updated at each step, and can result in a non-trivial speedup. NH 02/14
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function y = clenshaw_scl(x, c)
% Clenshaw scheme for scalar-valued functions.
bk1 = 0*x; 
bk2 = bk1;
x = 2*x;
n = size(c,1)-1;
for k = (n+1):-2:3
    bk2 = c(k) + x.*bk1 - bk2;
    bk1 = c(k-1) + x.*bk2 - bk1;
end
if ( mod(n, 2) )
    tmp = bk1;
    bk1 = c(2) + x.*bk1 - bk2;
    bk2 = tmp;
end
y = c(1) + .5*x.*bk1 - bk2;
end

function y = clenshaw_vec(x, c)
% Clenshaw scheme for array-valued functions.
x = repmat(x(:), 1, size(c, 2));
bk1 = zeros(size(x, 1), size(c, 2)); 
bk2 = bk1;
e = ones(size(x, 1), 1);
x = 2*x;
n = size(c, 1)-1;
for k = (n+1):-2:3
    bk2 = e*c(k,:) + x.*bk1 - bk2;
    bk1 = e*c(k-1,:) + x.*bk2 - bk1;
end
if ( mod(n, 2) )
    tmp = bk1;
    bk1 = e*c(2,:) + x.*bk1 - bk2;
    bk2 = tmp;
end
y = e*c(1,:) + .5*x.*bk1 - bk2;
end

References

  1. Hale N, Townsend A. A fast FFT-based discrete Legendre transform. IMA.NA, 2016, 36(4): 1670-1684.
  2. Hale N, Townsend A. A fast, simple, and stable Chebyshev–Legendre transform using an asymptotic formula. SISC, 2014, 36(1): A148-A167.
  3. Townsend A, Webb M, Olver S. Fast polynomial transforms based on Toeplitz and Hankel matrices. Math.Comp, 2018, 87(312): 1913-1934.
  4. Shen J, Wang Y, Xia J. Fast structured Jacobi-Jacobi transforms. Math.Comp, 2019, 88(318): 1743-1772.
  5. Ruiz-Antolin D, Townsend A. A nonuniform fast Fourier transform based on low rank approximation. SISC, 2018, 40(1): A529-A547.
  • Title: Clenshaw Algorithm
  • Author: Gypsophila
  • Created at : 2026-01-28 00:00:00
  • Updated at : 2026-01-29 14:03:02
  • Link: https://chenx.space/2026/01/28/Clenshaw/
  • License: This work is licensed under CC BY-NC-SA 4.0.
Comments
On this page
Clenshaw Algorithm
  1. Benchmark
  2. Remarks
  3. References