Replacing a scalar/symbol by an array

[RPKrijnen]

Dear community,

I would like to replace all the symbols in a given expression by arrays all of the same length. Doing this would basically turn a scalar equation into a vector valued equation. As far as I know, simply replacing all those symbols by corresponding array symbols isn't allowed.

Now the reason I would even want to do this is because I need to evaluate some expression a large number of times. In order to make this computation more efficient I figured I could substitute in all my constants (measurements) at once and lambdify this expression, or ideally transform it into a form acceptable by PyTorch.

Here is an example of the code:

import sympy as sm

a1, b1 = sm.symbols('a1 b1')
F = sm.MatrixSymbol('F', 3, 3)
C = F.T*F
I1 = sm.Matrix(C).trace()
W = a1/(2*b1)*(I1**(-1/3))
P = sm.derive_by_array(W, sm.Matrix(F))

This yields expressions of W and P that are dependent on the individual components of F. The thing is that I have multiple evaluations of F. I would therefore either like to replace the scalar element with an array of values or extend the dimension of F after performing these manipulations to the expression. I could then use subs() or xreplace() to substitute in my measurement values.

If anyone knows how to perform such an "illegal" manipulation that would be greatly appreciated. Other comments and criticisms are also more than welcome.

Edit: replaced item.expr with W and added sympy import to be explicit

[moorepants]

Your code fails for me:

In [1]: import sympy as sm

In [2]: %paste
a1, b1 = sm.symbols('a1 b1')
F = sm.MatrixSymbol('F', 3, 3)
C = F.T*F
I1 = sm.Matrix(C).trace()
W = a1/(2*b1)*(I1**(-1/3))
P = sm.derive_by_array(item.expr, sm.Matrix(F))

## -- End pasted text --
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[2], line 6
      4 I1 = sm.Matrix(C).trace()
      5 W = a1/(2*b1)*(I1**(-1/3))
----> 6 P = sm.derive_by_array(item.expr, sm.Matrix(F))

NameError: name 'item' is not defined

[RPKrijnen]

You are entirely correct. That was some bad copy pasting from my original code base. Replacing item.expr with W should fix it.
I also updated the code above:

import sympy as sm

a1, b1 = sm.symbols('a1 b1')
F = sm.MatrixSymbol('F', 3, 3)
C = F.T*F
I1 = sm.Matrix(C).trace()
W = a1/(2*b1)*(I1**(-1/3))
P = sm.derive_by_array(W, sm.Matrix(F))

[moorepants]

Does this do what you want:

In [7]: f = sm.lambdify([a1, b1] + list(F), (W, P))

In [8]: import numpy as np

In [11]: f(*np.ones(11))
Out[11]: 
(0.24037492838456825,
 array([[-0.01780555, -0.01780555, -0.01780555],
        [-0.01780555, -0.01780555, -0.01780555],
        [-0.01780555, -0.01780555, -0.01780555]]))

In [13]: f(*np.ones((11, 3)))
Out[13]: 
(array([0.24037493, 0.24037493, 0.24037493]),
 array([[[-0.01780555, -0.01780555, -0.01780555],
         [-0.01780555, -0.01780555, -0.01780555],
         [-0.01780555, -0.01780555, -0.01780555]],
 
        [[-0.01780555, -0.01780555, -0.01780555],
         [-0.01780555, -0.01780555, -0.01780555],
         [-0.01780555, -0.01780555, -0.01780555]],
 
        [[-0.01780555, -0.01780555, -0.01780555],
         [-0.01780555, -0.01780555, -0.01780555],
         [-0.01780555, -0.01780555, -0.01780555]]]))

[RPKrijnen]

Almost, seeing as I want to optimize over a1 and b1 I would prefer not to reevaluate all the components dependent on F. Something more along the lines of what I am looking for might be:

import sympy as sm
import numpy as np

a1, b1 = sm.symbols('a1 b1')
F = sm.MatrixSymbol('F', 3, 3)
C = F.T*F
I1 = sm.Matrix(C).trace()
W = a1/(2*b1)*(I1**(-1/3))
P = sm.derive_by_array(W, sm.Matrix(F))

# single measurement
f = sm.lambdify([a1, b1] + list(F), (W, P))
f_opt = lambda a,b : f(*np.append([a, b], np.ones(9))) 
f_opt(2.1, 4.2)

# multiple measurements
f_opt = lambda a,b : f(*np.append([np.repeat([a], 3), np.repeat([b], 3)], np.ones((9, 3)), axis=0))
f_opt(2.1, 4.2)

But as we increase the amount of measurements (np.ones) I expect to parse huge arrays to the function which might get inefficient.
My ideal scenario would be something that could precompute the function with the measurements without entering the parameter values a1 and b1 explicitly i.e. resulting in a function with only a1 and b1 as remaining variables.

But to be honest I am definitely no expert on how functions in Python work under the hood. So I might be envisioning issues here that don't exist.

[moorepants]

lambdify is using numpy under the hood. If you know that the numerical inputs to your functions are too large for numpy (i.e. working in memory), then you'll need to use other ways to evaluate the functions. But if your data fits in memory and works fine with NumPy, I don't see any issue with what you have written in your #multiple measurements statement.

[moorepants]

If you use ufuncify() instead of lambdify() it will fuse all operations in the numpy loop, for speed gains. The numexpr backend for lambdify also fuses the operations in the underlying loop. If you still want it to be faster then you will need to effectively inline the fixed numerical values from your measurements, by printing your own C/C++/Fortan/etc code from the SymPy expressions. lambdify() can't really do what you desire, as it isn't designed in a way that allows inline numerical data from arrays.

[RPKrijnen]

Thank you for your help! I will close this discussion as I have no further questions on this topic.