martinResearch / MatlabAutoDiff

A matlab implementation of forward automatic differentiation with operator overloading and sparse jacobians
BSD 3-Clause "New" or "Revised" License
39 stars 12 forks source link

pagemtimes #7

Closed jamesascham closed 2 years ago

jamesascham commented 2 years ago

Hello,

Thanks again for this great resource. I've managed to get it up and running with a rather large/complex algorithm. My J calcs are now an order of magnitude faster.

Is there the ability to add a pagemtimes method? I'm using multiprod (Paolo de Leva) currently which has an order of magnitude computational penalty.

Kind regards, James

martinResearch commented 2 years ago

You're welcome. Pagemtimes could be added. It should take a couple or hours to do so, but the implementation might not be very efficient in term of memory. What are the sizes of your arrays ?

martinResearch commented 2 years ago

I added support for pagemtimes. It took me a bit more time than initially expected. I don't have this function in my version of matlab (R2018) so I took multiprod to backport it (assuming the interface is the same). I don't support all different cases in pagemtimes (I am not doing the broadcasting). Hopefully that is enough for your use case.

jamesascham commented 2 years ago

Hi Martin,

Thanks so much!

Although, many of my cases are where ndims(x)~=ndims(y).

Both of: A = ones(3,3); B = ones(3,3,10); AB = pagemtimes(A,B); BA = pagemtimes(B,A);

I.e. each page of B multiplied by A, or A multiplied by each page of B.

Kind regards, James

jamesascham commented 2 years ago

It also seems that if A and B are both 3D AutoDiff instances, the code fails on line, v=repmat(reshape(y,1,size_y(1),size_y(2),prod(size_y(3:end))),size_x(1),1,1,1); as reshape(y,1,size_y(1),size_y(2),prod(size_y(3:end))) is a 4-D AutoDiff input to repmat.

martinResearch commented 2 years ago

Hi James, I added support for ndims(x)~=ndims(y) and fixed the problem when both inputs are AutoDiff instances. Please let me know if that works for you now.

jamesascham commented 2 years ago

Thanks Martin! There is a typo in line 604, should read y = y.* ones(new_size_y);

Otherwise it works perfectly!