Extend numpy.dot() to accept more than 2 arrays. (Trac #741)

[numpy-gitbot]

Original ticket http://projects.scipy.org/numpy/ticket/741 on 2008-04-17 by trac user joris, assigned to unknown.

if a,b,c,d are numpy arrays, than their matrix multiplication can be computed with:

dot(a, dot(b, dot(c, d)))

or by

array(mat(a) * mat(b) * mat(c) * mat(d))

More convenient would be dot(a,b,c,d). An easy way to implement this would be

def dot(*args):
return reduce(olddot, args)

where 'olddot' is the current implementation of dot that only takes 2 arguments. The suggestion above would not break any code, as dot(a,b) == olddot(a,b).

[numpy-gitbot]

@rkern wrote on 2008-04-17

This has been discussed before, and some implementations were experimented with. The main problem is that the associativity of the multiplications does matter if one is concerned about memory usage. Some intermediates are (much) bigger than others. However, the syntax of dot(a, b, c, d) does not clearly show this association, and this will make things difficult for some users to identify the cause of the poor performance. If the implementation weren't a trivial one-liner, we might have included it with strong warnings. But since it is a trivial one-liner, we don't feel too bad telling people that want it to just do the one-liner in their own code.

At least that's how I remember the conversation. I can't find the thread off-hand, though.

[numpy-gitbot]

Milestone changed to 1.2 by @rkern on 2008-04-17

[numpy-gitbot]

trac user joris wrote on 2008-04-17

The conversation Robert mentioned can be found in the numpy archives on gmane, with the keywords:[[BR]]

• extended dot() (Tim Hochberg) [[BR]]
• simple multi-arg wrapper for dot (Bil Baxter, Anne Archibald, Colin J. Williams)

The associativity argument of Tim boils down to:

1. to be efficient, an extended dot should be able to choose the order of operations. E.g. if A is 100x3, B is 3x100, and C is 100x1, then dot(dot(A,B), B) is less efficient than dot(A, dot(B,C)) since the former requires an intermediate 100x100 array.

2. the optimum order of operations is difficult to choose because dot() is intrinsically not-associative because it can do both matrix and dot multiplications: e.g. dot([1,2], dot([3,2], [1,1])) != dot(dot([1,2],[3,2]), [1,1])

Bill provides an interesting suggestion: no optimization provided, but giving the user a way to specify the order of multiplication:

mdot(a,(b,c),d)
mdot(a,(b,c,d))
mdot(a,b,c,d)

etc.

with

def mdot(*args):
    if len(args) == 1:
        return args[0]
    elif len(args) == 2:
        return _mdot(args[0], args[1])
    else:
        return _mdot(args[:-1], args[-1])


def _mdot(a,b):
    if isinstance(a,tuple):
        if len(a) > 1:
            a = mdot(*a)
        else:
            a = a[0]           
    if isinstance(b, tuple):
        if len(b) > 1:
            b = mdot(*b)
        else:
            b = b[0]

    return N.dot(a,b)

However, in its current form mdot(a,b,c,d) is about twice as slow as dot(a,dot(b,dot(c,d))), which is kind of show-stopper.

[numpy-gitbot]

@stefanv wrote on 2008-04-17

Is the optimisation prohibitively expensive for, say, anything up to 5 elements?

[numpy-gitbot]

@rkern wrote on 2008-04-17

I would say that the ambiguity would be the deal-breaker. There's no point in making it fast if it's wrong.

[numpy-gitbot]

@stefanv wrote on 2008-04-17

Right, because dot does broadcasting as well as matrix multiplication.

[numpy-gitbot]

Milestone changed to 1.3.0 by @cournape on 2008-08-13

[numpy-gitbot]

Milestone changed to Unscheduled by @cournape on 2009-03-02

[numpy-gitbot]

@mwiebe wrote on 2011-03-23

The "dot as a method" introduced in 1.5 improves the syntax available for many dots in a row, so closing the bug as wontfix.

The optimization ideas discussed might be served well with the deferred expression evaluation idea. If dot returns an unevaluated expression tree instead of an evaluated result, a complicated end-result expression can be analyzed for reordering that reduce intermediate value sizes, etc.