# Faster version of num2cell overloaded assignment

7 views (last 30 days)
Nathan Fitzgerald on 28 Mar 2019
Answered: Nathan Fitzgerald on 28 Mar 2019
I'm trying to write a class that overloads basic mathmatical operators without needing to overload every array manipulation function matlab has to offer. Are there any basic class archtetures that accomplish this? Please forgive the long-winded lead-in description below.
There are a few examples in the Matlab File Exchange of classes that attempt to implment automatic differentiation. The general concept is that you instantiate an object of a class that you can pass through a calculation that normally operates on doubles, and have that object automatically track the derivatives through the chain rule. The high level syntax is something like
z = myCalculation(x); % normal calculation where we would want to track derivatives.
% Assume x is a column vector like [1;2;3;4]
zAutoDiff = myCalculation( AutoDiff(x) ); % run the calculation again with derivative tracking
disp( jacobian(zAutoDiff) ); % the AutoDiff class has a method called jacobian that reports dzdx.
% If numel(x) == numel(z) == 4, then dzdx is a 4x4 matrix
The approach to this taken for the autodifferentiation examples in the file exchange (that I'm aware of) generally use the paragdim
classdef AutoDiff
properties
val
der
end
methods
function self = AutoDiff(x,d)
self.val = x
if nargin > 1
self.der = d
end
end
end
end
and then you overload the basic operator methods like "plus", "mtimes", "uminus", etc to update the "val" and "der" properties appropriately as they pass through a calculation.
At this point, I see a choice between 2 data storage architectures.
1. Instantiate a 1x1 AutoDiff object where "val" and "der" are matrices that contain all of the info. This requires the overloading of a TON of basic addressing functions like ctranspose, transpose, vertcat, numel, length, as well as subsref and subsasgn if you want statements like y(3) to return or operate on the 3rd element of val, not the third element of our 1x1 AutoDiff object (which, we've specified has only having 1 element in this storage approach). Forgetting to implement a function like permute leads to an error that will not cause an exception. The code will seem to run fine but the answer will be wrong, pretty dangerous and super hard to debug.
2. Instantiate the AutoDiff object as having the same size as the input "x", such that the "val" property is a scalar for every individual element of the object. This means that methods like transpose, permute, numel, size, etc do not have to be overloaded. However, manipulation of data in the class would need to use num2cell and comma separated lists to get things back into the structure:
tmp = num2cell(x); % x is a vector whose elements I want as individual AutoDiff.val scalars
[z(:).val] = tmp{:} % This works but is SLOW for large vectors
Is there some other template for creating this sort of class that:
• Doesn't require the class designer to rewrite basic array manipulation functions. I don't see way around this when the base object is 1x1, but we're trying to make it look like it has the size of the mxn internal property using a bunch of operator overloading.
• Doesn't require the increadibly slow num2cell to parse data into a comma-separated list? In my tests, the comma-separated assignment doesn't seem to be the bottleneck, just the "tmp = num2cell(x)" step. Can the [z(:).val] assignment happen a litte more directly?

Nathan Fitzgerald on 28 Mar 2019
I think I've come to an answer that I'm happy with. As I postulated in the very last question of the original post, there is a way that the [z(:).val] assignment can happen more directly. Instead of using num2cell to break x up into individual cells that could then be used to create a comma-separated list (e.g. below),
tmp = num2cell(x); % x is a vector whose elements I want as individual AutoDiff.val scalars
[z(:).val] = tmp{:} % This works but is SLOW for large vectors
I wrote a mex function that directly took the inputs x and dealt them to individual outputs
[self(1:numel(x)).val] = fastDealScalar(x);
where fastDealScalar.c is
#include "mex.h"
/* Output mxArray's should all be 1x1 */
#define ROWS 1
#define COLUMNS 1
#define ELEMENTS 1
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
mwSize index;
mwSize nElems;
double *y, *z;
/* define the number of elements in the input matrix */
y = mxGetPr(prhs[0]);
nElems = mxGetNumberOfElements(prhs[0]);
/* loop through the element of the input and assign them
* individually to the output */
for ( index = 0; index < nElems; index++ ) {
/* set the output pointer for this element */
plhs[index] = mxCreateNumericMatrix(ROWS, COLUMNS, mxDOUBLE_CLASS, mxREAL);
/* create a C pointer to a copy of the output matrix */
z = mxGetPr(plhs[index]);
/* copy a scalar from the input matrix to the output */
*(z) = *(y+index);
}
return;
}
(This was hacked together from matlab example mex files, so forgive the coding style.... I'm not very experienced in C)
So that bypasses the "num2cell" bottleneck and now my code seems to be performing well.