Символическое произведение матриц по ядрам в SymPy

Показываем, как в SymPy получить формулу ядра B(n,k) для произведения матриц: Sum по внутреннему индексу. Пример на треугольнике Паскаля, без явных матриц.

Когда матрица задана замкнутой формулой A(n, k), иногда требуется получить формулу для произведения матриц A · A, не формируя явные числовые матрицы. Цель — вывести символическое выражение для ядра произведения B(n, k) напрямую из определения A(n, k). На примере треугольника Паскаля — краткое руководство по тому, как сделать это в SymPy.

Постановка задачи

Здесь мы хотим получить именно формулу для ядра произведения, а не численно вычислять матрицы. Для наглядности покажем конечные усечения матриц, чтобы проверить результат. В примере используется матрица Паскаля A и её произведение на саму себя B. Для контраста добавлен покомпонентный (поточечный) продукт C — это как раз то, что нам не нужно.

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)


# базовое ядро 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


# ядро произведения B (то, что нужно) ######################################

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  # Для матриц оператор * — это матричное умножение.


# покомпонентный продукт C (не то, что нужно) ##############################

expr_c = expr_a * expr_a  # Для формул оператор * действует покомпонентно после вычисления.

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

Что здесь происходит на самом деле

Есть два разных вида умножения. Когда вы перемножаете объекты SymPy Matrix, оператор * выполняет матричное умножение, суммируя по общему индексу. Когда же вы перемножаете два скалярных выражения, оператор * просто перемножает их поэлементно. Если затем вычислить такое скалярное ядро в точке (n, k) и собрать числовую таблицу, получится покомпонентный продукт. Из‑за этого несоответствия expr_a * expr_a не является ядром матричного произведения.

Чтобы символически получить произведение матриц, заданных формулами, нужно явным образом ввести суммирование по внутреннему индексу, которое связывает строки левого ядра со столбцами правого. Иначе говоря, если A(n, k) и другое ядро D(n, k) задают две матрицы, то их ядро произведения в позиции (n, k) в этой схеме равно Sum(A(n, i) * D(i, k), i = 0..n).

Решение: закодировать правило умножения через суммирование

Реализуйте символическую композицию напрямую, введя суммирование по внутреннему индексу. Это повторяет правило матричного умножения и даёт нужное выражение ядра без материализации матриц.

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

Выражение expr_b_auto совпадает с ожидаемым ядром произведения. Если вычислить его на конечной сетке, получится та же матрица, что и при перемножении соответствующих усечений матрицы Паскаля.

Зачем это нужно

При работе с операторами, заданными ядрами A(n, k), часто удобно оставаться в чисто символическом мире. Замкнутая формула для ядра произведения сохраняет преобразования на алгебраическом уровне, позволяет избежать преждевременного перехода к числовым матрицам и упрощает проверку тождеств или передачу ядра дальше по символическому конвейеру.

Итоги

Если матрица задана ядром A(n, k), то её произведение с другой «ядерной» матрицей получается явным введением промежуточного индекса в виде суммирования. Для выражений SymPy, представляющих такие ядра, построение Sum по общему индексу даёт корректное символическое произведение, тогда как наивное перемножение выражений приводит к покомпонентному результату, который не соответствует матричному умножению.