Reverse-mode automatic differentiation for multivariable gradients

A minimal backpropagation-style engine: one forward pass builds a computation graph, then one backward pass computes the gradient with respect to every input in a single traversal – the technique underlying modern deep-learning frameworks, applied here to a small analytic function.

import math

from mathematicskit.calculus.systems.autodiff import Variable, gradient

Gradient of a simple multivariable function

f = lambda x, y: x * x * y + y.sin() * x
x0, y0 = 2.0, 0.5
grad = gradient(f, [x0, y0])
print(f"gradient at ({x0}, {y0}): {grad}")

# Closed-form check: df/dx = 2xy + sin(y), df/dy = x^2 + x*cos(y)
expected_dx = 2.0 * x0 * y0 + math.sin(y0)
expected_dy = x0**2 + x0 * math.cos(y0)
print(f"closed-form gradient:      [{expected_dx}, {expected_dy}]")

gradient at (2.0, 0.5): [2.479425538604203, 5.7551651237807455]
closed-form gradient:      [2.479425538604203, 5.7551651237807455]

Shared subexpressions accumulate gradient contributions correctly

x = Variable(3.0)
y = x * x + x  # used twice: y = x^2 + x, dy/dx = 2x + 1
y.backward()
print(f"d(x^2 + x)/dx at x=3: {x.grad} (expected {2 * 3.0 + 1})")

d(x^2 + x)/dx at x=3: 7.0 (expected 7.0)