mtimes and transpose are being treated as a single function when used together?

I'm trying to overload the mtimes function so each * operator will perform my own matrix multiplication (using mex) instead of matlab's built-in function.

I managed to overload the mtimes function when writing it in the same m-file with the calling function, and when I perform A*B I see that my function is being called (using debug message running from my mtimes function). So far so good. But when I perform A*B.' my mtimes is NOT called!. Neither when I perform A.'*B or A.'*B.'.

See the example code below, and the debug messages I see in Matlab console when I run it.

function testme
A = randn(3,3);
B = randn(3,3);

disp('start');
C1 = A*B;
C2 = A*B.';
C3 = A.'*B;
C4 = A.'*B.';
disp('end');
end

function out = mtimes(a,b)
out = a + b; % here should come my mex function. using '+' to let others copy-paste and check it for themselves without the mex being missing.
disp('been in my-mtimes');
end

output on console:

>> testme
start
been in my-mtimes
end

It looks like A*B.' is calling different functions than 'mtimes' following 'transpose'... Any idea which function is called on A*B.' or A.'*B? BTW wrapping B.' in brackets (B.') doesn't change the behavior. Neither using A*transpose(B)... None of them are entering my mtimes function... I guess this is some kind of optimization that matlab does to save the transpose MIPS, but there must be a function that does that, which I can't find... Any help would be appreciated!

Thanks, Koby

[Bruno Luong]

Yes MATLAB optimizes in the product by using the right LAPACK/BLAS (?) without explicit transpose of the arguments.

Can you avoid the optimization by split in two half

Bt = B.';
C = A*Bt;

[Koby Goldberg]

Thanks Bruno, but splitting the line is not an option to me, because I must not change the code... That's why I was looking for a way to overload the * operator, so the original code would be the same, but the actual multiplication function would be my new mex function...

If mtimes is not called in this case, which function does? If I knew, I could overload that one as well with a mex function.

[Matt J]

If mtimes is not called in this case, which function does? If I knew, I could overload that one as well with a mex function.

Even if you knew, I don't think the interpreter would let you shadow it. The fundamental data types in Matlab have rules all their own.

[Bruno Luong]

I'm afraid those optimization are done under the hood and no equivalent function form exists.

[Bruno Luong]

"Neither using A*transpose(B)... None of them are entering my mtimes function..."

If you put the statement in a script it will however call MTIMES. Reason: The script does not have the jit accelerator that does this kind of hidden optimization. I did not check but might be if you put a relevant code in a script and call it with RUN command, even within the function, you might avoid the JIT ACC to kick in.

[Walter Roberson]

This was a good question; this is behavior that is not clearly documented in the class hierarchy charts.

[Matt J]

BTW wrapping B.' in brackets (B.') doesn't change the behavior. Neither using A*transpose(B)... None of them are entering my mtimes function...

But this works: A*(B.'+0)

[James Tursa]

It all depends on what syntaxes the parser recognizes in order to call the combined BLAS function. I have never seen any documentation on this, nor have I done any testing to see how this might have changed from release to release.

[Koby Goldberg]

Thanks Matt, nice trick with adding 0 ! I verified, it indeed works :)

Thanks!

Answer 1

[Matt J]

To circumvent this, my recommendation would be to overload mtimes in a subclass of type double.

classdef myclass < double                   

        methods
          function obj=myclass(data)      

              obj=obj@double(data);

          end

          function out = mtimes(a,b)
              out = myclass( a + b );
              disp('been in my-mtimes');
          end

        end

end

[Koby Goldberg]

Thanks Matt, that's an interesting answer. I guess when performing transpose on those objects would perform the Matlab's built-in transpose function, because this class inherits from double.

The only drawback is that I have to define all the matrices as this class type upfront. And I have tons of matrices...

So if possible, I would like to return to my original question: which function runs when A*B.' is performed?

[Matt J]

So if possible, I would like to return to my original question: which function runs when A*B.' is performed?

I don't know, but I'm 99.99% certain that it is not accessible.

Why does it matter that you have to pre-convert lots of input matrices to myclass? Put them in a cell array and do a loop, or call cellfun,

argumentsMyClass =  cellfun(@myclass, arguments,'uni',0 );

[James Tursa]

"which function runs when A*B.' is performed?"

For full floating point matrices, a BLAS library function that does the multiply and transpose all in one operation. The transpose is not physically formed first, rather the indexing used in the multiply takes care of the transpose operation implicitly. The library code is highly optimized for your machine and is multi-threaded. The specific function called depends on the class and complexity of A and B, whether or not B is the same matrix as A, whether these are 2D matrices or 1D vectors, etc. One can write mex code that links with these libraries and calls these functions directly (there are several examples of this in the FEX).

The library routine called could be DGEMM or DSYRK or CGEMM or CSYRK or ...

The MATLAB parser needs to recognize the syntax for the classes involved in order for these combined operation functions to get called. Otherwise the transpose is explicitly performed first followed by the matrix multiply. In addition to A*B.', it recognizes A*B', A'*B, A.'*B, etc.

If sparse matrices are involved, then different library code gets called.

[Koby Goldberg]

Thanks guys for all the info and good suggestions. So, bottom line, I guess it is impossible to override the *' pair, since this optimization is performed by the interpreter under the hood. Anyway, you've given me good workaround ideas!

Thanks, Koby

[Matt J]

You're welcome, but please Accept-click the answer if you consider the matter addressed.