# Fast Analytic Sampling Approximation from Cauchy Kernel

Authors: Youfa Li; Jing Shang; Honglei Yang; Gengrong Zhang; Shouzhi Yang
Journal: Journal of Function Spaces
 (2016)

---

## Abstract

The paper aims at establishing a fast numerical algorithm forB
k
(
f
), where f is any function in the Hardy space H
2
(
T
d
) and k is the scale level. Here, B
k
(
f
) is an approximation to f we recently constructed by applying the multiscale transform to the Cauchy kernel. We establish the matrix expression of B
k
(
f
) and find that it has the structure of a multilevel Hankel matrix. Based on the structure, a fast numerical algorithm is established to compute B
k
(
f
). The computational complexity is given. A numerical experiment is carried out to check the efficiency of our algorithm.


---

## Body

## 1. Introduction

Approximation to a function (a time-continuous signal) by its samples is the heart of modern applied mathematics and engineering and has attracted much attention. Readers are referred to [1–5] for just a few references. For any function f
∈
L
2
(
T
d
) with T
=
[
0,2
π
), it can be reconstructed by its Fourier coefficients as follows: (1)
f
=
∑
k
∈
Z
d
f
^
k
e
i
k
·
.Here, L
2
(
T
d
) is the space of 2
π
Z
d-periodic and square integrable functions, equipped with the inner product (2)
f
,
g
=
1
2
π
d
∫
T
d
f
e
i
t
g
e
i
t
¯
d
t
,
∀
f
,
g
∈
L
2
T
d
,where t
≔
(
t
0
,
…
,
t
d
-
1
)
∈
T
d, e
i
t
≔
(
e
i
t
0
,
…
,
e
i
t
d
-
1
), and d
t
≔
d
t
0
⋯
d
t
d
-
1. Truncating the series in (1) leads to the classical Fourier sampling approximation. Specifically, for sufficiently large N, f can be approximated by its Fourier coefficients {
f
^
(
k
)
:
k
l
2
≤
N
} as follows:(3)
f
≈
∑
k
l
2
≤
N
f
^
k
e
i
k
·
,where ·
l
2 is the Euclid norm of a vector in R
d. From the point of view of feature characterization, (3) can be used to characterize some features such as instantaneous phase and frequency by those of the linear Fourier atoms {
e
i
k
·
:
k
∈
Z
d
}. Note that each linear Fourier atom e
i
k
·
,
k
∈
Z
d has the linear phase and constant frequency. Therefore, e
i
k
· is time-stable. Then, the characterization from (3) implies that, whether f is time-stable or not, its features are characterized by those of the linear Fourier atoms. However, when f is not time-stable, the effectiveness of characterization by the linear atoms in (3) is inferior to that by nonlinear Fourier atoms [6–8].Recently, an interesting sampling approximation from Takenaka-Malmquist (T-M) system, a special class of nonlinear Fourier atoms, has attracted much attention in the literature; for example, see [9, 10] and the references therein. Note that the T-M system is substantially generated from the Schmidt orthogonalization of Cauchy kernel 1
/
(
1
-
z
), where z
∈
D
≔
{
z
∈
C
,
|
z
|
<
1
}
. In this sense, we say that the sampling approximation from the T-M system is substantially from the Cauchy kernel. Implementing the multiscale transform on the tensor product of the Cauchy kernel, Li and Qian [11] constructed the analytic sampling approximation to any function in the Hardy space H
2
(
T
d
), where H
2
(
T
d
), a subspace of L
2
(
T
d
), is defined by (4)
H
2
T
d
≔
f
∈
L
2
T
d
:
f
^
n
=
0
,
∀
n
∈
Z
d
∖
N
0
dwith N
0 being the set of nonnegative integers. Then, using the partial Hilbert transform, [11] generalized the approximation to L
2
(
T
d
). A subsequent unsolved problem of the approximation in [11] is that the fast numerical algorithm has not been established. The aim of this present paper is to solve the problem by the fast computation theory of multilevel circulant matrices.


## 2. Fast Algorithm for Multiscale Analytic Sampling Approximation

### 2.1. Multiscale Analytic Sampling Approximation from Cauchy Kernel

Before introducing the approximation result in [11], we give the definition of the multiscale transform on H
2
(
T
d
). For any k
∈
N
0, define the k-scale transform T
k
m
:
H
2
(
T
d
)
→
H
2
(
T
d
) associated with the shift parameter m
∈
Z
d by (5)
T
k
m
f
=
f
e
i
·
-
2
π
2
-
k
m
,
∀
f
∈
H
2
T
d
.As for the parameter m, we just need to focus on the case of m
∈
L
k due to the periodicity of f, where (6)
L
k
≔
m
0
,
…
,
m
d
-
1
:
m
ν
∈
0,1
,
…
,
2
k
-
1
,
ν
∈
0,1
,
…
,
d
-
1
.Suppose that {
ϕ
k
}
k
=
0
∞ is a sequence of 2
I
d-nonstationary refinement functions in H
2
(
T
d
); namely, they satisfy(7)
ϕ
k
=
∑
m
∈
L
k
P
k
m
T
k
m
ϕ
0
,
∀
k
∈
N
0
,where {
P
k
(
m
)
}
m
∈
L
k belongs to S
(
2
k
), the space of 2
k
Z
d-periodic complex-valued sequences. From the perspective of the Fourier coefficient, it is easy to check that (7) is equivalent to(8)
ϕ
k
^
n
=
P
k
^
n
ϕ
0
^
n
,
∀
n
∈
N
0
d
,where P
k
^
(
n
)
≔
∑
m
∈
L
k
‍
P
k
(
m
)
e
-
i
2
π
2
-
k
n
·
m. Related to ϕ
k, the so-called k-scale projection operator B
k is defined by(9)
B
k
f
=
∑
m
∈
L
k
f
,
T
k
m
ϕ
k
T
k
m
ϕ
k
,
∀
f
∈
H
2
T
d
.Defining ϕ
0
(
e
i
t
) to be ∏
ν
=
0
d
-
1
‍
1
/
1
-
a
¯
ν
e
i
t
ν, the tensor product of the Cauchy kernel on the unit disc D, Li and Qian [11] gave the expression of B
k
(
f
) using the analytic samples of f and estimated the error f
-
B
k
(
f
)
2 concretely.Lemma 1 (see [11]).
Let(10)
ϕ
0
e
i
t
≔
∏
ν
=
0
d
-
1
1
1
-
a
¯
ν
e
i
t
ν
,where e
i
t
=
(
e
i
t
0
,
…
,
e
i
t
d
-
1
), t
=
(
t
0
,
…
,
t
d
-
1
)
∈
T
d, and (
a
0
,
…
,
a
d
-
1
)
∈
D
d. Construct ϕ
k by(11)
ϕ
k
^
m
=
P
k
^
m
ϕ
0
^
m
,
∀
m
∈
N
0
d
,
k
≥
1
,where {
P
k
^
(
m
)
}
m
∈
N
0
d belongs to S
(
2
k
) such that(12)
ϕ
k
^
m
=
2
-
k
/
2
.Then, for any f
∈
H
2
(
T
d
), B
k
(
f
) defined in (9) can be expressed by the k-scale analytic samples {
f
(
a
0
e
i
2
π
2
-
k
m
0
′
′
+
m
0
,
…
,
a
d
-
1
e
i
2
π
2
-
k
(
m
d
-
1
′
′
+
m
d
-
1
)
)
} as follows:(13)
B
k
f
=
∑
m
∈
L
k
∑
m
′
∈
L
k
∑
m
′
′
∈
L
k
P
k
m
′
P
¯
k
m
′
′
f
a
0
e
i
2
π
2
-
k
m
0
′
′
+
m
0
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
′
′
+
m
d
-
1
ϕ
0
e
i
·
-
2
π
2
-
k
m
+
m
′
,where {
P
k
(
m
)
}
m
∈
L
k is the inverse Fourier transform of {
P
k
^
(
m
)
}
m
∈
L
k. Moreover,(14)
f
-
B
k
f
2
2
≤
2
∑
m
∈
L
k
∑
q
∈
N
0
d
∖
0
f
^
m
+
2
k
q
2
+
9
f
2
2
+
4
f
^
0
2
∑
ν
=
0
d
-
1
a
ν
2
k
2
∏
ν
=
0
d
-
1
1
-
a
ν
2
k
2
,where 0
≔
(
0
,
…
,
0
).


### 2.2. Mathematical Materials on Multilevel Circulant Matrix

Following [12–14], we will define a multilevel circulant matrix. We begin with some denotations. For any n
∈
N and n
≔
(
n
0
,
…
,
n
d
-
1
)
∈
N
d, let N
n
≔
{
0,1
,
…
,
n
-
1
} and N
n
≔
N
n
0
×
⋯
×
N
n
d
-
1, where × is the Cartesian product. An n
×
n matrix is referred to as a 1-level matrix of order n. Recursively, matrix C
n is a d-level matrix of order n if it consists of n
0
2 matrices of (
d
-
1
)-level of order (
n
1
,
…
,
n
d
-
1
). To point to the entries of C
n, denote C
n by [
c
i
,
j
]
i
,
j
∈
N
n. A d-level matrix C
n
=
[
c
i
,
j
]
i
,
j
∈
N
n of order n is referred to as a d-level circulant matrix if(15)
c
i
,
j
=
c
i
0
-
j
0
m
o
d
n
0
,
…
,
i
d
-
1
-
j
d
-
1
m
o
d
n
d
-
1
,
∀
i
,
j
∈
N
n
,where i
≔
(
i
0
,
…
,
i
d
-
1
) and j
≔
(
j
0
,
…
,
j
d
-
1
).Ann
×
1 column (or 1
×
n row) vector is referred to as a 1-level vector of order n. Following the definition of a multilevel matrix, vector V
n is a d-level vector of order n if it consists of n
0 vectors of (
d
-
1
)-level of order (
n
1
,
…
,
n
d
-
1
). To point to the entries of V
n, denote V
n by (
V
n
(
i
)
)
i
∈
N
n
.A multilevel circulant matrix has a nice structure as shown in the following lemma.Lemma 2 (see [12–14]).
Suppose thatC
n is a d-level matrix of order n and c
n is the transpose of the first row of C
n, where n
=
(
n
0
,
…
,
n
d
-
1
). Then, C
n is a d-level circulant matrix if and only if(16)
C
n
=
1
Π
n
Φ
n
∗
diag
⁡
Φ
n
∗
c
n
Φ
n
,where Π
n
≔
∏
ν
=
0
d
-
1
‍
n
ν, Φ
n
≔
F
n
0
⊗
⋯
⊗
F
n
d
-
1, and Φ
n
∗ is the transpose conjugate of Φ
n. Here, F
n
ν is the Fourier matrix of order n
ν and ⊗ is the Kronecker product of matrices.Following (15), a d-level matrix H
n
≔
[
h
i
,
j
]
i
,
j
∈
N
n of order n is referred to as a d-level Hankel matrix if (17)
h
i
,
j
=
h
i
0
+
j
0
,
…
,
i
d
-
1
+
j
d
-
1
,
∀
i
,
j
∈
N
n
.Since it can be converted to a d-level circulant matrix by a transform to be given in Lemma 3, the Hankel matrix with the additional property(18)
h
i
,
j
=
h
i
0
+
j
0
m
o
d
n
0
,
…
,
i
d
-
1
+
j
d
-
1
m
o
d
n
d
-
1is crucial for establishing the fast algorithm in Section 2.3.Lemma 3.
Define ad-level matrix of order n by (19)
D
n
≔
E
n
0
⊗
E
n
1
⊗
⋯
⊗
E
n
d
-
1
,where (20)
E
n
v
≔
0
⋯
0
0
1
0
⋯
0
1
0
0
⋯
1
0
0
⋮
⋰
⋮
⋮
⋮
1
⋯
0
0
0
n
v
×
n
v
.If H
n is a d-level Hankel matrix of order n satisfying (18), then(21)
H
~
n
≔
H
n
D
n
=
h
~
i
,
j
i
,
j
∈
N
nis a d-level circulant matrix.Proof.
From (21) and (18), we arrive at(22)
h
~
i
,
j
=
h
i
,
n
-
1
-
j
=
h
i
0
-
1
-
j
0
m
o
d
n
0
,
…
,
i
d
-
1
-
1
-
j
d
-
1
m
o
d
n
d
-
1
,where 1
≔
(
1
,
…
,
1
). By (22), if i
′
≔
(
i
0
′
,
…
,
i
d
-
1
′
) and j
′
≔
(
j
0
′
,
…
,
j
d
-
1
′
) satisfy(23)
i
0
-
1
-
j
0
m
o
d
n
0
=
i
0
′
-
1
-
j
0
′
m
o
d
n
0
,
…
,
i
d
-
1
-
1
-
j
d
-
1
m
o
d
n
0
=
i
d
-
1
′
-
1
-
j
d
-
1
′
m
o
d
n
d
-
1
,then(24)
h
~
i
,
j
=
h
i
,
n
-
1
-
j
=
h
i
′
,
n
-
1
-
j
′
=
h
~
i
′
,
j
′
.Recall that (23) is equivalent to(25)
i
0
-
j
0
m
o
d
n
0
=
i
0
′
-
j
0
′
m
o
d
n
0
,
…
,
i
d
-
1
-
j
d
-
1
m
o
d
n
0
=
i
d
-
1
′
-
j
d
-
1
′
m
o
d
n
d
-
1
,which together with (24) leads to (26)
h
~
i
,
j
=
h
~
i
0
-
j
0
m
o
d
n
0
,
…
,
i
d
-
1
-
j
d
-
1
m
o
d
n
d
-
1
.The proof is concluded.Note 1.
It is clear thatD
n
-
1
=
D
n. Then, it follows from (21) that a d-level Hankel matrix H
n satisfying (18) can be written as(27)
H
n
=
1
Π
n
Φ
n
∗
diag
⁡
Φ
n
∗
h
~
n
Φ
n
D
n
,where h
~
n is the transpose of the first row of H
~
n given in (21).


### 2.3. Fast Multiscale Analytic Sampling Approximation

This subsection aims at developing a fast algorithm to compute the numerical values ofB
k
(
f
) at e
i
t
m
,
x, where(28)
t
m
,
x
≔
x
+
2
π
2
-
k
m
=
x
0
+
2
π
2
-
k
m
0
,
…
,
x
d
-
1
+
2
π
2
-
k
m
d
-
1
,
x
≔
(
x
0
,
…
,
x
d
-
1
) is any fixed point on R
d, and m
=
(
m
0
,
…
,
m
d
-
1
)
∈
L
k. For convenient narration, define a column vector a
k
,
x by(29)
a
k
,
x
≔
B
k
f
e
i
t
m
,
x
m
∈
L
kand a row vector 2
k by 2
k
≔
(
2
k
,
…
,
2
k
)
.Lemma 4.
Using the matrix notation,B
k
(
f
) in (13) can be expressed by(30)
B
k
f
=
B
2
k
P
2
k
T
A
2
k
P
¯
2
k
,where A
2
k is a matrix defined by (31)
A
2
k
≔
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
+
m
d
-
1
′
m
,
m
′
∈
L
k
,
P
2
k
≔
(
P
k
,
m
′
)
m
′
∈
L
k is a column vector, and B
2
k is a matrix of functions given by(32)
B
2
k
≔
ϕ
0
e
i
·
-
2
π
2
-
k
m
+
m
′
m
,
m
′
∈
L
k
.Moreover, A
2
k and B
2
k are both d-level Hankel matrices with property (18).Proof.
For anym
∈
L
k, it follows from (7) that(33)
T
k
m
ϕ
k
=
T
k
m
∑
m
′
∈
L
k
P
k
,
m
′
T
k
m
′
ϕ
0
=
∑
m
′
∈
L
k
P
k
,
m
′
T
k
m
+
m
′
ϕ
0
=
∑
m
′
∈
L
k
P
k
,
m
′
ϕ
0
e
i
·
-
2
π
2
-
k
m
+
m
′
.Then, the column vector (
T
k
m
ϕ
k
)
m
∈
L
k can be expressed by(34)
T
k
m
ϕ
k
m
∈
L
k
=
B
2
k
P
2
k
.By (33) and the Cauchy integral formula, for any m
∈
L
k, we have (35)
f
,
T
k
m
ϕ
k
=
f
,
∑
m
′
∈
L
k
P
k
,
m
′
ϕ
0
e
i
·
-
2
π
2
-
k
m
+
m
′
=
1
2
π
d
∫
T
d
f
e
i
t
∑
m
′
∈
L
k
P
k
,
m
′
ϕ
0
e
i
t
-
2
π
2
-
k
m
+
m
′
¯
d
t
0
⋯
d
t
d
-
1
=
∑
m
′
∈
L
k
P
k
,
m
′
¯
1
2
π
i
d
∫
T
d
f
e
i
t
∏
ν
=
0
d
-
1
1
e
i
t
ν
-
a
ν
e
i
2
π
2
-
k
m
ν
+
m
ν
′
d
e
i
t
0
⋯
d
e
i
t
d
-
1
=
∑
m
′
∈
L
k
P
k
,
m
′
¯
1
2
π
i
d
-
1
∫
T
d
-
1
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
,
e
i
t
1
,
…
,
e
i
t
d
-
1
∏
ν
=
1
d
-
1
1
e
i
t
ν
-
a
ν
e
i
2
π
2
-
k
m
ν
+
m
ν
′
d
e
i
t
1
⋯
d
e
i
t
d
-
1
⋮
=
∑
m
′
∈
L
k
P
k
,
m
′
¯
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
+
m
d
-
1
′
.Therefore, the column vector (
〈
f
,
T
k
m
ϕ
k
〉
)
m
∈
L
k can be rewritten as(36)
f
,
T
k
m
ϕ
k
m
∈
L
k
=
A
2
k
P
¯
2
k
.Now, the proof of (30) is concluded by (9), (34), and (36).
For anym
=
(
m
0
,
m
1
,
…
,
m
d
-
1
) and m
′
=
(
m
0
′
,
m
1
′
,
…
,
m
d
-
1
′
)
∈
L
k, it is straightforward to check that (37)
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
+
m
d
-
1
′
=
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
m
o
d
2
k
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
+
m
d
-
1
′
m
o
d
2
k
,and then A
2
k is a d-level Hankel matrix satisfying (18). Similarly, B
2
k is also a d-level Hankel matrix with property (18).The following lemma is crucial for investigating the matrix expression ofa
k
,
x.Lemma 5.
For anyk
∈
N
0, define a matrix of functions(38)
K
2
k
,
a
≔
1
1
-
a
¯
e
i
·
-
2
π
2
-
k
ξ
+
η
ξ
,
η
∈
N
2
k
,where a
∈
D. Then, for any x
∈
R and m
∈
N
2
k, the value of K
2
k
,
a at x
+
2
π
2
-
k
m satisfies(39)
K
2
k
,
a
x
+
2
π
2
-
k
m
=
K
2
k
,
a
x
S
2
k
m
,where (40)
S
2
k
≔
O
2
k
-
1,1
I
2
k
-
1
1
O
1
,
2
k
-
1
.Proof.
We first prove by induction onm that(41)
S
2
k
m
=
O
2
k
-
m
,
m
I
2
k
-
m
I
m
O
m
,
2
k
-
m
,
m
∈
N
2
k
.Suppose (41) holds with m being replaced by m
-
1. Then, (42)
S
2
k
m
=
S
2
k
S
2
k
m
-
1
=
O
2
k
-
1,1
I
2
k
-
1
1
O
1
,
2
k
-
1
O
2
k
-
m
+
1
,
m
-
1
I
2
k
-
m
+
1
I
m
-
1
O
m
-
1
,
2
k
-
m
+
1
=
O
2
k
-
1,1
I
2
k
-
1
1
O
1
,
2
k
-
1
O
1
,
m
-
1
1
O
1
,
2
k
-
m
O
2
k
-
m
,
m
-
1
O
2
k
-
m
,
1
I
2
k
-
m
I
m
-
1
O
m
-
1,1
O
m
-
1
,
2
k
-
m
=
O
2
k
-
m
,
m
-
1
O
2
k
-
m
,
1
I
2
k
-
m
I
m
-
1
O
m
-
1,1
O
m
-
1
,
2
k
-
m
O
1
,
m
-
1
1
O
1
,
2
k
-
m
=
O
2
k
-
m
,
m
I
2
k
-
m
I
m
O
m
,
2
k
-
m
.Therefore, (41) holds for any m
∈
N
2
k.
For anyξ
,
η
∈
N
2
k, we derive from (41) that (43)
K
2
k
,
a
x
S
2
k
m
ξ
,
η
=
K
2
k
,
a
x
ξ
,
2
k
-
m
+
η
,
0
≤
η
≤
m
-
1
;
K
2
k
,
a
x
ξ
,
η
-
m
,
m
-
1
<
η
≤
2
k
-
1
;
=
1
1
-
a
¯
e
i
x
+
2
π
2
-
k
m
-
2
π
2
-
k
ξ
+
η
=
K
2
k
,
a
x
+
2
π
2
-
k
m
ξ
,
η
.The proof of (39) is concluded.Using (39), we will give the matrix expression of a
k
,
x defined in (29).Theorem 6.
Let{
A
2
k
,
B
2
k
,
P
2
k
} and S
2
k be as in Lemmas 4 and 5, respectively. Define a d-level matrix G
2
k
≔
(
g
m
′
,
m
)
m
′
,
m
∈
L
k by(44)
g
m
′
,
m
=
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
P
2
k
m
′
,where m
=
(
m
0
,
…
,
m
d
-
1
). Then, a
k
,
x in (29) can be expressed by(45)
a
k
,
x
=
G
2
k
B
2
k
t
0
,
x
A
2
k
P
¯
2
k
.Proof.
Using (41) to directly compute (44) gives us that (46)
g
m
′
,
m
=
P
2
k
m
0
+
m
0
′
m
o
d
2
k
,
…
,
m
d
-
1
+
m
d
-
1
′
m
o
d
2
k
,from which we arrive at (47)
g
m
′
,
m
=
g
m
0
+
m
0
′
m
o
d
2
k
,
…
,
m
d
-
1
+
m
d
-
1
′
m
o
d
2
k
,where m
′
=
(
m
0
′
,
…
,
m
d
-
1
′
). Hence, G
2
k is a d-level Hankel matrix with additional property (18).
It follows from (10), (32), (38), and (39) that the value of B
2
k at t
m
,
x defined in (28) can be computed by(48)
B
2
k
t
m
,
x
=
K
2
k
,
a
0
x
0
+
2
π
2
-
k
m
0
⊗
⋯
⊗
K
2
k
,
a
d
-
1
x
0
+
2
π
2
-
k
m
d
-
1
=
K
2
k
,
a
0
x
0
S
2
k
m
0
⊗
⋯
⊗
K
2
k
,
a
d
-
1
x
0
S
2
k
m
d
-
1
=
K
2
k
,
a
0
x
0
⊗
⋯
⊗
K
2
k
,
a
d
-
1
x
0
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
=
B
2
k
t
0
,
x
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
.It is deduced from (30) and (48) that(49)
B
k
f
t
m
,
x
=
B
2
k
t
m
,
x
P
2
k
T
A
2
k
P
¯
2
k
=
B
2
k
t
0
,
x
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
P
2
k
T
A
2
k
P
¯
2
k
=
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
P
2
k
T
B
2
k
T
t
0
,
x
A
2
k
P
¯
2
k
.By (49), the column vector a
k
,
x defined in (29) can be expressed by (50)
a
k
,
x
=
G
2
k
T
B
2
k
T
t
0
,
x
A
2
k
P
¯
2
k
.Since both B
2
k
(
t
0
,
x
) and G
2
k are d-level Hankel matrices, they are both symmetric. Hence, the proof of (45) is concluded.Note 2.
SinceG
2
k, B
2
k
(
t
0
,
x
), and A
2
k are all d-level Hankel matrices with property (18), by (27) and (45), a
k
,
x can be factorized as(51)
a
k
,
x
=
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
g
~
2
k
Φ
2
k
D
2
k
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
b
~
2
k
Φ
2
k
D
2
k
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
a
~
2
k
Φ
2
k
D
2
k
P
¯
2
k
,where g
~
2
k, b
~
2
k, and a
~
2
k are the transposes of the first rows of G
~
2
k
≔
G
2
k
D
2
k, B
~
2
k
(
t
0
,
x
)
≔
B
2
k
(
t
0
,
x
)
D
2
k, and A
~
2
k
≔
A
2
k
D
2
k, respectively.Based on Note2, a fast algorithm for a
k
,
x will be established as follows.Algorithm 7.
Letϕ
0, {
P
k
,
m
}
m
∈
L
k, and {
P
^
k
,
m
}
m
∈
L
k be as in Lemma 1. According to (51), a
k
,
x defined in (29) can be computed by the following six steps.(1)
Byϕ
0
^
(
m
)
=
∏
ν
=
0
d
-
1
‍
a
¯
ν
m
ν, (11), and (12), compute {
P
^
k
,
m
}
m
∈
L
k, where m
=
(
m
0
,
…
,
m
d
-
1
).
(2)
By implementing IFFT on{
P
^
k
,
m
}
m
∈
L
k, we compute P
2
k.
(3)
SinceΦ
2
k
=
F
2
k
⊗
⋯
⊗
F
2
k and F
2
k is the Fourier matrix of order 2
k, Φ
2
k
∗
g
~
2
k, Φ
2
k
∗
a
~
2
k, and Φ
2
k
∗
b
~
2
k are computed by IFFTs.
(4)
Consider(52)
P
2
k
→
D
2
k
D
2
k
P
¯
2
k
→
F
F
T
Φ
2
k
D
2
k
P
¯
2
k
→
diag
⁡
Φ
2
k
∗
a
~
2
k
diag
⁡
Φ
2
k
∗
a
~
2
k
Φ
2
k
D
2
k
P
¯
2
k
→
I
F
F
T
a
k
,
x
1
≔
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
a
~
2
k
Φ
2
k
D
2
k
P
¯
2
k
.
(5)
Consider(53)
a
k
,
x
1
→
D
2
k
D
2
k
a
k
,
x
1
→
F
F
T
Φ
2
k
D
2
k
a
k
,
x
1
→
diag
⁡
Φ
2
k
∗
b
~
2
k
diag
⁡
Φ
2
k
∗
b
~
2
k
Φ
2
k
D
2
k
a
k
,
x
1
→
I
F
F
T
a
k
,
x
2
≔
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
b
~
2
k
Φ
2
k
D
2
k
a
k
,
x
1
.
(6)
Consider(54)
a
k
,
x
2
→
D
2
k
D
2
k
a
k
,
x
2
→
F
F
T
Φ
2
k
D
2
k
a
k
,
x
2
→
diag
⁡
Φ
2
k
∗
b
~
2
k
diag
⁡
Φ
2
k
∗
g
~
2
k
Φ
2
k
D
2
k
a
k
,
x
2
→
I
F
F
T
a
k
,
x
=
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
g
~
2
k
Φ
2
k
D
2
k
a
k
,
x
2
.
Computational Complexity. It is easy to check that directly computing a
k
,
x through (13) costs O
(
2
3
d
k
) operations. In Algorithm 7, however, FFT and IFFT are used for three and seven times, respectively, which cost O
(
2
k
d
·
d
k
) operations. Meanwhile, the complexity of other operations such as the multiplication of D
2
k and P
¯
2
k is O
(
2
k
d
). Therefore, the computational complexity of Algorithm 7 is O
(
2
k
d
·
d
k
).Note 3.
So far, on computinga
k
,
x, we have not found any numerical algorithm better than Algorithm 7. On the other hand, we notice that Algorithm 7 holds for equally spaced points on T
d. In another occasion, we will establish the fast algorithm for the case of nonequally spaced points.


## 2.1. Multiscale Analytic Sampling Approximation from Cauchy Kernel

Before introducing the approximation result in [11], we give the definition of the multiscale transform on H
2
(
T
d
). For any k
∈
N
0, define the k-scale transform T
k
m
:
H
2
(
T
d
)
→
H
2
(
T
d
) associated with the shift parameter m
∈
Z
d by (5)
T
k
m
f
=
f
e
i
·
-
2
π
2
-
k
m
,
∀
f
∈
H
2
T
d
.As for the parameter m, we just need to focus on the case of m
∈
L
k due to the periodicity of f, where (6)
L
k
≔
m
0
,
…
,
m
d
-
1
:
m
ν
∈
0,1
,
…
,
2
k
-
1
,
ν
∈
0,1
,
…
,
d
-
1
.Suppose that {
ϕ
k
}
k
=
0
∞ is a sequence of 2
I
d-nonstationary refinement functions in H
2
(
T
d
); namely, they satisfy(7)
ϕ
k
=
∑
m
∈
L
k
P
k
m
T
k
m
ϕ
0
,
∀
k
∈
N
0
,where {
P
k
(
m
)
}
m
∈
L
k belongs to S
(
2
k
), the space of 2
k
Z
d-periodic complex-valued sequences. From the perspective of the Fourier coefficient, it is easy to check that (7) is equivalent to(8)
ϕ
k
^
n
=
P
k
^
n
ϕ
0
^
n
,
∀
n
∈
N
0
d
,where P
k
^
(
n
)
≔
∑
m
∈
L
k
‍
P
k
(
m
)
e
-
i
2
π
2
-
k
n
·
m. Related to ϕ
k, the so-called k-scale projection operator B
k is defined by(9)
B
k
f
=
∑
m
∈
L
k
f
,
T
k
m
ϕ
k
T
k
m
ϕ
k
,
∀
f
∈
H
2
T
d
.Defining ϕ
0
(
e
i
t
) to be ∏
ν
=
0
d
-
1
‍
1
/
1
-
a
¯
ν
e
i
t
ν, the tensor product of the Cauchy kernel on the unit disc D, Li and Qian [11] gave the expression of B
k
(
f
) using the analytic samples of f and estimated the error f
-
B
k
(
f
)
2 concretely.Lemma 1 (see [11]).
Let(10)
ϕ
0
e
i
t
≔
∏
ν
=
0
d
-
1
1
1
-
a
¯
ν
e
i
t
ν
,where e
i
t
=
(
e
i
t
0
,
…
,
e
i
t
d
-
1
), t
=
(
t
0
,
…
,
t
d
-
1
)
∈
T
d, and (
a
0
,
…
,
a
d
-
1
)
∈
D
d. Construct ϕ
k by(11)
ϕ
k
^
m
=
P
k
^
m
ϕ
0
^
m
,
∀
m
∈
N
0
d
,
k
≥
1
,where {
P
k
^
(
m
)
}
m
∈
N
0
d belongs to S
(
2
k
) such that(12)
ϕ
k
^
m
=
2
-
k
/
2
.Then, for any f
∈
H
2
(
T
d
), B
k
(
f
) defined in (9) can be expressed by the k-scale analytic samples {
f
(
a
0
e
i
2
π
2
-
k
m
0
′
′
+
m
0
,
…
,
a
d
-
1
e
i
2
π
2
-
k
(
m
d
-
1
′
′
+
m
d
-
1
)
)
} as follows:(13)
B
k
f
=
∑
m
∈
L
k
∑
m
′
∈
L
k
∑
m
′
′
∈
L
k
P
k
m
′
P
¯
k
m
′
′
f
a
0
e
i
2
π
2
-
k
m
0
′
′
+
m
0
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
′
′
+
m
d
-
1
ϕ
0
e
i
·
-
2
π
2
-
k
m
+
m
′
,where {
P
k
(
m
)
}
m
∈
L
k is the inverse Fourier transform of {
P
k
^
(
m
)
}
m
∈
L
k. Moreover,(14)
f
-
B
k
f
2
2
≤
2
∑
m
∈
L
k
∑
q
∈
N
0
d
∖
0
f
^
m
+
2
k
q
2
+
9
f
2
2
+
4
f
^
0
2
∑
ν
=
0
d
-
1
a
ν
2
k
2
∏
ν
=
0
d
-
1
1
-
a
ν
2
k
2
,where 0
≔
(
0
,
…
,
0
).


## 2.2. Mathematical Materials on Multilevel Circulant Matrix

Following [12–14], we will define a multilevel circulant matrix. We begin with some denotations. For any n
∈
N and n
≔
(
n
0
,
…
,
n
d
-
1
)
∈
N
d, let N
n
≔
{
0,1
,
…
,
n
-
1
} and N
n
≔
N
n
0
×
⋯
×
N
n
d
-
1, where × is the Cartesian product. An n
×
n matrix is referred to as a 1-level matrix of order n. Recursively, matrix C
n is a d-level matrix of order n if it consists of n
0
2 matrices of (
d
-
1
)-level of order (
n
1
,
…
,
n
d
-
1
). To point to the entries of C
n, denote C
n by [
c
i
,
j
]
i
,
j
∈
N
n. A d-level matrix C
n
=
[
c
i
,
j
]
i
,
j
∈
N
n of order n is referred to as a d-level circulant matrix if(15)
c
i
,
j
=
c
i
0
-
j
0
m
o
d
n
0
,
…
,
i
d
-
1
-
j
d
-
1
m
o
d
n
d
-
1
,
∀
i
,
j
∈
N
n
,where i
≔
(
i
0
,
…
,
i
d
-
1
) and j
≔
(
j
0
,
…
,
j
d
-
1
).Ann
×
1 column (or 1
×
n row) vector is referred to as a 1-level vector of order n. Following the definition of a multilevel matrix, vector V
n is a d-level vector of order n if it consists of n
0 vectors of (
d
-
1
)-level of order (
n
1
,
…
,
n
d
-
1
). To point to the entries of V
n, denote V
n by (
V
n
(
i
)
)
i
∈
N
n
.A multilevel circulant matrix has a nice structure as shown in the following lemma.Lemma 2 (see [12–14]).
Suppose thatC
n is a d-level matrix of order n and c
n is the transpose of the first row of C
n, where n
=
(
n
0
,
…
,
n
d
-
1
). Then, C
n is a d-level circulant matrix if and only if(16)
C
n
=
1
Π
n
Φ
n
∗
diag
⁡
Φ
n
∗
c
n
Φ
n
,where Π
n
≔
∏
ν
=
0
d
-
1
‍
n
ν, Φ
n
≔
F
n
0
⊗
⋯
⊗
F
n
d
-
1, and Φ
n
∗ is the transpose conjugate of Φ
n. Here, F
n
ν is the Fourier matrix of order n
ν and ⊗ is the Kronecker product of matrices.Following (15), a d-level matrix H
n
≔
[
h
i
,
j
]
i
,
j
∈
N
n of order n is referred to as a d-level Hankel matrix if (17)
h
i
,
j
=
h
i
0
+
j
0
,
…
,
i
d
-
1
+
j
d
-
1
,
∀
i
,
j
∈
N
n
.Since it can be converted to a d-level circulant matrix by a transform to be given in Lemma 3, the Hankel matrix with the additional property(18)
h
i
,
j
=
h
i
0
+
j
0
m
o
d
n
0
,
…
,
i
d
-
1
+
j
d
-
1
m
o
d
n
d
-
1is crucial for establishing the fast algorithm in Section 2.3.Lemma 3.
Define ad-level matrix of order n by (19)
D
n
≔
E
n
0
⊗
E
n
1
⊗
⋯
⊗
E
n
d
-
1
,where (20)
E
n
v
≔
0
⋯
0
0
1
0
⋯
0
1
0
0
⋯
1
0
0
⋮
⋰
⋮
⋮
⋮
1
⋯
0
0
0
n
v
×
n
v
.If H
n is a d-level Hankel matrix of order n satisfying (18), then(21)
H
~
n
≔
H
n
D
n
=
h
~
i
,
j
i
,
j
∈
N
nis a d-level circulant matrix.Proof.
From (21) and (18), we arrive at(22)
h
~
i
,
j
=
h
i
,
n
-
1
-
j
=
h
i
0
-
1
-
j
0
m
o
d
n
0
,
…
,
i
d
-
1
-
1
-
j
d
-
1
m
o
d
n
d
-
1
,where 1
≔
(
1
,
…
,
1
). By (22), if i
′
≔
(
i
0
′
,
…
,
i
d
-
1
′
) and j
′
≔
(
j
0
′
,
…
,
j
d
-
1
′
) satisfy(23)
i
0
-
1
-
j
0
m
o
d
n
0
=
i
0
′
-
1
-
j
0
′
m
o
d
n
0
,
…
,
i
d
-
1
-
1
-
j
d
-
1
m
o
d
n
0
=
i
d
-
1
′
-
1
-
j
d
-
1
′
m
o
d
n
d
-
1
,then(24)
h
~
i
,
j
=
h
i
,
n
-
1
-
j
=
h
i
′
,
n
-
1
-
j
′
=
h
~
i
′
,
j
′
.Recall that (23) is equivalent to(25)
i
0
-
j
0
m
o
d
n
0
=
i
0
′
-
j
0
′
m
o
d
n
0
,
…
,
i
d
-
1
-
j
d
-
1
m
o
d
n
0
=
i
d
-
1
′
-
j
d
-
1
′
m
o
d
n
d
-
1
,which together with (24) leads to (26)
h
~
i
,
j
=
h
~
i
0
-
j
0
m
o
d
n
0
,
…
,
i
d
-
1
-
j
d
-
1
m
o
d
n
d
-
1
.The proof is concluded.Note 1.
It is clear thatD
n
-
1
=
D
n. Then, it follows from (21) that a d-level Hankel matrix H
n satisfying (18) can be written as(27)
H
n
=
1
Π
n
Φ
n
∗
diag
⁡
Φ
n
∗
h
~
n
Φ
n
D
n
,where h
~
n is the transpose of the first row of H
~
n given in (21).


## 2.3. Fast Multiscale Analytic Sampling Approximation

This subsection aims at developing a fast algorithm to compute the numerical values ofB
k
(
f
) at e
i
t
m
,
x, where(28)
t
m
,
x
≔
x
+
2
π
2
-
k
m
=
x
0
+
2
π
2
-
k
m
0
,
…
,
x
d
-
1
+
2
π
2
-
k
m
d
-
1
,
x
≔
(
x
0
,
…
,
x
d
-
1
) is any fixed point on R
d, and m
=
(
m
0
,
…
,
m
d
-
1
)
∈
L
k. For convenient narration, define a column vector a
k
,
x by(29)
a
k
,
x
≔
B
k
f
e
i
t
m
,
x
m
∈
L
kand a row vector 2
k by 2
k
≔
(
2
k
,
…
,
2
k
)
.Lemma 4.
Using the matrix notation,B
k
(
f
) in (13) can be expressed by(30)
B
k
f
=
B
2
k
P
2
k
T
A
2
k
P
¯
2
k
,where A
2
k is a matrix defined by (31)
A
2
k
≔
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
+
m
d
-
1
′
m
,
m
′
∈
L
k
,
P
2
k
≔
(
P
k
,
m
′
)
m
′
∈
L
k is a column vector, and B
2
k is a matrix of functions given by(32)
B
2
k
≔
ϕ
0
e
i
·
-
2
π
2
-
k
m
+
m
′
m
,
m
′
∈
L
k
.Moreover, A
2
k and B
2
k are both d-level Hankel matrices with property (18).Proof.
For anym
∈
L
k, it follows from (7) that(33)
T
k
m
ϕ
k
=
T
k
m
∑
m
′
∈
L
k
P
k
,
m
′
T
k
m
′
ϕ
0
=
∑
m
′
∈
L
k
P
k
,
m
′
T
k
m
+
m
′
ϕ
0
=
∑
m
′
∈
L
k
P
k
,
m
′
ϕ
0
e
i
·
-
2
π
2
-
k
m
+
m
′
.Then, the column vector (
T
k
m
ϕ
k
)
m
∈
L
k can be expressed by(34)
T
k
m
ϕ
k
m
∈
L
k
=
B
2
k
P
2
k
.By (33) and the Cauchy integral formula, for any m
∈
L
k, we have (35)
f
,
T
k
m
ϕ
k
=
f
,
∑
m
′
∈
L
k
P
k
,
m
′
ϕ
0
e
i
·
-
2
π
2
-
k
m
+
m
′
=
1
2
π
d
∫
T
d
f
e
i
t
∑
m
′
∈
L
k
P
k
,
m
′
ϕ
0
e
i
t
-
2
π
2
-
k
m
+
m
′
¯
d
t
0
⋯
d
t
d
-
1
=
∑
m
′
∈
L
k
P
k
,
m
′
¯
1
2
π
i
d
∫
T
d
f
e
i
t
∏
ν
=
0
d
-
1
1
e
i
t
ν
-
a
ν
e
i
2
π
2
-
k
m
ν
+
m
ν
′
d
e
i
t
0
⋯
d
e
i
t
d
-
1
=
∑
m
′
∈
L
k
P
k
,
m
′
¯
1
2
π
i
d
-
1
∫
T
d
-
1
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
,
e
i
t
1
,
…
,
e
i
t
d
-
1
∏
ν
=
1
d
-
1
1
e
i
t
ν
-
a
ν
e
i
2
π
2
-
k
m
ν
+
m
ν
′
d
e
i
t
1
⋯
d
e
i
t
d
-
1
⋮
=
∑
m
′
∈
L
k
P
k
,
m
′
¯
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
+
m
d
-
1
′
.Therefore, the column vector (
〈
f
,
T
k
m
ϕ
k
〉
)
m
∈
L
k can be rewritten as(36)
f
,
T
k
m
ϕ
k
m
∈
L
k
=
A
2
k
P
¯
2
k
.Now, the proof of (30) is concluded by (9), (34), and (36).
For anym
=
(
m
0
,
m
1
,
…
,
m
d
-
1
) and m
′
=
(
m
0
′
,
m
1
′
,
…
,
m
d
-
1
′
)
∈
L
k, it is straightforward to check that (37)
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
+
m
d
-
1
′
=
f
a
0
e
i
2
π
2
-
k
m
0
+
m
0
′
m
o
d
2
k
,
…
,
a
d
-
1
e
i
2
π
2
-
k
m
d
-
1
+
m
d
-
1
′
m
o
d
2
k
,and then A
2
k is a d-level Hankel matrix satisfying (18). Similarly, B
2
k is also a d-level Hankel matrix with property (18).The following lemma is crucial for investigating the matrix expression ofa
k
,
x.Lemma 5.
For anyk
∈
N
0, define a matrix of functions(38)
K
2
k
,
a
≔
1
1
-
a
¯
e
i
·
-
2
π
2
-
k
ξ
+
η
ξ
,
η
∈
N
2
k
,where a
∈
D. Then, for any x
∈
R and m
∈
N
2
k, the value of K
2
k
,
a at x
+
2
π
2
-
k
m satisfies(39)
K
2
k
,
a
x
+
2
π
2
-
k
m
=
K
2
k
,
a
x
S
2
k
m
,where (40)
S
2
k
≔
O
2
k
-
1,1
I
2
k
-
1
1
O
1
,
2
k
-
1
.Proof.
We first prove by induction onm that(41)
S
2
k
m
=
O
2
k
-
m
,
m
I
2
k
-
m
I
m
O
m
,
2
k
-
m
,
m
∈
N
2
k
.Suppose (41) holds with m being replaced by m
-
1. Then, (42)
S
2
k
m
=
S
2
k
S
2
k
m
-
1
=
O
2
k
-
1,1
I
2
k
-
1
1
O
1
,
2
k
-
1
O
2
k
-
m
+
1
,
m
-
1
I
2
k
-
m
+
1
I
m
-
1
O
m
-
1
,
2
k
-
m
+
1
=
O
2
k
-
1,1
I
2
k
-
1
1
O
1
,
2
k
-
1
O
1
,
m
-
1
1
O
1
,
2
k
-
m
O
2
k
-
m
,
m
-
1
O
2
k
-
m
,
1
I
2
k
-
m
I
m
-
1
O
m
-
1,1
O
m
-
1
,
2
k
-
m
=
O
2
k
-
m
,
m
-
1
O
2
k
-
m
,
1
I
2
k
-
m
I
m
-
1
O
m
-
1,1
O
m
-
1
,
2
k
-
m
O
1
,
m
-
1
1
O
1
,
2
k
-
m
=
O
2
k
-
m
,
m
I
2
k
-
m
I
m
O
m
,
2
k
-
m
.Therefore, (41) holds for any m
∈
N
2
k.
For anyξ
,
η
∈
N
2
k, we derive from (41) that (43)
K
2
k
,
a
x
S
2
k
m
ξ
,
η
=
K
2
k
,
a
x
ξ
,
2
k
-
m
+
η
,
0
≤
η
≤
m
-
1
;
K
2
k
,
a
x
ξ
,
η
-
m
,
m
-
1
<
η
≤
2
k
-
1
;
=
1
1
-
a
¯
e
i
x
+
2
π
2
-
k
m
-
2
π
2
-
k
ξ
+
η
=
K
2
k
,
a
x
+
2
π
2
-
k
m
ξ
,
η
.The proof of (39) is concluded.Using (39), we will give the matrix expression of a
k
,
x defined in (29).Theorem 6.
Let{
A
2
k
,
B
2
k
,
P
2
k
} and S
2
k be as in Lemmas 4 and 5, respectively. Define a d-level matrix G
2
k
≔
(
g
m
′
,
m
)
m
′
,
m
∈
L
k by(44)
g
m
′
,
m
=
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
P
2
k
m
′
,where m
=
(
m
0
,
…
,
m
d
-
1
). Then, a
k
,
x in (29) can be expressed by(45)
a
k
,
x
=
G
2
k
B
2
k
t
0
,
x
A
2
k
P
¯
2
k
.Proof.
Using (41) to directly compute (44) gives us that (46)
g
m
′
,
m
=
P
2
k
m
0
+
m
0
′
m
o
d
2
k
,
…
,
m
d
-
1
+
m
d
-
1
′
m
o
d
2
k
,from which we arrive at (47)
g
m
′
,
m
=
g
m
0
+
m
0
′
m
o
d
2
k
,
…
,
m
d
-
1
+
m
d
-
1
′
m
o
d
2
k
,where m
′
=
(
m
0
′
,
…
,
m
d
-
1
′
). Hence, G
2
k is a d-level Hankel matrix with additional property (18).
It follows from (10), (32), (38), and (39) that the value of B
2
k at t
m
,
x defined in (28) can be computed by(48)
B
2
k
t
m
,
x
=
K
2
k
,
a
0
x
0
+
2
π
2
-
k
m
0
⊗
⋯
⊗
K
2
k
,
a
d
-
1
x
0
+
2
π
2
-
k
m
d
-
1
=
K
2
k
,
a
0
x
0
S
2
k
m
0
⊗
⋯
⊗
K
2
k
,
a
d
-
1
x
0
S
2
k
m
d
-
1
=
K
2
k
,
a
0
x
0
⊗
⋯
⊗
K
2
k
,
a
d
-
1
x
0
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
=
B
2
k
t
0
,
x
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
.It is deduced from (30) and (48) that(49)
B
k
f
t
m
,
x
=
B
2
k
t
m
,
x
P
2
k
T
A
2
k
P
¯
2
k
=
B
2
k
t
0
,
x
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
P
2
k
T
A
2
k
P
¯
2
k
=
S
2
k
m
0
⊗
⋯
⊗
S
2
k
m
d
-
1
P
2
k
T
B
2
k
T
t
0
,
x
A
2
k
P
¯
2
k
.By (49), the column vector a
k
,
x defined in (29) can be expressed by (50)
a
k
,
x
=
G
2
k
T
B
2
k
T
t
0
,
x
A
2
k
P
¯
2
k
.Since both B
2
k
(
t
0
,
x
) and G
2
k are d-level Hankel matrices, they are both symmetric. Hence, the proof of (45) is concluded.Note 2.
SinceG
2
k, B
2
k
(
t
0
,
x
), and A
2
k are all d-level Hankel matrices with property (18), by (27) and (45), a
k
,
x can be factorized as(51)
a
k
,
x
=
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
g
~
2
k
Φ
2
k
D
2
k
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
b
~
2
k
Φ
2
k
D
2
k
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
a
~
2
k
Φ
2
k
D
2
k
P
¯
2
k
,where g
~
2
k, b
~
2
k, and a
~
2
k are the transposes of the first rows of G
~
2
k
≔
G
2
k
D
2
k, B
~
2
k
(
t
0
,
x
)
≔
B
2
k
(
t
0
,
x
)
D
2
k, and A
~
2
k
≔
A
2
k
D
2
k, respectively.Based on Note2, a fast algorithm for a
k
,
x will be established as follows.Algorithm 7.
Letϕ
0, {
P
k
,
m
}
m
∈
L
k, and {
P
^
k
,
m
}
m
∈
L
k be as in Lemma 1. According to (51), a
k
,
x defined in (29) can be computed by the following six steps.(1)
Byϕ
0
^
(
m
)
=
∏
ν
=
0
d
-
1
‍
a
¯
ν
m
ν, (11), and (12), compute {
P
^
k
,
m
}
m
∈
L
k, where m
=
(
m
0
,
…
,
m
d
-
1
).
(2)
By implementing IFFT on{
P
^
k
,
m
}
m
∈
L
k, we compute P
2
k.
(3)
SinceΦ
2
k
=
F
2
k
⊗
⋯
⊗
F
2
k and F
2
k is the Fourier matrix of order 2
k, Φ
2
k
∗
g
~
2
k, Φ
2
k
∗
a
~
2
k, and Φ
2
k
∗
b
~
2
k are computed by IFFTs.
(4)
Consider(52)
P
2
k
→
D
2
k
D
2
k
P
¯
2
k
→
F
F
T
Φ
2
k
D
2
k
P
¯
2
k
→
diag
⁡
Φ
2
k
∗
a
~
2
k
diag
⁡
Φ
2
k
∗
a
~
2
k
Φ
2
k
D
2
k
P
¯
2
k
→
I
F
F
T
a
k
,
x
1
≔
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
a
~
2
k
Φ
2
k
D
2
k
P
¯
2
k
.
(5)
Consider(53)
a
k
,
x
1
→
D
2
k
D
2
k
a
k
,
x
1
→
F
F
T
Φ
2
k
D
2
k
a
k
,
x
1
→
diag
⁡
Φ
2
k
∗
b
~
2
k
diag
⁡
Φ
2
k
∗
b
~
2
k
Φ
2
k
D
2
k
a
k
,
x
1
→
I
F
F
T
a
k
,
x
2
≔
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
b
~
2
k
Φ
2
k
D
2
k
a
k
,
x
1
.
(6)
Consider(54)
a
k
,
x
2
→
D
2
k
D
2
k
a
k
,
x
2
→
F
F
T
Φ
2
k
D
2
k
a
k
,
x
2
→
diag
⁡
Φ
2
k
∗
b
~
2
k
diag
⁡
Φ
2
k
∗
g
~
2
k
Φ
2
k
D
2
k
a
k
,
x
2
→
I
F
F
T
a
k
,
x
=
1
Π
2
k
Φ
2
k
∗
diag
⁡
Φ
2
k
∗
g
~
2
k
Φ
2
k
D
2
k
a
k
,
x
2
.
Computational Complexity. It is easy to check that directly computing a
k
,
x through (13) costs O
(
2
3
d
k
) operations. In Algorithm 7, however, FFT and IFFT are used for three and seven times, respectively, which cost O
(
2
k
d
·
d
k
) operations. Meanwhile, the complexity of other operations such as the multiplication of D
2
k and P
¯
2
k is O
(
2
k
d
). Therefore, the computational complexity of Algorithm 7 is O
(
2
k
d
·
d
k
).Note 3.
So far, on computinga
k
,
x, we have not found any numerical algorithm better than Algorithm 7. On the other hand, we notice that Algorithm 7 holds for equally spaced points on T
d. In another occasion, we will establish the fast algorithm for the case of nonequally spaced points.


## 3. Experiment

To check the efficiency of Algorithm7, a numerical computation experiment on the function(55)
f
e
i
t
=
e
i
t
1
1
+
0.3
e
i
t
0
1.5
+
0.4
e
i
t
0
+
2
t
1
,
t
=
t
0
,
t
1
∈
T
2
,will be carried out. Table 1 shows the approximation error ratios (56)
r
a
t
i
o
1
≔
b
k
,
0
-
a
k
,
0
l
2
b
k
,
0
l
2
,
r
a
t
i
o
2
≔
R
e
b
k
,
0
-
a
k
,
0
l
2
R
e
b
k
,
0
l
2
,
r
a
t
i
o
3
≔
Im
⁡
b
k
,
0
-
a
k
,
0
l
2
Im
⁡
b
k
,
0
l
2
,and the times T
1 and T
2 corresponding to different choices of k, a
0, and a
1, where b
k
,
0
=
(
f
(
e
i
t
m
,
0
)
)
m
∈
L
k and T
1 and T
2 are the compute running times cost by directly computing a
k
,
0 and Algorithm 7, respectively. The data on the time cost confirms that Algorithm 7 is faster than direct computing. Therefore, the result in this experiment matches the computational complexity analysis in Section 2. On the other hand, as for the approximation accuracy, as the scale level k increases, the approximation efficiency becomes better. This coincides with the approximation error estimation given in (14).Table 1
The error ratios (i.e., ratio 1, ratio 2, and ratio 3) corresponding to differentk, a
0, and a
1. T
1 and T
2 are the time cost by direct computing and Algorithm 7, respectively.
k
a
0
a
1
ratio 1
ratio 2
ratio 3
T
1
T
2
2
0.01
0.32
0.0704
0.0733
0.0676
0.107 s
0.086 s
3
0.08
0.16
0.0051
0.0051
0.0051
0.214 s
0.103 s
4
0.44
0.41
2.5618
×
1
0
-
5
2.5652
×
1
0
-
5
2.5584
×
1
0
-
5
4.052 s
0.116 s



---

*Source: 5780718-2016-03-15.xml*