Code covered by the BSD License  

Highlights from
Matrix Products Expressed in Terms of Individual Operands

Be the first to rate this file! 1 Download (last 30 days) File Size: 11.5 KB File ID: #29498

Matrix Products Expressed in Terms of Individual Operands

by

 

23 Nov 2010 (Updated )

A class representing products of matrices, internally storing/manipulating them separately.

| Watch this File

File Information
Description

This submission defines a class for representing products of matrices (or of any objects that know how to multiply) when it is more efficient to store and manipulate the matrices separately.

 USAGE:
 
   P=ProdCascade({A1,A2,...,An})
 
This creates an object P representing the matrix product A1*A2*...*An. Internally, however, the object stores/manipulates the matrices {A1,A2,...,An} individually, which can be more efficient memory-wise than when the product is fully evaluated. Furthermore, when evaluating a matrix-vector product A1*A2*...*An*x, it can be more efficient speed-wise to exploit associativity and compute the result using n successive products.
 
A simple example of this idea is the outer product u*v.' of column vectors u and v. Clearly, it is more efficient to multiply the outer product with a vector x as u*(v.'*x) than as (u*v.')*x. The following gives a demonstration of how the class exploits this.
 
      u=rand(4000,1);
      v=rand(4000,1);
      x=rand(4000,1);
  
      Pmat=u*v.';
      P=ProdCascade({u,v.'});
  
  
      tic
         y1=Pmat*x;
         z1=Pmat.'*y1;
      toc;
      %Elapsed time is 0.050854 seconds.

  

      tic
         y2=P*x;
         z2=P.'*y2;
      toc
      %Elapsed time is 0.000457 seconds.

 
 

       PercentError =100*norm([y2;z2]-[y1;z1])/norm([y1;z1]), %=7.8336e-014

  

      >>whos Pmat P
  
            Name Size Bytes Class Attributes
  
            P 4000x4000 64244 ProdCascade
            Pmat 4000x4000 128000000 double
 
 
Other situations where it can be better to store a matrix decomposed into a product would include, for example, when a large full matrix has a sparse decomposition.
 
Methods overloaded by the class include mtimes (*) ,tranpose (.') , ctranpose (') , det, and inv, which exploit the fact that these operations more or less distribute across matrix products. Also, some methods allow the object P to be manipulated as if it were the cell array {A1,...,An} that generated the object initially. For example, the expression A1=P{1} will extract the first operand of the product.
 
The ProdCascade class can be useful not just for holding products of matrices, but also of any operator objects that know how to multiply. Here is an example that uses my MatrixObj class
 
  http://www.mathworks.com/matlabcentral/fileexchange/26611-on-the-fly-definition-of-custom-matrix-objects
 
to represent a frequency-domain filtering operation as a multiplication with ProdCascade object, P.

       %%data
        N=2^14;
        x=rand(N,1);
        LowPass=ones(N,1); LowPass(10:end-9)=0;
  
       %%Make FFT operator
        Q=MatrixObj;
        Q.Ops.mtimes=@(Q,x) fft(x);
        Q.Trans.mtimes=@(Q,x) ifft(x);
  
        %%Make filtering operator, P*x should be equivalent to ifft(fft(x).*LowPass)

        L=spdiags(LowPass,0,speye(N));
        P=ProdCascade({Q',L,Q});
   
        isequal( P*x,ifft(fft(x).*LowPass) )%=1
  
 
 DISCAIMER: Error checking is never done to see whether the operators in a ProdCascade are compatible for successive multiplication.

MATLAB release MATLAB 7.11 (R2010b)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Updates
01 Dec 2010

Modified the description page. No new code to download.

26 Dec 2010

Edited the copywrite info. No new code.

Contact us