SymPy: lambdified dot() with (3,n)-array

Question

I have a lambdifyd sympy function with a dot product in it, e.g.,

import numpy as np
import sympy


class dot(sympy.Function):
    pass


x = sympy.Symbol('x')
a = sympy.Matrix([1, 1, 1])
f = dot(x, a)

ff = sympy.lambdify((x), f, modules='numpy')

x = np.random.rand(3)
print(ff(x))  # okay

(Weirdly, the custom dot declaration works. Don't know why exactly, but never mind that. If there's a better solution, let me know.)

I'd now like to call ff with a bunch of vectors at once, so I went

x = np.random.rand(3, 10)
print(ff(x))

No good!

ValueError: shapes (3,10) and (3,1) not aligned: 10 (dim 1) != 3 (dim 0)

Alright, so I'll somehow have to transpose the first argument of dot. Using the good old .T on sympy.Symbol('x') isn't legal though.

Any hints on how to do mass dot-products from a lambdified sympy expression?

Answer 1

You're doing multiple odd things, but I can't tell how much of this is due to an oversimplified MCVE.

First, a bit more elegant definition of your function:

import sympy as sym    

x = sym.Symbol('x')
a = sym.Matrix([1, 1, 1])
dot = sym.Function('dot')
f = dot(x, a)

ff = sym.lambdify(x, f, modules='numpy')

The reason this and your original kludge works is because all you have to achieve is to have something that says "dot". Once you have that, lambdify will substitute np.dot to that part of your symbol.

Now, just for completeness, here's how I would have done this instead:

import numpy as np

a = np.array([[1],[1],[1]])
ff = lambda x,a=a: np.dot(x,a)

I know this might not be an option in your actual problem, but my experience is that if something can be done without symbolic math, then it's worth doing it like that.

Now, for your error. The error is pretty clear, and so is the math. You defined a function which for any input x computes x*a with a 3d column vector a. As the error suggests, this makes sense in a very limited number of cases. It would make sense if both operands were 3-element 1d arrays, in this case the scalar product would be returned. However, since one of your operands is fixed to shape (3,1), np.dot only performs matrix multiplication (and for vector input, returns a 1-element 1d array instead of a scalar). Due to your definition, it only works with matrices that can be multiplied from the right by a, i.e. matrices of shape (N,3). Clearly this is not the case for your input.

What you should do, is transpose x on the numerical side:

x = np.random.rand(3,10)
print(ff(x.T))

This will enter an array of shape (10,3) into the function, which then multiplies it by one of shape (3,1), resulting in a 2d array of shape (10,1): a column vector, each row containing the scalar product if the given input vector with a.

The other option is to swap the definition of your function:

f = dot(a.T,x)
ff = sym.lambdify(x, f, modules='numpy')
# or
a = np.array(1,1,1)  # transpose of previous a
ff = lambda x,a=a: np.dot(a,x)

Both will create a function that multiply an array of shape (1,3) from the right with the input. Then your input x of shape (3,10) is directly compatible; the output will be a 1d array of the 10 scalar products. In this formulation your