I am in a situation where I need to work with parameterized matrices. For e.g., say I start with two matrices A and B,
A = [1 2] B = [a b]
[3 4] [5 6]
Here matrix B is parameterized with the variables a and b. I might at some point need to combine the matrices, say using matrix multiplication, to get AB = C:
C = AB = [a+10 b+12 ]
[3a+20 3b+24]
I would now like to store matrix C with those variables in there and have the ability to provide values for the variables and evaluate the elements. I.e. matrix C would be a matrix parameterized by the variables a and b.
Instead of generating C on the fly any time I need it by providing the parameters to B and doing all the matrix multiplications, I would like to store the general structure of the matrix C by doing the matrix multiplications once, with C inheriting any of the variables of its constituents. I am hoping this will save run time since the structure of matrix C is cached, i.e. I don't need to do the matrix operations every time I change a parameter.
In my case, there maybe many such matrices, each with their own variable elements and I might need to combine them arbitrarily (matrix products, Kronecker products etc.). I was wondering if there is an established way to achieve this in Python. It would also be helpful if the solution was performant as the there will be many matrix operations in the code I am trying to run. Unfortunately, I have to stick with Python for now, so solutions in other languages would be less helpful.
I have tried achieving this result using the sympy module, and I think I have a very crude setup going that seems to work, but I wanted to know if there was a more canonical way of doing this in Python. Also, I am not sure if sympy would be the most performant, and I am not sure if such a system can be obtained using numpy. A simplified version of my sympy code is given below for reference:
import sympy as s
from sympy import Matrix
from sympy.physics.quantum import TensorProduct as tp
A = Matrix([[1,2],
[3,4]])
a, b = s.symbols('a b')
B = Matrix([[a,b],
[5,6]])
# To obtain C, parameterized by a and b
C = A*B
display(C)
# Finally, to evaluate C with parameters a=3 and b=0.5
C_eval = C.evalf(subs={a:3, b:0.5})
display(C_eval)
# Or, to generate and evaluate a different combination of A and B
D = tp(A,B)
display(D)
D_eval = D.evalf(subs={a:3, b:0.5})
display(D_eval)
I could see here several solutions.
SymPy great for this! It allows you cache expressions like matrix multiplications, but it's really primarily a symbolic library, so not well optimised for performance as other numeric libraries.
NumPy is extremely fast and you could combine it with numexpr:
import numpy as np
import numexpr as ne
A = np.array([[1, 2], [3, 4]])
B_expr = 'a * B1 + b * B2'
B1 =np.array([[1, 0],[0, 5]])
B2 = np.array([[0, 1],[0, 6]])
AB_expr = 'A1 * (a * B1 + b * B2)'
A1 = A
a_val = 3
b_val = 0.5
ne.evaluate(AB_expr, local_dict={'a':a_val,'b':b_val, 'A1':A1, 'B1': B1, 'B2': B2})
SymEngine is C++ based symbolic library and its much faster than SymPy. You could use it similar as SymPy by python library
import symengine as se
a, b = se.symbols('a b')
A = se.Matrix([[1, 2], [3, 4]])
B =se.Matrix([[a, b],[5, 6]])
C = A * B
C_eval = C.subs({a: 3, b: 0.5})
So if you need couples symbolic logic and caching you could use SymPy or SymEngine
If you want more performance, you could use NumPy + numexpr
Also, for very complex operations you could read about Theano