Symbolic matrix multiplication in SymPy: derive the product kernel B(n,k) from A(n,k) using Pascal's triangle

Learn to compute symbolic matrix product kernels in SymPy: compose A(n,k) via summation, avoid element-wise pitfalls, and verify with a Pascal matrix example.

When a matrix is defined by a closed-form formula A(n, k), sometimes you need the formula for the matrix product A · A without ever materializing any concrete matrices. The goal is to get a symbolic expression for the product kernel B(n, k) directly from the definition of A(n, k). Using Pascal’s triangle as a running example, below is a compact guide to doing exactly that with SymPy.

Problem setup

The intent is to compute the formula for the product kernel, not to evaluate numeric matrices. For illustration only, we’ll also show finite truncations of the matrices to verify the result. The example uses the Pascal matrix A and its product with itself B. An element-wise product C is included to highlight what we do not want.

import sympy as sp


def expr_to_grid(expr, dim):
    grid = sp.zeros(dim, dim)
    for _r in range(dim):
        for _c in range(_r + 1):
            grid[_r, _c] = expr.subs(row, _r).subs(col, _c).doit()
    return grid


row = sp.Symbol('n', integer=True, nonnegative=True)
col = sp.Symbol('k', integer=True, nonnegative=True)
mid = sp.Symbol('i', integer=True, nonnegative=True)


# base kernel A ############################################################

expr_a = sp.binomial(row, col)

arr_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 expr_to_grid(expr_a, 5) == arr_a


# product kernel B (desired) ###############################################

expr_b = sp.Sum(
    sp.binomial(row, mid) * sp.binomial(mid, col),
    (mid, 0, row)
)

arr_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 expr_to_grid(expr_b, 5) == arr_b

assert arr_b == arr_a * arr_a  # For matrices, * is matrix multiplication.


# element-wise product C (not the target) ##################################

expr_c = expr_a * expr_a  # For formulas, * acts element-wise after evaluation.

arr_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 expr_to_grid(expr_c, 5) == arr_c

What’s really going on

There are two different multiplications at play. When you multiply SymPy Matrix objects, the operator * performs matrix multiplication, summing over the shared index. When you multiply two scalar expressions, the operator * simply multiplies them pointwise. If you later evaluate such a scalar kernel at (n, k) to build a numeric grid, you end up with an element-wise product. That mismatch is why expr_a * expr_a does not represent the matrix product kernel.

To obtain the product of formula-defined matrices symbolically, you need an explicit summation over the inner index that joins rows of the left kernel with columns of the right kernel. In other words, if A(n, k) and another kernel D(n, k) define two matrices, then their product kernel at position (n, k) is Sum(A(n, i) * D(i, k), i = 0..n) in this setup.

Solution: encode the multiplication rule as a summation

Implement the symbolic composition directly by introducing a summation over the inner index. This mirrors the matrix multiplication rule and yields the desired kernel expression without materializing matrices.

def compose_kernels(left_expr, right_expr):
    return sp.Sum(
        left_expr.subs(col, mid) * right_expr.subs(row, mid),
        (mid, 0, row)
    )

expr_b_auto = compose_kernels(expr_a, expr_a)

assert expr_b_auto == expr_b

The expression expr_b_auto matches the expected product kernel. Evaluating it to a finite grid reproduces the same matrix you get from multiplying the corresponding Pascal truncations.

Why this matters

When working with operators defined by kernels A(n, k), it is common to reason purely in the symbolic domain. Having the product kernel as a closed-form expression keeps transformations algebraic, avoids premature evaluation to numeric matrices, and makes it straightforward to verify identities or pass the kernel further down a symbolic pipeline.

Takeaways

If a matrix is defined by a kernel A(n, k), its product with another kernelized matrix is obtained by inserting the intermediate index explicitly as a summation. For SymPy expressions that represent such kernels, building a Sum over the shared index produces the correct symbolic product, while naive multiplication of the expressions produces an element-wise result that does not correspond to matrix multiplication.