How to Build a Second-Order Partial Derivatives Matrix in SymPy and Collapse Trivial Zeros

Learn why is_constant() marks Derivative(x, y) as constant in SymPy and how to build a second-order partials matrix with zero-collapsing via .expr checks.

Building a matrix of second-order partials in SymPy looks straightforward until everything in the result suddenly reports as constant. If you tried to prune trivially zero terms by calling is_constant() and saw a wall of True, you’ve run into how SymPy models derivatives of independent symbols. Let’s unpack what happens and how to construct a usable “2nd order Jacobian” with clean zero-collapsing.

Problem setup

The goal is to assemble a matrix of first-order partials for a set of symbols and then stack their derivatives again with respect to each symbol, producing a second-order structure. Here is a minimal version using ASCII-only names and the same program logic.

from sympy import symbols, Matrix, zeros, Derivative
sym_names = ["Delta", "Gamma", "v"]
def make_vars():
    return symbols(" ".join(sym_names))
def first_order_partials():
    grid = []
    vars_seq = make_vars()
    for i in range(len(sym_names)):
        left = vars_seq[i]
        row = []
        for j in range(len(sym_names)):
            right = vars_seq[j]
            if left == right:
                row.append(1)
            else:
                row.append(Derivative(left, right))
        grid.append(row)
    return vars_seq, Matrix(grid)
def second_order_partials():
    vars_seq, jmat = first_order_partials()
    out = zeros(len(vars_seq)**2, len(vars_seq))
    nrows, ncols = jmat.shape
    block_idx = 0
    for k in range(len(vars_seq)):
        denom = vars_seq[k]
        for r in range(nrows):
            for c in range(ncols):
                num = jmat[r, c]
                d2 = Derivative(num, denom)
                out[block_idx * len(vars_seq) + r, c] = d2
        block_idx += 1
    return out

After building this, trying to filter by is_constant() at the element level leads to everything appearing constant.

Why everything looks constant

The key point is that a partial derivative of one independent symbol with respect to another is zero. That is exactly what SymPy encodes:

from sympy import symbols, Derivative, diff
x, y = symbols('x y')
d = Derivative(x, y)
# Shows an unevaluated derivative
print(d)
# SymPy knows this is constant (it is equal to 0)
print(d.is_constant())
# Evaluating the derivative produces 0
print(d.doit())
print(diff(x, y))

So Derivative(x, y).is_constant() returns True because the derivative equals 0. The same logic applies to Derivative(1, x), which also reduces to 0. If your first-order matrix contains either constants or derivatives of independent symbols, every second derivative of those entries is still constant, hence the uniform True.

Getting non-constant derivatives

If you expect nonzero derivatives, the thing being differentiated must actually depend on the variable. In SymPy that means using a function, not a bare symbol:

from sympy import symbols, Function, Derivative
x = symbols('x')
f = Function('f')
u = Derivative(f(x), x)
print(u)          # stays unevaluated
print(u.doit())   # still a symbolic derivative of f(x)
print(u.is_constant())  # False

This is the only way to represent an expression that truly varies with x rather than an independent symbol standing by itself. There is no way to make a set of symbols all be functions of each other.

Collapsing trivially zero terms correctly

If your intent is to keep derivatives unevaluated but eliminate those that are guaranteed to be zero, you should not test the Derivative object itself with is_constant(). Instead, check whether the expression being differentiated is constant and act accordingly. This pattern removes things like Derivative(1, x) while leaving Derivative(x, y) in its unevaluated form if you wish to keep it symbolic. When you do want to collapse to zero, you can set those entries to 0 where the derivative’s .expr is constant.

from sympy import symbols, Matrix, zeros, Derivative
sym_names = ["Delta", "Gamma", "v"]
def make_vars():
    return symbols(" ".join(sym_names))
def first_order_partials():
    grid = []
    vars_seq = make_vars()
    for i in range(len(sym_names)):
        left = vars_seq[i]
        row = []
        for j in range(len(sym_names)):
            right = vars_seq[j]
            if left == right:
                row.append(1)
            else:
                row.append(Derivative(left, right))
        grid.append(row)
    return vars_seq, Matrix(grid)
def second_order_partials_collapsed():
    vars_seq, jmat = first_order_partials()
    out = zeros(len(vars_seq)**2, len(vars_seq))
    nrows, ncols = jmat.shape
    block_idx = 0
    for k in range(len(vars_seq)):
        denom = vars_seq[k]
        for r in range(nrows):
            for c in range(ncols):
                d2 = Derivative(jmat[r, c], denom)
                # Remove trivially zero derivatives (expression is constant)
                out_idx = block_idx * len(vars_seq) + r
                out[out_idx, c] = 0 if d2.expr.is_constant() else d2
        block_idx += 1
    return out

With this approach, anything like Derivative(1, Delta) will become 0 immediately. Entries whose differentiated expression is not constant are preserved as symbolic derivatives.

Why this matters

Understanding what SymPy considers “constant” prevents misleading diagnostics. is_constant() on a Derivative reports whether the derivative’s value is independent of variables, not whether the inner expression is variable. That’s why Derivative(x, y).is_constant() is True: the derivative is zero. If you need to decide whether a derivative is trivially zero before evaluation, you must inspect the expression being differentiated via di.expr.is_constant().

Takeaways

When assembling higher-order derivative matrices in SymPy, be clear about the distinction between independent symbols and functions of symbols. If you expect variability, make it explicit with Function. If your goal is to keep symbolic structure but drop guaranteed zeros, check the derivative’s .expr rather than the derivative object itself, and replace those entries with 0. This small shift aligns the code with SymPy’s semantics and yields a cleaner, more informative second-order matrix.