Indeed, I'm happy to report that the numerical errors were not due to your package, after all.
Interestingly, the errors I was having were not due to MMX giving me the pseudo-inverse solution (which is what I wanted), but rather they were due to MATLAB's pinv() function, whose output I passed to MMX. In particular, I found out that MATLAB does not use symmetric-matrix-specific algorithms for functions such as pinv, svd, or sqrt, which may result in output matrices having non-zero imaginary terms. Of course, this did mess up the MMX output (as noted earlier), but the input was wrong to begin with.
After going around that problem, MMX worked like a charm!
It appears to me that your package suffers from numerical artifacts, at least with respect to 'BACKSLASH' computations.
To check, I ran something like this:
% A:3x3x1850, B:3x1x1850
C = mmx( 'backslash', A, B );
D = zeros( size( C ) );
for k = 1 : size(D,3)
D(:,:,k) = A(:,:,k) \ B(:,:,k);
E = sum( (C - D).^2, 1 );
stem( permute( E, [1 3 2] ) )
It is readily visible that E is zero everywhere except for 4-5 indices, where the error is pretty large.
Unfortunately, I haven't been able to consruct a MWE with rand() data or something similar. Please PM me if you want me to send you the data that gave me the errors.
I hope that you have the time to work on this issue (assuming that I haven't done something silly), as I think that this is a very handy package to have! I'd be willing to help if I can.