Wavelet Toolbox

Fast Multiplication of Large Matrices

This section illustrates matrix-vector multiplication in the wavelet domain.

• The problem is

• let m be a dense matrix of large size (n, n). We want to perform a large number, L, of multiplications of m by vectors v.

• The idea is

• Stage 1: (executed once) Compute the matrix approximation, sm, at a suitable level k. The matrix will be assimilated with an image.

  Stage 2: (executed L times) divided in the following three steps:

  Compute vector approximation.

  Compute multiplication in wavelet domain.

  Reconstruct vector approximation.

It is clear that when sm is a sufficiently good approximation of m, the error with respect to ordinary multiplication can be small. This is the case in the first example below where m is a magic square. Conversely, when the wavelet representation of the matrix m is dense the error will be large (for example, if all the coefficients have the same order of magnitude). This is the case in the second example below where m is two-dimensional Gaussian white noise. THe figure in Example 1 compares for n = 512, the number of flops required by wavelet based method and by ordinary method versus L.

The function flops is used in the following examples to count floating point operations. This function was available in earlier MATLAB versions and is now obsolete. So, the instructions that call this function will not give the displayed results anymore. But, these examples still illustrate the advantage of the wavelet method in matrix multiplication.

Example 1: Effective Fast Matrix Multiplication

• n = 512; 
  lev = 5; 
  wav = 'db1';
  
  % Wavelet based matrix multiplication by a vector: 
  % a "good" example 
  % Matrix is magic(512) Vector is (1:512)
  
  m = magic(n); 
  v = (1:n)'; 
  [Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wav);
  
  % ordinary matrix multiplication by a vector. 
  flops(0);
  p = m * v; 
  flomv = flops
  
  flomv = 
         524288
  
  % Compute matrix approximation at level 5. 
  flops(0) 
  sm = m;
  for i = 1:lev 
      sm = dyaddown(conv2(sm,Lo_D),'c'); 
      sm = dyaddown(conv2(sm,Lo_D'),'r'); 
  end 
  flmapp = flops
  
  flmapp = 
          2095154
  
  % The three steps: 
  % 1. Compute vector approximation. 
  % 2. Compute multiplication in wavelet domain. 
  % 3. Reconstruct vector approximation.
  
  
  flops(0) 
  sv = v; 
  for i = 1:lev, sv = dyaddown(conv(sv,Lo_D)); end 
  sp = sm * sv; 
  for i = 1:lev, sp = conv(dyadup(sp),Lo_R); end 
  sp = wkeep(sp,length(v)); 
  flwmv = flops
  
  flwmv = 
          9058
  
  % Plot ordinary versus wavelet based m*v flops in loglog.
  
  
  

Figure 3-1: Wavelet Based Matrix Multiplication by a Vector

• % Relative square norm error in percent when using wavelets. 
  rnrm = 100 * (norm(p-sp)/norm(p))
  
  rnrm = 
          2.9744e-06
  

Example 2: Ineffective Fast Matrix Multiplication

The commands used are the same as in Example 1, but applied to a new matrix m.

• % Wavelet based matrix multiplication by a vector: 
  % a "bad" example with a randn matrix.
  % Change the matrix m of example1 using:
  randn('state',0);
  m = randn(n,n);
  

Then, you obtain

• % Relative square norm error in percent 
  rnrm = 100 * (norm(p-sp)/norm(p))
  
  rnrm = 
          99.2137
  

Compressing Images

Wavelets in Action: Examples and Case Studies

© 1994-2005 The MathWorks, Inc.