5.0

5.0 | 1 rating Rate this file 22 downloads (last 30 days) File Size: 55.51 KB File ID: #25969

Efficient Object-Oriented Kronecker Product Manipulation

by Matt J

 

28 Nov 2009 (Updated 23 Feb 2010)

Code covered by the BSD License  

A class for efficient manipulation of N-fold Kronecker products in terms of their operands only.

Download Now | Watch this File

File Information
Description

This is a class for efficiently representing and manipulating N-fold Kronecker products of matrices (or of objects that behave like matrices) in terms of their operands only.
 
Given matrices {A,B,C,D,...} and a scalar s, an object M of this class can be used to represent
 
    Matrix = s * A kron B kron C kron D kron ... (Eq. 1)
 
where "A kron B" denotes kron(A,B), the Kronecker product of A and B. Internally, however, M stores the operands {s,A,B,C,D,...} separately, which is typically far more byte-compact than numerically expanding out the RHS of Eq. 1.
 
Furthermore, many mathematical manipulations of Kronecker products are more efficient when done in terms of {s,A,B,C,D,...} separately than when done with the explicit numerical form of M as a matrix. The class overloads a number of methods and math operators in a way that exploits the Kronecker product structure accordingly.
 
Among these methods/operators are: mtimes (*), times (.*) , tranpose (.') , ctranpose (') , rdivide (./), ldivide (.\), mldivide (\), mrdivide (/), inv, pinv, power, mpower, norm, sum, cond, eig, svd, abs, nnz, orth, chol, lu, qr, full, sparse, ...
 
Some restrictions apply to these overloads. In particular, bi-operand math operations involving two KronProd objects, e.g. M1*M2, typically require the operands of each KronProd to be of compatible sizes. However, I find these restrictions to be satisfied often in applications.

Consult "help KronProd/methodname" for more info on each method. Optionally also, read krontest.m for demonstrations of their use.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 EXAMPLE #1:
 
 A primary application of this class is to efficiently perform separable tensorial operations, i.e., where a linear transform is applied to all columns of an array, then all rows, and so on.
 
The following example is of a separable transformation of a 3D array X that transforms all of its columns via multiplication with a non-square matrix A, then transforms all rows by multiplication with B, then finally transforms all 3rd-dimensional axes by multiplication with C. Two approaches to this are compared. The first approach uses kron(). The second uses the KronProd class. Other operations are also shown for illustration purposes.
 
Notice the orders of magnitude reduction both in CPU time and in memory consumption, using the KronProd object.
 
  
      %DATA
      m=25; n=15; p=40;
      mm=16; nn=n; pp=10;
      A=rand(mm,m); B=pi*eye(n); C=rand(pp,p);
      s=4; % a scalar
      X=rand(m,n,p);
  
  
  
  
      %METHOD I: based on kron()
      tic;
  
         Matrix = s*kron(C,kron(B,A));
  
         y1 = Matrix*X(:); %The tensorial transformation
             y1=reshape(y1,[mm,nn,pp]);
  
         z1 = Matrix.'*y1(:);
  
         w1 = Matrix.'\z1;
 
      toc;
      %Elapsed time is 78.729007 seconds.
  
  
  
  
  
      %METHOD II: based on KronProd object
      tic;
  
         Object = KronProd({A,pi,C},[1 2 3],[m,n,p],s); %equivalent to Matrix above
  
           y2 = Object*X;
                  % This operation could also have been implemented
                  % as y2=reshape( Object*X(:) , [mm,nn,pp]);
                    
           z2 = Object.'*y1;
  
           w2 = Object.'\z1;
  
      toc
      % Elapsed time is 0.003958 seconds.
  
  
 
  
      %%%ERROR ANALYSIS
      PercentError=@(x,y) norm(x(:)-y(:),2)/norm(x(:),'inf')*100;
  
          PercentError(y1,y2), % = 3.0393e-012
  
          PercentError(size(y1),size(y2)), % = 0
  
          PercentError(z1,z2), % = 1.3017e-012
          PercentError(w1,w2), % = 4.3409e-011
  
  
 
      %%%MEMORY FOOTPRINT
  
     >> whos Matrix Object
  
  
             Name Size Bytes Class Attributes
       
             Matrix 2400x15000 288000000 double
             Object 2400x15000 8102 KronProd
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
EXAMPLE #2: As a more practical example, the KronProd class is very useful in conjunction with the following tool for signal interpolation/reconstruction:

http://www.mathworks.com/matlabcentral/fileexchange/26292-regular-control-point-interpolation-matrix-with-boundary-conditions

An example involving 2D signal reconstruction using cubic B-splines is provided in the file Example2D.m at the above link.

MATLAB release MATLAB 7.9 (2009b)
Zip File Content  
Other Files
@KronProd/abs.m,
@KronProd/cellfun.m,
@KronProd/chol.m,
@KronProd/compress.m,
@KronProd/cond.m,
@KronProd/ctranspose.m,
@KronProd/display.old,
@KronProd/eig.m,
@KronProd/end.m,
@KronProd/full.m,
@KronProd/inv.m,
@KronProd/kron.m,
@KronProd/KronProd.m,
@KronProd/KronProd.old,
@KronProd/ldivide.m,
@KronProd/loadobj.m,
@KronProd/lu.m,
@KronProd/mldivide.m,
@KronProd/mpower.m,
@KronProd/mrdivide.m,
@KronProd/mtimes.m,
@KronProd/nnz.m,
@KronProd/norm.m,
@KronProd/numel.m,
@KronProd/orth.m,
@KronProd/pinv.m,
@KronProd/power.m,
@KronProd/private/allm.m,
@KronProd/private/anym.m,
@KronProd/private/blockraster.m,
@KronProd/private/cp1.m,
@KronProd/private/cp1s.m,
@KronProd/private/dimprod.m,
@KronProd/private/encell.m,
@KronProd/private/enrow.m,
@KronProd/private/equaldims.m,
@KronProd/private/horzcp.m,
@KronProd/private/is1x1.m,
@KronProd/private/iscol.m,
@KronProd/private/isnum.m,
@KronProd/private/isrow.m,
@KronProd/private/issquare.m,
@KronProd/private/isvar.m,
@KronProd/private/isvec.m,
@KronProd/private/license.txt,
@KronProd/private/parse_sizevec.m,
@KronProd/private/procShape.m,
@KronProd/private/quadformop.m,
@KronProd/private/row0s.m,
@KronProd/private/row1s.m,
@KronProd/private/rowc.m,
@KronProd/private/sameobj.m,
@KronProd/private/summ.m,
@KronProd/private/trimobjset.m,
@KronProd/private/vernum.m,
@KronProd/qr.m,
@KronProd/quadform.m,
@KronProd/rank.m,
@KronProd/rdivide.m,
@KronProd/sameobj.m,
@KronProd/size.m,
@KronProd/sparse.m,
@KronProd/subsref.m,
@KronProd/sum.m,
@KronProd/summ.m,
@KronProd/svd.m,
@KronProd/times.m,
@KronProd/transpose.m,
@KronProd/uminus.m,
@KronProd/uplus.m,
GettingStarted.txt,
krontest.m,
krontestR2006a.m,
license.txt,
UpdateNotes02222010.txt
Tags for This File  
Everyone's Tags
Tags I've Applied
Add New Tags Please login to tag files.
Comments and Ratings (3)
12 Jan 2010 Sean  
22 Feb 2010 Oleg Komarov

You may consider to edit the description since it's really unreadable. No way I'm gonna get past all those blocks of comments and those funny spaced examples.

23 Feb 2010 Matt J

Oleg, I fixed some of the ragged text justification and uncommented some of the text. Not sure if those were the unreadabilities that you were talking about, though. The Description editor is a little unwieldy and gives me very limited control of spacing.

Please login to add a comment or rating.
Updates
30 Nov 2009

Minor typo fixes in description of package

15 Dec 2009

*New class methods: times, rdivide ldivide, rank, orth, svd, pinv, qr, lu, chol, loadobj.

*Expansions of old methods: eig, mldivide, mtimes, mrdivide

*See UpdateNotes.txt in package for more

23 Dec 2009

Bug fix in qr() method (see UpdateNotes.txt). Also added relevant test of the fix to krontest.m

08 Jan 2010

Added additional example to Description section with link to a related FEX submission.

11 Jan 2010

*Overloaded cellfun() for the KronProd class

*Added simpler options for the domainsize parameter in the constructor syntax KronProd(opset,opdims,domainsizes,...)

*See UpdateNotes for details

22 Feb 2010

*Overloaded kron() for the KronProd class.
*Reimplemented the class definition using classdef file.
*loadobj now always does a construct-on-load.
*See UpdateNotes for more details.

23 Feb 2010

At the suggestion of a user, neatened up the Description section

Tag Activity for this File
Tag Applied By Date/Time
kronecker Matt J 30 Nov 2009 12:17:09
oop Matt J 30 Nov 2009 12:17:09
tensor Matt J 30 Nov 2009 12:17:09
class Matt J 18 Dec 2009 10:39:13
data type Matt J 18 Dec 2009 10:39:13
compact Matt J 18 Dec 2009 10:40:06

Contact us at files@mathworks.com