Efficient Object-Oriented Kronecker Product Manipulation
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)})
Now we can work with this as if the full expansion you'd get from the kron function.
nzeros = nnz(kp)
Or a matrix operation where I get back KronProds.
[u, s, v] = svd(kp)
Matt has overloaded subsref so when you index into the KronProd, you get back a KronProd or just matrix.
kkp = kp(1)
km = kp{1}
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
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)
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!
Comments
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.