File Exchange

image thumbnail

Efficient Object-Oriented Kronecker Product Manipulation

version 1.10.0.0 (55.3 KB) by Matt J
A class for efficient manipulation of N-fold Kronecker products in terms of their operands only.

9 Downloads

Updated 05 Aug 2010

View License

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.

Cite As

Matt J (2020). Efficient Object-Oriented Kronecker Product Manipulation (https://www.mathworks.com/matlabcentral/fileexchange/25969-efficient-object-oriented-kronecker-product-manipulation), MATLAB Central File Exchange. Retrieved .

Comments and Ratings (22)

Matt J

Hi Cara,
The reason for that is because the element-wise multiplication of two Kronecker products results in a matrix that is also a Kronecker product, provided the matrices involved are the correct size. However, this is not true for either addition or the log10() operation. They do not preserve Kronecker product structure.

Whether you can gain anything from KronProd for your application depends on what you plan to do with the resulting matrix log10(kron(B,A)). For example, if you plan to multiply with an ordinary vector, x, you can use do so in the factored form,
A1 = KronProd({log10(A),ones(size(B))},[1,2]);
B1 = KronProd({ones(size(A)),log10(B)},[1,2]);
Cx = A1*x+B1*x;

Hi Matt,
I have an application where I need to take log10 after the kronecker product of matrices A and B. Originally, I tried:

C1 = KronProd({A,B},[1,2]);
C2 = log10(C1);

and got the error that log10 is not a defined function for inputs of type KronProd. Instead, I tried using the rule that log(A*B) = log(A) + log(B):

A1 = KronProd({log10(A),ones(size(B))},[1,2]);
B1 = KronProd({ones(size(A)),log10(B)},[1,2]);
C1 = A1+B1;

but plus / + is also not supported by KronProd. I'm surprised there is no support for element-wise addition as there is for element-wise multiplication. Unless I am missing something?

Thanks for your time.

Matt J

Hi AP,

Sorry I didn't respond sooner. For some reason, the File Exchange doesn't send prompt notifications anymore...

When using KronProd operators, it is important to read the help doc to be aware of some of the restrictions.

>> help KronProd.times
times method for KronProd.

Y=times(M,X) or Y=M.*X

At least one of M or X is a KronProd object. The other may be a scalar in which
case Y is an appropriately scaled KronProd object.

If M and X are both KronProd objects, then the operation is attempted on corresponding
operands of M.opset and X.opset, and will only work if the sizes of the operands
are compatible. E.g., if

M1 represents (A kron B)
M2 represents (C kron D)

then M.*X tries to construct a KronProd representing (A.*C kron B.*D). If
the operands are sized compatibly with this, full(Y) is equivalent to full(M).*full(X).

****

In your case, the operation doesn't work because the matrices in {ones(1,p),A,ones(N,1)} are not the same size as the corresponding matrices in {ones(q,1),X',ones(1,r).

Here is an example of a supported .* operation with KronProd objects:

>> Object = KronProd({rand(3), rand(4)}) .* KronProd({ones(3), randn(4)})

AP

Hi Matt

The KronProd object makes me surprised. Now I am doing with the multiplication of two matrices from Kronecker product in a large-scale setting but I am not sure that I make something wrong. My KronProd object cannot multiply with other matrices. My example is shown as follows:

n=5; N=200; p=18*n; q=7*n; r=7*n;
H=kron(kron(ones(N,1),A),ones(1,p)).*repmat(kron(X',ones(q,1)),1,r); %this is what i need

but what I try is that

Object1 = KronProd({ones(1,p),A,ones(N,1)},[1,2,3]);
Object2 = KronProd({ones(q,1),X',ones(1,r)},[1,2,3]);
Object3 = Object1.*Object2;

I hope to see that Object3 will be equal to H but it has some errors.

"A.*B where A and B are both KronProds requires that A and B have operands of matching sizes"

So I still try another way by

Object3 = Object1.*.*repmat(kron(X',ones(q,1)),1,r)

But it still has some errors,

"Undefined times operation involving KronProd"

In this case, do you have any comments about creating a matrix H using KronProd? I may make something wrong but I cannot solve it. Thank you in advance.

Matt J

Hi Elad,

Sorry I didn't see this sooner. As you mentioned, full(Object) is the way to obtain the full array form of the Kronecker product. It might also make sense to do sparse(Object), if your input operands are themselves sparse matrices. You are correct, however, that doing so would defeat the purpose and benefits of the object representation. I would point out, however, that KronProd objects are capable of doing other things besides multiplication, see or example,

>> methods KronProd

So, if you elaborate on what you're trying to do, I might b able to make some recommendations.

Elad Harel

Is it possible to extract the full matrix or is the speedup only occur when multiplying the output? For instance, I can use full(Object), but this really slows down the code. I suppose the speedup arises because this attempts to take advantage of the sparsity of the data and if the entire tensor is extracted those savings are lost, but I want to make sure I'm not doing something wrong. To be explicit, I want the full array kron(C,kron(B,A)) and not just the product of this array with some other matrix (e.g. X in your example).

Matt J

@bazrafshan,

Presumably you would do something like

ifft( conv(fft(x),fft(y)))

However, I don't see the computational motivation of doing so. x.*y is surely the fastest way to implement x.*y. I also don't see the connection between this question and KronProd.

Hi Matt

How does it work for this case?

Suppose we have the following matrices:

x = [3 5 4; 7 6 1; -1 2 0];

y = [2 7 1; 2 -3 2; 5 6 9];

How can I get the same results as x.*y using convolution in frequency domain?

Thanks in advance

Matt J

Well, if all else fails, you can always use pcg(). The KronProd representation will facilitate the matrix multiplication operations.

Manny

Thanks for your answer, Matt.

The matrix F with which I am working is the transition matrix of state-space model:

x_{t+1} = F*x_t + e_t,

so I am not sure it is always diagonalizable. The problem with which I am working is indeed a system of equations:

(I-kron(F,F))*vec(P) = vec(Q),

where F has dimension close to 350x350.

Matt J

Manny,

There may be options, but it depends on the specifics of F and what exactly you're trying to do with the inverse.

For example, if F is diagonalizable (e.g., symmetric or possessing distinct eigenvalues) and if you are really trying to solve the equations

(I-kron(F,F))*x=y

for x, then I might do the following,

K=KronProd({F},[1,1]);

[V,D]=eig(K);
d=1-full(cellfun(@diag,D));

x = V*((V\y)./d);

Manny

I need to obtain inv(I-kron(F,F)), where F is a large matrix. Is there any script in the codes that could allow an efficient treatment of that (as an object)?

Manny

WALEEM

Thanks for sharing this code, Matt. With the help of your code, I have been able to resolve an estimation issue I have had (for the past 13 months) with a model. In the process of estimating the model, I had to deal with a 4^10-by-4^10 (8796Gb) matrix in a several thousands loop. This almost meant the end of the project.

But I stumbled upon your code, and was able to complete the construction of the model.

I can't thank you enough for this. I will be sharing my code when it is all done, and I will duely reference your code.

Thank you.

Matt J

Hi Roemer,

Thanks for you feedback. Regarding concatenating KronProd objects: The concatenation operation that you've described isn't possible precisely as you've requested it, because the concatenation of 2 Kronecker products is not, in general, mathematically equivalent to a 3rd Kronecker product. Because it's not possible mathematically, it is equally not possible in software.

However, I did create the following FEX package to allow one to combine simpler operators into more complicated operators, which might be equivalent to what you're asking for

http://www.mathworks.com/matlabcentral/fileexchange/26611-on-the-fly-definition-of-custom-matrix-objects

Here's a simplified example where I concatenate kron(eye(3),eye(3)) with itself

>> A=eye(3); B=eye(3); K=KronProd({A,B});
>> Kcat=MatrixObj; Kcat.Ops.mtimes=@(A,x)[K*x;K*x];
>> Kcat*(1:9)'

ans =

1
2
3
4
5
6
7
8
9
1
2
3
4
5
6
7
8
9

The mixture of the MatrixObj class with other classes has had some problems in past MATLAB versions, unfortunately, so how well it works for you could depends on what version you have.

This is a fantastic piece of work. The amount of speed-up and memory efficiency has been a tremendous help already. There are calculations I have to do on often enormous arrays resulting from Kronecker products, easily needing up to 400GB of RAM/memory-mapped space. Using your OOP implementation is of huge help to me and my (PhD) project, and it's a nice way to learn how to do some OOP myself. Thanks!

I have one question/request, I assume you might not have the time, but maybe you could give me some pointers on how to do it. Maybe it's very easy, and I'm just missing it.
I would like to be able to do:

KronProdC = cat(1,KronProdA,KronProdB)

I.e., concatenating KronProd objects along dimension 1 (or 2). Concatenation would occur over 1st/2nd dimension of the kronecker product. Maybe it's not necessary to expand the class, but use a trick to modify existing KronProd objects, but my intuition for this is not very trained.
Being able to do this, I could extend your work into a Kathri-Rao product, which would speed up a large set of problems and algorithms using the PARAFAC/TuckerX models. These multi-dimensional decomposition models are used a lot in chemometrics, and are gaining a lot of popularity in neuroscience the past years (my field). Both usually use Matlab for these calculations, and this would speedup results for a lot of people, many could benefit. So, I got very excited when I saw and tested your work :).
Any help is greatly appreciated,
Roemer

Danielle

Wow! No more out of memory errors and so fast!

Omar

This is what I need :-)

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.

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.

Sean

Updates

1.10.0.0

In constructor KronProd(opset,opinds, ...) introduced default value for opinds of 1:length(opset).

This allows simple 1-argument contructor calls such as
KronProd({C,B,A}), which is equivalent to kron(A,kron(B,C))

1.9.0.0

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

1.8.0.0

*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.

1.7.0.0

*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

1.6.0.0

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

1.5.0.0

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

1.3.0.0

*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

1.1.0.0

Minor typo fixes in description of package

MATLAB Release Compatibility
Created with R2009b
Compatible with any release
Platform Compatibility
Windows macOS Linux