Code covered by the BSD License  

Highlights from
Matrix Products Expressed in Terms of Individual Operands

from Matrix Products Expressed in Terms of Individual Operands by Matt J
A class representing products of matrices, internally storing/manipulating them separately.

ProdCascade(opsequence)
function C=ProdCascade(opsequence)
%ProdCascade class constructor.
%
%by Matt Jacobson, 
%Copyright, The University of Michigan, 2005
%
%USAGE:
%
%  P=ProdCascade({A1,A2,...,An})
%
%This creates and 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 evaluate a product 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;
%     toc;
%     %Elapsed time is 0.024294 seconds.
% 
%     tic
%       y2=P*x;
%     toc
%     %Elapsed time is 0.000664 seconds.
% 
%     >>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 low-pass filtering operation as a multiplication with 
%ProdCascade object, P.
%
%       N=2^14;
%       x=rand(N,1);
%       LowPass=ones(N,1);  LowPass(10:end-9)=0;
% 
%       Q=MatrixObj; 
%       Q.Ops.mtimes=@(Q,x) fft(x);
%       Q.Trans.mtimes=@(Q,x) ifft(x);
% 
% 
%       L=spdiags(LowPass,0,speye(N));
%       P=ProdCascade({Q',L,Q});
%        %Operation P*x should be equivalent to ifft(fft(x).*LowPass)
% 
%       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. 

try
 superiorto('MatrixObj');
end

C.opsequence={};

if nargin==0

 C = class(C,'ProdCascade');

elseif isa(opsequence,'ProdCascade')

 C=opsequence;

else

 if ~isrow(opsequence), 
  error('input must be a row array');
 end

%  nonscalars=opsequence(~cellfun(@isscalar, opsequence));
%  for ii=2:length(nonscalars)
%      if ~isscalar(nonscalars{ii}) && size(nonscalars{ii},1)~=size(nonscalars{ii-1},2)
%          error(['Cascade members are not all of dimensions compatible for multiplication'])
%      end
%  end
 
 
 C.opsequence=opsequence;
 C = class(C,'ProdCascade');

end

Contact us