I have formulas defining matrices.
The result I want is the formula defining their matrix product.
At no point do I want actual matrices. Those shown below are just for illustration.
The examples I use are Pascal's triangle A and its matrix product with itself B.
(The element-wise product C is also shown, but not what I want.)
I want to calculate sym_b
from sym_a
:
sym_b = some_magic(sym_a, sym_a)
In the code below I have written sym_b
myself.
import sympy as sp
def formula_to_matrix(formula, size):
matrix = sp.zeros(size, size)
for _n in range(size):
for _k in range(_n + 1):
matrix[_n, _k] = formula.subs(n, _n).subs(k, _k).doit()
return matrix
n = sp.Symbol('n', integer=True, nonnegative=True)
k = sp.Symbol('k', integer=True, nonnegative=True)
i = sp.Symbol('i', integer=True, nonnegative=True)
# initial matrix A #########################################################
sym_a = sp.binomial(n, k)
mat_a = sp.Matrix([
[1, 0, 0, 0, 0],
[1, 1, 0, 0, 0],
[1, 2, 1, 0, 0],
[1, 3, 3, 1, 0],
[1, 4, 6, 4, 1]
])
assert formula_to_matrix(sym_a, 5) == mat_a
# matrix product B (correct) ###############################################
sym_b = sp.Sum(
sp.binomial(n, i) * sp.binomial(i, k),
(i, 0, n)
)
mat_b = sp.Matrix([
[ 1, 0, 0, 0, 0],
[ 2, 1, 0, 0, 0],
[ 4, 4, 1, 0, 0],
[ 8, 12, 6, 1, 0],
[16, 32, 24, 8, 1]
])
assert formula_to_matrix(sym_b, 5) == mat_b
assert mat_b == mat_a * mat_a # For matrices the * is matrix multiplication.
# element-wise product C (wrong) ###########################################
sym_c = sym_a * sym_a # For formulas the * is element-wise multiplication.
mat_c = sp.Matrix([
[1, 0, 0, 0, 0],
[1, 1, 0, 0, 0],
[1, 4, 1, 0, 0],
[1, 9, 9, 1, 0],
[1, 16, 36, 16, 1]
])
assert formula_to_matrix(sym_c, 5) == mat_c
I was hoping, that SymPy provides a function for this purpose.
But apparently I have to write it myself, using the rules of matrix multiplication.
def product_of_matrix_formulas(formula_a, formula_b):
return sp.Sum(
formula_a.subs(n, n).subs(k, i) * formula_b.subs(n, i).subs(k, k),
(i, 0, n)
)
sym_b_calculated = product_of_matrix_formulas(sym_a, sym_a)
assert sym_b_calculated == sym_b