Fast Numerical Algorithms using Wavelet Transforms

This demonstration is based on the paper Fast wavelet transforms and numerical algorithms by G. Beylkin, R. Coifman, and V. Rokhlin and the lecture notes Wavelets and Fast Numerical ALgorithms by G. Beylkin.

Assume we have a matrix $M \in \mathbb{R}^{n \times n}$, and a vector $x \in \mathbb{R}^n$, there are two ways to compute $y = Mx$:

Use the regular matrix multiplication, which takes $O(n^2)$ time.
Transform both $M$ and $x$ to their standard/nonstandard forms $\tilde M$ and $\tilde x$ respectively. Compute the multiplication of $\tilde y = \tilde M \tilde x$, then find the inverse of $\tilde y$.If the matrix $\tilde M$ is sparse, this can take $O(n)$ time.

Examples and Comparisons

In our following examples, we compare the methods described in the previous section. We focus on two important metrics: time taken and accuracy of the algorithm.

We start off by setting up the matrix M, the vector x, and the wavelet wt. As this algorithm is more efficient when $n$ is large, we set $n=2048$.

using Wavelets, WaveletsExt, Random, LinearAlgebra, BenchmarkTools

# Construct matrix M, vector x, and wavelet wt
function data_setup(n::Integer, random_state = 1234)
    M = zeros(n,n)
    for i in axes(M,1)
        for j in axes(M,2)
            M[i,j] = i==j ? 0 : 1/abs(i-j)
        end
    end

    Random.seed!(random_state)
    x = randn(n)

    return M, x
end

M, x = data_setup(2048)
wt = wavelet(WT.haar)

The following code block computes the regular matrix multiplication.

y₀ = M * x

The following code block computes the nonstandard wavelet multiplication. Note that we will only be running 4 levels of wavelet transform, as running too many levels result in large computation time, without much improvement in accuracy.

L = 4
NM = mat2sparseform_nonstd(M, wt, L)
y₁ = nonstd_wavemult(NM, x, wt, L)

The following code block computes the standard wavelet multiplication. Once again, we will only be running 4 levels of wavelet transform here.

SM = mat2sparseform_std(M, wt, L)
y₂ = std_wavemult(SM, x, wt, L)

Performance Tip

Running more levels of transform (i.e. setting large values of L) does not necessarily result in more accurate approximations. Setting L to be a reasonably large value (e.g. 3 or 4) not only reduces computational time, but could also produce better results.

We assume that the regular matrix multiplication produces the true results (bar some negligible machine errors), and we compare our results from the nonstandard and standard wavelet multiplications based on their relative errors ($\frac{||\hat y - y_0||_2}{||y_0||_2}$).

julia> relativenorm(y₁, y₀)       # Comparing nonstandard wavelet multiplication with true value0.0015112400000964752
julia> relativenorm(y₂, y₀)       # Comparing standard wavelet multiplication with true value0.0010877481648021538

Lastly, we compare the computation time for each algorithm.

julia> @benchmark $M * $x                          # Benchmark regular matrix multiplicationBenchmarkTools.Trial: 2129 samples with 1 evaluation.
 Range (min … max):  1.638 ms …  10.458 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.248 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.320 ms ± 445.544 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

             ▂██▇▇▆▂▁
  ▂▂▂▂▃▂▂▃▄▅█████████▇▆▅▄▃▃▃▃▂▂▃▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▁▂▂▂▁▂▂▂▂▁▂▁▂ ▃
  1.64 ms         Histogram: frequency by time        4.04 ms <

 Memory estimate: 16.12 KiB, allocs estimate: 1.
julia> @benchmark nonstd_wavemult($NM, $x, $wt, L) # Benchmark nonstandard wavelet multiplicationBenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  269.601 μs …  12.535 ms  ┊ GC (min … max): 0.00% … 87.61%
 Time  (median):     290.101 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   312.287 μs ± 338.284 μs  ┊ GC (mean ± σ):  2.69% ±  2.50%

     ▁▇██▇▆▅▄▃▁▁ ▁▁▂▂▂▂▂▁▁▁▁                                    ▂
  ▄▅▅███████████████████████████▇▇▆▇▆▆▆▇▆▆▆▆▆▅▇▇▇▆▆▇▅▆▆▅▅▅▅▅▅▄▄ █
  270 μs        Histogram: log(frequency) by time        451 μs <

 Memory estimate: 129.19 KiB, allocs estimate: 28.
julia> @benchmark std_wavemult($SM, $x, $wt, L)    # Benchmark standard wavelet multiplicationBenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  198.100 μs …  10.477 ms  ┊ GC (min … max): 0.00% … 89.49%
 Time  (median):     206.701 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   220.090 μs ± 210.798 μs  ┊ GC (mean ± σ):  1.66% ±  1.80%

  ▅▇██▇▇▆▅▄▂▁     ▁▁▁▁▁▁▁▁                                      ▂
  ████████████▇▇▇▇████████▇███▇▇▇▇▇▇▇▇▇▆▆▅▆▆▆▃▅▆▆▅▆▆▅▆▅▅▅▅▄▃▃▅▅ █
  198 μs        Histogram: log(frequency) by time        335 μs <

 Memory estimate: 65.14 KiB, allocs estimate: 13.

As we can see above, the nonstandard and standard wavelet multiplication methods are significantly faster than the regular method and provides reasonable approximations to the true values.

Performance Tips

This method should only be used when the dimensions of $M$ and $x$ are large. Additionally, this would be more useful when one aims to compute $Mx_0, Mx_1, \ldots, Mx_k$ where $M$ remains constant throughout the $k$ computations. In this case, it is more efficient to use

NM = mat2sparseform_nonstd(M, wt)
y = nonstd_wavemult(NM, x, wt)

instead of

y = nonstd_wavemult(M, x, wt)