Wavelet Compression of Integral Operators
This numerical tour explores approximation of full matrices by sparse matrices over wavelets coefficients.
Contents
Installing toolboxes and setting up the path.
You need to download the following files: signal toolbox and general toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.
For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.
Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.
Execute this line only if you are using Matlab.
getd = @(p)path(p,path); % scilab users must *not* execute this
Then you can add the toolboxes to the path.
getd('toolbox_signal/'); getd('toolbox_general/');
Integral Operators
A 1D integral operator is discretized as a (n,n) 2D matrix that has many non zero entries (full matrix). It is thus very different from differential operators that are very sparse, with only O(n) non zero entries.
A typical integral operator is the Hilbert transform, that corresponds to the convolution with the kernel 1/(x-y), here with zero boundary conditions outside the domain.
% size of the domain n = 1024; % kernel [Y,X] = meshgrid(1:n,1:n); K = 1 ./ max(abs(X-Y),1);
This kernel can be displayed as a smooth image, that decays away from the diagonal.
clf; subplot(2,1,1); imageplot(K); subplot(2,1,2); plot(K(n/2,:)); axis('tight'); title('Central row');
We can compute a signal and the action of the kernel on the signal.
x = load_signal('piece-regular', n);
y = K*x;
We display both.
clf; subplot(2,1,1); plot(x); axis('tight'); title('x'); subplot(2,1,2); plot(y); axis('tight'); title('y=K*x');
From the knowledge of y one can solve the system A*x=y to retrive x. This can be done with a direct method, for instance using the operator \, that uses a LU decomposition.
x1 = K \ y; disp(['Error (should be 0) = ' num2str(norm(x-x1)/norm(x)) '.']);
Error (should be 0) = 7.3282e-13.
Standard Wavelet Representation of Operators
A wavelet representation of an operator corresponds to a change of basis, by replacing the Dirac pointwise basis by a wavelet basis.
The wavelet representation works for orthogonal wavelets (otherwise one needs to use both primal and dual wavelets, which complicates the computations). We select the Haar transform.
options.filter = 'haar';
If one denotes by W the matrix of the pointwise wavelet transform and W' its transpose (which is also its inverse because the transform is orthogonal), the matrix in the transformed domain is represented formally as KW = W'*K*W.
Computing W'*K*W corresponds to applying the 1D wavelet transform to all the columns of K and then to all the rows. This is called the separable 2D wavelet transform. This should not be confounded with the isotropic 2D wavelet transform, which does not corresponds to a change of basis.
KW can be computed by applying the fast wavelet transform, which results in an O(n^2) algorithm. The result is called the standard wavelet representation of the operator, which should not be confounded with the non standard representation.
Jmin = 1; options.separable = 1; KW = perform_wavelet_transf(K, Jmin, +1, options);
One can display the resulting compressed operator.
clf; plot_wavelet(KW, Jmin, options);
Wavelet Domain Approximation of Operators
Since the kernel K is smooth, most of its wavelet coefficients are small. It can thus be well approximated by only a few coefficients.
Exercice 1: (the solution is exo1.m) Compute a best M=.05*n^2 terms approximation KWT of KW by thresholding KW using a well chosen threshold T.
exo1;
The approximated kernel KT is obtained by applying the inverse wavelet transform.
KT = perform_wavelet_transf(KWT, Jmin, -1, options);
Compare original and approximation.
clf; imageplot(K, 'Original', 1,2,1); imageplot(K - KT, 'Original - Approximation', 1,2,2);
Important: the computation of the approximated kernel KT is done here just for display purpose. In practical situations, one only works over the transformed domain (in particular one solve equations y=K*x over the transformed domain), so KT is not needed explicitely.
Approximation Accuracy Study
Here we study how the quality of the solution xT of the system KT*xT=y decreases when the number of coefficients M of the approximated matrix K decreases.
The system is solved over the wavelet domain.
% transform the observations yW = perform_wavelet_transf(y, Jmin, +1, options); % solve the system xTW = KWT \ yW; % transform back the solution xT = perform_wavelet_transf(xTW, Jmin, -1, options);
Display.
clf; subplot(2,1,1); plot(x); axis('tight'); title('x'); subplot(2,1,2); plot(xT); axis('tight'); title(['Approximation, M/n^2=' num2str(M/n^2, 2) ]);
Exercice 2: (the solution is exo2.m) Study the approximation error decay with the increase of the number of coefficients.
exo2;
Linear System Inversion Speed-up
Decreasing the value of the number of coefficients M used to approximate the operator deteriorates the quality of the solution, but it sparsifies the matrix (many entries are 0). Using a sparse linear solver takes advantage of this sparsity and reduces the computation time.
Turn the matrix into a sparse one.
KWT = sparse(KWT);
We can compare the speed taken by spacial inversion of the system y=K*x and wavelet domain inverse yW=KW*xW.
tic; x1 = K \ y; t1 = toc(); tic; x1 = KWT \ yW; t2 = toc(); disp(['Speed up of wavelet inversion : ' num2str(1-t2/t1) '.']);
Speed up of wavelet inversion : -0.29934.
Exercice 3: (the solution is exo3.m) Compute the speed up for an increasing value of the number of coefficients M. You can average over several trials to get an accurate value for the speed up.
exo3;
Spacially Varying Convolution
A spacially varying convolution is a blurring operator that is not translation invariant. The convolution can be accelerated over the spacial domain.
Such a blurring operator is called "Foveation". Its wavelet representation is studied in
Chang, Mallat, Yap, Wavelet Foveation, Applied and Computational Harmonic Analysis, Volume 9, Number 3, October 2000 , pp. 312-335(24)
Load a high resolution signal.
n = 2048;
x = load_signal('piece-polynomial', n);
At each spacial location a width of a Gaussian kernel is defined. We use here a slowly varying kernel, where the resolution is finer in the center.
smin = 1; smax = .03*n; sigma = [linspace(1,0,n/2)';linspace(0,1,n/2)']; sigma = smin + (smax-smin)*sigma.^2;
The kernel matrix is defined by assigning a Gaussian kernel to each row.
K = zeros(n); t = 1:n; for i=1:n K(i,:) = exp( -(t-i).^2 / (2*sigma(i)^2) ); K(i,:) = K(i,:)/sum(K(i,:)); end
Display the kernel.
v = sort(K(:)); v = v(round(linspace(1,n^2,12))); clf; subplot(1,2,1); imageplot(K); axis on; title('K (graylevels)'); subplot(1,2,2); a = contour(K, v); axis ij; axis square; title('K (levelsets)');
Blurring of a signal.
y = K*x;
Display the action of the kernel.
clf; subplot(2,1,1); plot(x); axis('tight'); title('x'); subplot(2,1,2); plot(y); axis('tight'); title('y=K*x');
For later use, store the average application time.
t0 = 0; ntrials = 50; for i=1:ntrials tic; y = K*x; t0 = t0 + toc(); end t0 = t0 / ntrials;
Exercice 4: (the solution is exo4.m) Compute the approximation with M=.05*n^2 coefficients of K over the wavelet domain.
exo4;
Exercice 5: (the solution is exo5.m) Compute the application of the operator yT=KT*xW over the wavelet domain using compressed operator. Compute the speed up, taking into acount the sparse matrix multiplication and the forward/inverse 1D wavelet transforms (not the 2D matrix transform and the thresholding, which are considered as pre-processings).
exo5;
