Matrix Products Expressed in Terms of Individual Operands

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

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.