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:

  1. Compute vector approximation.
  2. Compute multiplication in wavelet domain.
  3. 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 floating point operations (flops) required by wavelet based method and by ordinary method versus L.

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. 
  p = m * v; 
  
  % The number of floating point operations used is 524,288
  
  % Compute matrix approximation at level 5. 
  sm = m;
  for i = 1:lev 
      sm = dyaddown(conv2(sm,Lo_D),'c'); 
      sm = dyaddown(conv2(sm,Lo_D'),'r'); 
  end 
  
  % The number of floating point operations used is 2,095,154
  
  % The three steps: 
  % 1. Compute vector approximation. 
  % 2. Compute multiplication in wavelet domain. 
  % 3. Reconstruct vector approximation.
  
  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)); 
  
  % Now, the number of floating point operations used is 9058
  
  % 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
  

Provide feedback about this page

Compressing Images

Wavelets in Action: Examples and Case Studies

Recommended Products

• MATLAB
• Signal Processing Toolbox
• Optimization Toolbox
• Statistics Toolbox
• Control System Toolbox

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit