The algorithm is described in many places, see e.g., [3,2,13]. The basic idea is to represent both the matrix and the vector in a wavelet basis and multiply the transformed matrix and transformed vector. This algorithm requires access to the intermediate results of the low-pass filter H (which are normally discarded in a wavelet transform). The function fwt1ns returns a vector which contains also these blocks, hence the name non-standard.
The matrix is transformed with fwt2 as usual, the vector with fwt1ns, the multiplication is done by nsmult, and the product is inverse transformed with ifwt1ns. See also the demo wavmultd (and nsexampl for a source of test matrices).
The output vector of fwt1ns is partitioned as follows: Let H be the low-pass filter corresponding to the scaling function followed by down-sampling, G the high-pass filter corresponding to the wavelet (followed by down-sampling). Then (h and g contain the filter coefficients)
w = fwt1ns(f,h,g)
yields
Example: The product of
Aij=1/|i-j|,
,
Aii=0, with
a random vector.
[h,g] = wavecoef('coi',30);
[i,j] = meshgrid(0:31);
A = 1./(i-j); A(1:33:32*32) = zeros(1,32);
WA = fwt2(A,h,g);
f = rand(32,1);
wf = fwt1ns(f,h,g);
p = ifwt1ns(nsmult(WA,wf),h,g);
norm(p-A*f)
The error was
9.81109e-15
.
Note that fwt1ns followed by ifwt1ns is not the identity. There must be an intermediate multiplication by the transform of the identity matrix (see help nsexampl).
In the demo wavmultd the user can choose which operator to study and specify at what level the wavelet coefficients are truncated. Figure 10 shows the transformed matrix on the left. Only those matrix elements that are larger than the specified (relative) truncation level are displayed and used when computing the product. The graphs on the right are, from top to bottom, the original vector (chosen randomly), the product when there is no truncation, product computed with the truncated matrix, and the error between the last two.