File Exchange Pick of the Week

Our best user submissions

Efficient Object-Oriented Kronecker Product Manipulation

Sean's pick this week is KronProd by Matt J.

Background

This one is old and very dear to my heart. Some time in the 2009 era, as I was starting grad school, I needed help calculating subpixel resolution displacement vectors between two CT image volumes. The displacement vector calculation was one of the main underlying pieces of most of my grad research. This was before imregcorr existed, though it still does not support 3D which was my requirement.
I'd managed to get a 2D variant working, but needed 3D. Matt replied to my newsreader post to discuss the algorithm and its extension in 3D, and pointed to this to be more efficient than using kron directly. Performance was important as I'd have to compute over four billion displacement vectors over the next year.
Matt has provided many algorithms and overloaded a bunch of operators to work with the KronProd class and configured it with custom indexing and size. In 2009, I didn't understand how custom indexing with subsref/subsasgn/size overloading worked, so it seemed like magic. With the knowledge of how those work now, and the performance implications in that era, I probably would've pulled the core algorithm into a function for speed since my use was stateless, i.e. a new KronProd each time.
Here are a few simple operations with a KronProd. We'll create one from a three arrays.
kp = KronProd({magic(3), [randi(10,[3 4]), zeros(3,1)], randn(9, 15)})
kp =
81×225 KronProd array with properties: opset: {[3×3 double] [3×5 double] [9×15 double]} opinds: [1 2 3] numops: 3 eyemask: [0 0 0] domainset: [3 5 15] maxdim: 3 scalarcoeff: 1 scalarcumprods: [1 1 1] scalarmask: [0 0 0] domainsizes: [3 5 15] rangesizes: [3 3 9]
Now we can work with this as if the full expansion you'd get from the kron function.
nzeros = nnz(kp)
nzeros = 14580
Or a matrix operation where I get back KronProds.
[u, s, v] = svd(kp)
u =
81×81 KronProd array with properties: opset: {3×1 cell} opinds: [1 2 3] numops: 3 eyemask: [0 0 0] domainset: [3×1 double] maxdim: 3 scalarcoeff: 1 scalarcumprods: [1 1 1] scalarmask: [0 0 0] domainsizes: [3 3 9] rangesizes: [3 3 9]
s =
81×225 KronProd array with properties: opset: {3×1 cell} opinds: [1 2 3] numops: 3 eyemask: [0 0 0] domainset: [3×1 double] maxdim: 3 scalarcoeff: 1 scalarcumprods: [1 1 1] scalarmask: [0 0 0] domainsizes: [3 5 15] rangesizes: [3 3 9]
v =
225×225 KronProd array with properties: opset: {3×1 cell} opinds: [1 2 3] numops: 3 eyemask: [0 0 0] domainset: [3×1 double] maxdim: 3 scalarcoeff: 1 scalarcumprods: [1 1 1] scalarmask: [0 0 0] domainsizes: [3 5 15] rangesizes: [3 5 15]
Matt has overloaded subsref so when you index into the KronProd, you get back a KronProd or just matrix.
kkp = kp(1)
kkp =
3×3 KronProd array with properties: opset: {[3×3 double]} opinds: 1 numops: 1 eyemask: 0 domainset: 3 maxdim: 1 scalarcoeff: 1 scalarcumprods: 1 scalarmask: 0 domainsizes: 3 rangesizes: 3
km = kp{1}
km = 3×3
8 1 6 3 5 7 4 9 2
He did not overload subsasgn though, so we can't do the opposite.
try
kp{1} = rand(3)
catch ME
fprintf(2, ME.message)
end
Unable to perform assignment because brace indexing is not supported for variables of this type.
He's also overloaded size so the appearance of size is that of the full Kronecker product. Under the hood, however, the KronProd is actually a scalar, it just appears to be the expanded size.
sz = size(s)
sz = 1×2
81 225
Starting in MATLAB R2021b, we now have an easier and more formal means for controlling indexing and creating scalar objects using modular indexing.
Right now subsref overloads (), {}, and dot(.). It overloads dot, just because it has to with the old approach, so we don't need to do that with the newer indexing. So the class definition line would need to be:
classdef KronProd < matlab.mixin.indexing.RedefinesParen & matlab.mixin.indexing.RedefinesBrace
Now, we'd need to implement the abstract methods for parenAssign, parenReference, braceAssign, braceReference. These methods would then use similar logic to what's used today, in a protected methods block, something along these lines:
function B = braceReference(obj, indexOp)
inds = M.opinds(indexOp(1).Indices{1});
B=M.opset(inds);
if isnum(inds)
B=B{:};
end
end

Comments

Matt, if we ever meet in real life, I owe you a beer!
Give it a try and let us know what you think here or leave a comment for Matt.
|
  • print

Comments

To leave a comment, please click here to sign in to your MathWorks Account or create a new one.