Skip to content

Commit 1cba88a

Browse files
jagotjagot
authored
Transforms between bases (#63)
* Added orthogonality trait * Implemented transforms between bases * Documentation for basis transforms * Added tests of basis transforms * Improved docs * Fixed tests on Julia 1.6 * Bump version * Test skipped line
1 parent a68a60e commit 1cba88a

15 files changed

+629
-2
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "CompactBases"
22
uuid = "2c0377a8-7469-4ebd-be0f-82e501f20078"
33
authors = ["Stefanos Carlström <stefanos.carlstrom@gmail.com>"]
4-
version = "0.3.12"
4+
version = "0.3.13"
55

66
[deps]
77
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"

‎docs/make.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,7 @@ makedocs(;
3535
]
3636
],
3737
],
38+
"Basis transforms" => "basis_transforms.md"
3839
],
3940
repo="https://github.com/JuliaApproximation/CompactBases.jl/blob/{commit}{path}#L{line}",
4041
sitename="CompactBases.jl",

‎docs/plots.jl

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -269,6 +269,75 @@ function diagonal_operators()
269269
savedocfig("diagonal_operators")
270270
end
271271

272+
function basis_transforms()
273+
function compare_bases(x, bases, functions)
274+
m = length(bases)
275+
n = length(functions)
276+
277+
cs = [B \ f.(axes(B,1))
278+
for B in bases, f in functions]
279+
es = extrema.(cs)
280+
ymins = [minimum(first.(es[i,:])) for i = 1:m]
281+
ymaxs = [maximum(last.(es[i,:])) for i = 1:m]
282+
283+
cfigure("basis comparison", figsize=(10,16)) do
284+
for i = 1:m
285+
B = bases[i]
286+
a = axes(B,1)
287+
χ = B[x,:]
288+
xl = [mean_position(x, view(χ, :, j)) for j in 1:size(B,2)]
289+
csubplot(m,n+1,(i,1), nox=i<m) do
290+
plot(x, χ)
291+
i == m && xlabel(L"x")
292+
reltext(0.1, 0.8, string('A'+i-1))
293+
end
294+
for j = 1:n
295+
f = functions[j]
296+
c = cs[i,j]
297+
csubplot(2m,n+1,(2i-1,j+1), nox=true, noy=j<n) do
298+
plot(x, χ*c)
299+
plot(x, f.(x), "k--")
300+
margins(0,0.1)
301+
ylim(-0.2,1.25)
302+
axes_labels_opposite(:y)
303+
i == m && xlabel(L"x")
304+
end
305+
csubplot(2m,n+1,(2i,j+1), nox=i<m, noy=j<n) do
306+
plot(xl, c, ".-")
307+
for i′ = 1:m
308+
i′ == i && continue
309+
A = bases[i′]
310+
d = basis_transform(B, A, cs[i′,j])
311+
plot(xl, d, ".--", label=string('A'+i′-1))
312+
end
313+
xlim(x[1],x[end])
314+
dy = (ymaxs[i]-ymins[i])
315+
ylim(ymins[i]-0.1dy, ymaxs[i]+0.1dy)
316+
legend(loc=2)
317+
axes_labels_opposite(:y)
318+
i == m && xlabel(L"x")
319+
end
320+
end
321+
end
322+
end
323+
end
324+
325+
A = BSpline(LinearKnotSet(7, 0, 2, 11))
326+
B = FEDVR(range(0.0, 2.5, 5), 6)[:,2:end-1]
327+
C = StaggeredFiniteDifferences(0.03, 0.3, 0.01, 3.0)
328+
D = FiniteDifferences(range(0.25,1.75,length=31))
329+
E = BSpline(LinearKnotSet(9, -1, 3, 15))[:,2:end]
330+
331+
gauss = x -> exp(-(x-1.0).^2/(2*0.2^2))
332+
u = x -> (0.5 < x < 1.5)*one(x)
333+
334+
x = range(-1.5, stop=3.5, length=1001)
335+
336+
compare_bases(x, (A,B,C,D,E), (gauss,u))
337+
338+
savedocfig("basis_transforms")
339+
end
340+
272341
macro echo(expr)
273342
println(expr)
274343
:(@time $expr)
@@ -285,3 +354,4 @@ include("bspline_plots.jl")
285354
include("fd_plots.jl")
286355
@echo densities()
287356
@echo diagonal_operators()
357+
@echo basis_transforms()

‎docs/src/basis_transforms.md

Lines changed: 221 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,221 @@
1+
# Transforms between different bases
2+
3+
Sometimes, we have a function ``f`` expanded over the basis functions
4+
of one basis ``A``,
5+
```math
6+
f_A = A\vec{c}
7+
```
8+
and wish to re-express as an expansion over the
9+
basis functions of another basis ``B``,
10+
```math
11+
f_B = B\vec{d}.
12+
```
13+
We thus want to solve ``f_A = f_B`` for ``\vec{d}``; this is very
14+
easy, and follows the same reasoning as in [Solving
15+
equations](@ref). We begin by projecting into the space ``B`` by using
16+
the projector ``BB^H``:
17+
```math
18+
BB^HA\vec{c} =
19+
BB^HB\vec{d}.
20+
```
21+
This has to hold for any ``B``, as long as the basis functions of
22+
``B`` are linearly independent, and is then equivalent to
23+
```math
24+
B^HA\vec{c} =
25+
B^HB\vec{d},
26+
```
27+
the solution of which is
28+
```math
29+
\vec{d} =
30+
(B^HB)^{-1}
31+
B^HA
32+
\vec{c}
33+
\defd
34+
S_{BB}^{-1}
35+
S_{BA}
36+
\vec{c}.
37+
```
38+
39+
[`basis_overlap`](@ref) computes ``S_{BA}=B^HA`` and
40+
[`basis_transform`](@ref) the full transformation matrix
41+
``(B^HB)^{-1}B^HA``; sometimes it may be desirable to perform the
42+
transformation manually in two steps, since the combined matrix may be
43+
dense whereas the constituents may be relatively sparse. For
44+
orthogonal bases, ``B^HB`` is naturally diagonal.
45+
46+
There is no requirement that the two bases span the same intervals,
47+
i.e. `axes(A,1)` and `axes(B,1)` need not be fully overlapping. Of
48+
course, if they are fully disjoint, the transform matrix will be
49+
zero. Additionally, some functions may not be well represented in a
50+
particular basis, e.g. the step function is not well approximated by
51+
B-splines (of order ``>1``), but poses no problem for
52+
finite-differences, except at the discontinuity. To judge if a
53+
transform is faithful thus amount to comparing the transformed
54+
expansion coefficients with those obtained by expanding the function
55+
in the target basis directly, instead of comparing reconstruction
56+
errors.
57+
58+
## Examples
59+
60+
We will now illustrate transformation between B-splines, FEDVR, and
61+
various finite-differences:
62+
```julia
63+
julia> A = BSpline(LinearKnotSet(7, 0, 2, 11))
64+
BSpline{Float64} basis with typename(LinearKnotSet)(Float64) of order k = 7 on 0.0..2.0 (11 intervals)
65+
66+
julia> B = FEDVR(range(0.0, 2.5, 5), 6)[:,2:end-1]
67+
FEDVR{Float64} basis with 4 elements on 0.0..2.5, restricted to elements 1:4, basis functions 2..20 ⊂ 1..21
68+
69+
julia> C = StaggeredFiniteDifferences(0.03, 0.3, 0.01, 3.0)
70+
Staggered finite differences basis {Float64} on 0.0..3.066973760326056 with 90 points with spacing varying from 0.030040496962651875 to 0.037956283408020486
71+
72+
julia> D = FiniteDifferences(range(0.25,1.75,length=31))
73+
Finite differences basis {Float64} on 0.2..1.8 with 31 points spaced by Δx = 0.05
74+
75+
julia> E = BSpline(LinearKnotSet(9, -1, 3, 15))[:,2:end]
76+
BSpline{Float64} basis with typename(LinearKnotSet)(Float64) of order k = 9 on -1.0..3.0 (15 intervals), restricted to basis functions 2..23 ⊂ 1..23
77+
```
78+
79+
We use these bases to expand two test functions:
80+
```math
81+
f_1(x) = \exp\left[
82+
-\frac{(x-1)^2}{2\times0.2^2}
83+
\right]
84+
```
85+
and
86+
```math
87+
f_2(x) = \theta(x-0.5) - \theta(x-1.5),
88+
```
89+
where ``\theta(x)`` is the [Heaviside step
90+
function](https://en.wikipedia.org/wiki/Heaviside_step_function).
91+
92+
```julia
93+
julia> gauss(x) = exp(-(x-1.0).^2/(2*0.2^2))
94+
gauss (generic function with 1 method)
95+
96+
julia> u(x) = (0.5 < x < 1.5)*one(x)
97+
u (generic function with 1 method)
98+
```
99+
100+
If for example, we first expand ``f_1`` in B-splines (``A``), and then wish to
101+
project onto the staggered finite-differences (``C``), we do the
102+
following:
103+
104+
```julia
105+
julia> c = A \ gauss.(axes(A, 1))
106+
17-element Vector{Float64}:
107+
-0.0003757237742430358
108+
0.001655433792947186
109+
-0.003950561899347059
110+
0.00714774775995981
111+
-0.011051108016208545
112+
0.0173792655408227
113+
0.04090160980771093
114+
0.6336420285181446
115+
1.3884690430556483
116+
0.6336420285181449
117+
0.04090160980771065
118+
0.01737926554082302
119+
-0.011051108016208311
120+
0.007147747759958865
121+
-0.00395056189934631
122+
0.0016554337929467365
123+
-0.00037572377424269145
124+
125+
julia> S_CA = basis_transform(C, A)
126+
90×17 SparseArrays.SparseMatrixCSC{Float64, Int64} with 1530 stored entries:
127+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
128+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
129+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
130+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
131+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
132+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
133+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
134+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
135+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
136+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
137+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
138+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
139+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
140+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
141+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
142+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
143+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
144+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
145+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
146+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
147+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
148+
⣿⣿⣿⣿⣿⣿⣿⣿⡇
149+
⠛⠛⠛⠛⠛⠛⠛⠛⠃
150+
151+
julia> S_CA*c
152+
90-element Vector{Float64}:
153+
3.805531055898412e-5
154+
1.6890846360940058e-5
155+
-4.5403814886370545e-5
156+
-4.121425310424971e-5
157+
1.2245471085512564e-5
158+
6.482348697494439e-5
159+
8.896755347867799e-5
160+
9.703404969825733e-5
161+
0.0001290654548866452
162+
0.00023543152355538763
163+
0.0004663267687150117
164+
⋮
165+
0.0
166+
0.0
167+
0.0
168+
0.0
169+
0.0
170+
0.0
171+
0.0
172+
0.0
173+
0.0
174+
0.0
175+
176+
julia> basis_transform(C, A, c) # Or transform without forming matrix
177+
90-element Vector{Float64}:
178+
3.80553105589841e-5
179+
1.6890846360940065e-5
180+
-4.540381488637052e-5
181+
-4.121425310424972e-5
182+
1.2245471085512601e-5
183+
6.482348697494439e-5
184+
8.896755347867799e-5
185+
9.703404969825729e-5
186+
0.00012906545488664526
187+
0.0002354315235553877
188+
0.00046632676871501176
189+
⋮
190+
0.0
191+
0.0
192+
0.0
193+
0.0
194+
0.0
195+
0.0
196+
0.0
197+
0.0
198+
0.0
199+
0.0
200+
```
201+
202+
In the plots below, we show each basis (labelled ``A-E``), and their
203+
respective reconstructions of ``f_1`` and ``f_2`` as solid red lines,
204+
with the true functions shown as dashed black lines. For each
205+
basis–function combination, shown are also the expansion coefficients
206+
as red points joined by solid lines, as well as the expansion
207+
coefficients transformed from the other bases, as labelled in the
208+
legends. As an example, we see that first expanding in
209+
finite-differences ``D`` and then transforming to B-splines ``A`` is
210+
numerically unstable, leading to very large coefficients at the edges
211+
of ``D``. Going _to_ finite-differences is always stable, even though
212+
the results may not be accurate.
213+
214+
![Transforms between bases](../figures/basis_transforms.svg)
215+
216+
## Reference
217+
218+
```@docs
219+
basis_overlap
220+
basis_transform
221+
```

‎docs/src/overview.md

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,9 @@ julia> S = B'B
7474
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.1 ⋅ ⋅
7575
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.1 ⋅
7676
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 0.1
77+
78+
julia> orthogonality(B)
79+
CompactBases.Orthogonal()
7780
```
7881

7982
In contrast, a non-orthogonal basis such as B-splines has a _banded_
@@ -111,6 +114,9 @@ julia> S = B'B
111114
⋅ ⋅ 9.79816e-7 0.00036737 0.0100953 0.0349662 0.0464947 0.0246054 0.00347002
112115
⋅ ⋅ ⋅ 1.65344e-6 0.000811839 0.00747795 0.0246054 0.0331746 0.0139286
113116
⋅ ⋅ ⋅ ⋅ 1.32275e-5 0.000365961 0.00347002 0.0139286 0.0222222
117+
118+
julia> orthogonality(B)
119+
CompactBases.NonOrthogonal()
114120
```
115121

116122
where the bandwidth depends on the B-spline order.

‎src/CompactBases.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,7 @@ using RecipesBase
5050
vandermonde(B::Basis) = B[locs(B),:]
5151

5252
include("distributions.jl")
53+
include("orthogonality.jl")
5354
include("restricted_bases.jl")
5455
include("types.jl")
5556
include("unit_vectors.jl")
@@ -78,6 +79,8 @@ include("inner_products.jl")
7879
include("densities.jl")
7980
include("linear_operators.jl")
8081

82+
include("basis_transforms.jl")
83+
8184

8285
"""
8386
centers(B)

0 commit comments

Comments
 (0)