Re: Algorithmic Differentiation in Octave

[Brad Bell]

On 01/27/2017 08:41 AM, Olaf Till wrote:
> On Fri, Jan 27, 2017 at 04:46:22AM -0700, Brad Bell wrote:
> > On 01/25/2017 01:59 AM, Olaf Till wrote:
> > > After a look at the source, it seems to me that the interface to cppad
> > > only overloads scalar operations, not matrix multiplication and not
> > > element-wise matrix operations (like e.g. 'matrix1 ./ matrix2' in
> > > Octave). This would be a rather severe limitation for application in
> > > Octave...
> > >
> > > Olaf
> > Does John's method for extending to matrix operations; see
> > http://lists.gnu.org/archive/html/octave-maintainers/2017-01/msg00235.html
> > solve your objection above ?
> Maybe I see it wrong, but I had seen JWEs comment more as an
> explanation to you of the way(s) used in Octave to overload e.g. a
> matrix multiplication, not as a complete sketch of extending your
> interface to matrix operations...
>
> (To avoid misunderstandings: Octave performs matrix multiplications
> ('a * b') or element-wise matrix operations (as 'a .* b') already for
> its builtin matrix types, without any overloading.)
>
> Your interface obviously already creates the type a_double (in a
> different way?). Do you mean to change your code to make swig produce
> m-files similar the one shown by JWE?
It is my intention to have a developer, for each target language, who 
creates a natural interface for that target language. See the heading 
Purpose on
    http://www.seanet.com/~bradbell/cppad_swig/cppad_swig.htm
I would provide any needed support and additional connections to the 
Cppad in order to make this possible and fast.
> The core problem still is how to make your code record matrix
> multiplications and element-wise matrix operations. Even if it should
> be slow only at recording time, I wouldn't like a solution which
> involves splitting into element operations at the Octave interpreter
> level. As for the other solution mentioned by you (at the swig code
> level with C-code) the details are not obvious to me either. But I
> never used swig and am not familiar with it.
Given that Octave, like python with numpy, has a way to create matrix 
operations from scalar ones, I do not think the recording will be slow.

The solutions, with atomic functions, and matrix AD operations are much 
more complicated. They may give some benefit in the long run, but I do 
not think they are a good idea for a first implementation.

> In any case, whether the splitting into element operations takes place
> in Octave or in the swig interface, the result, after recording and
> running, must be translated back into Jakobians and Hessians in the
> matrix- or array-representatin of Octave. This is an additional
> problem, which AFAICS can't be reasonably solved. So maybe it would be
> more reasonable to support matrix operations directly in the cppad
> library?
I do not see any problem here. The Jacobians and Hessians are returned 
as simple vectors and special purpose *.m files  could put them in any 
desired form.

> But all these matrix issues could be implemented in a simple forward
> method in pure interpreted Octave code in a rather straightforward
> way... and consider that probably in most large user functions the
> path of computation depends on the parameters...
>
> Olaf
Forward, tapeless mode, has advantages. In the simple case, of only 
needing first order derivatives, the complex step method provides this 
with no extra work (in Octave) and is probably as fast as possible.

On the other hand, optimizing the tape, computing sparsity patterns and 
sparse derivatives, reverse mode, and just function evaluation with C++ 
speed, has advantages also. For example, reverse mode can compute the 
derivative of a scalar value function with respect to an arbitrary 
number of variables in a small multiple of the number of floating 
operations for the original function.