In Matlab 2017b Mathworks introduced the decomposition class, which precomputes a decomposition of a linear system:
>> A = rand(10000, 10000); >> b = rand(10000, 1); >> tic; x = A \ b; toc Elapsed time is 5.049433 seconds. >> tic; dA = decomposition(A); toc Elapsed time is 4.916148 seconds. >> tic; x2 = dA \ b; toc Elapsed time is 0.039961 seconds. >> nnz(x ~= x2) ans = 0
As can be seen, once the decomposition object is created, this allows solving the system at least two orders of magnitude faster with the exact same results.
The decomposition class is as smart as mldivide (backslash operator for solving Ax = b) and mrdivide (forward slash operator for solving xA = b) and follows the same logic as can be seen in the flow charts for dense or sparse matrices. That means it recognizes special kinds of system matrices and performs a suitable decomposition that exploits the matrix structure.
An idea that quickly comes to mind is using a decomposition object inside of a parfor or e.g. functions like blockproc with 'UseParallel' enabled.
When attempting to do so, you will quickly be disappointed:
Warning: Loading a decomposition is not supported. The loaded object is the decomposition of an empty matrix. > In decomposition.loadobj (line 673) In parallel.internal.pool.deserialize (line 29) In parallel.internal.pool.deserializeFunction (line 12) In remoteParallelFunction (line 19)
and some subsequent errors like
Error using \ (line 391) Matrix dimensions must agree.
Matlab's parallel pool works by replicating everything in the workspace for each of the workers. This replication is performed by serialization followed by an unserialization inside of each worker, i.e. converting all objects into a bytestream as if they were written to disk in a mat file. In fact, those are the same methods that are used when calling save() or load().
For classes, it is possible to overload the methods saveobj() and loadobj() to customize the (un)serialization. The usual pattern in these two methods is to store all properties of an object in the fields of a struct, or vice versa, to reconstruct an object from such a struct.
Unfortunately, these methods are missing implementations in the decomposition class:
methods (Hidden)
...
function s = saveobj(~)
s = struct;
warning(message('MATLAB:decomposition:SaveNotSupported'));
end
end
methods (Hidden, Static)
function obj = loadobj(~)
warning(message('MATLAB:decomposition:LoadNotSupported'));
obj = decomposition;
end
end
Luckily, it is straight forward to implement the (un)serialization for decomposition objects:
methods (Hidden)
...
function s = saveobj(obj, ~)
mc = metaclass(obj);
props = {mc.PropertyList.Name}';
s = struct;
for ii = 1 : numel(props)
s.(props{ii}) = obj.(props{ii});
end
end
end
methods (Hidden, Static)
function obj = loadobj(s)
fns = fieldnames(s);
obj = decomposition([]);
for ii = 1 : numel(fns)
obj.(fns{ii}) = s.(fns{ii});
end
end
end
Now, when attempting to load a serialized decomposition object, or when using it in parfor and the like, we get:
Warning: While loading an object of class 'decomposition': > In class 'decomposition', no set method is defined for Dependent property 'Datatype'. A Dependent property needs a set method to assign its value.
Again, this is easily fixed by adding a suitable set method:
methods
...
function f = set.Datatype(f, dt)
assert(any(strcmpi(dt, {'single', 'double'})));
f.IsSingle = strcmpi(dt, 'single');
end
...
end
You can add this method right after the get.Datatype() method. You might have to adapt the file permissions of matlab/toolbox/matlab/matfun/decomposition.m to be able to save the file.