1 core.autodiff

core.autodiff

Automatic differentiation for symbolic expressions.

Implements symbolic differentiation using the chain rule, producing gradient expressions that can be compiled for fast evaluation.

Supports native gradient rules for vector expressions (VectorSum, LinearCombination, DotProduct) with O(1) coefficient lookup for scalability to n=10,000+ variables.

1.1 Classes

Name	Description
SparsityPattern	Describes which gradient elements are structurally non-zero.
VectorExpressionPattern	Classification of patterns in gradient expressions for optimization.
VectorGradientPattern	Represents a gradient of the form ∇f(x) = Ax + b.

1.1.1 SparsityPattern

core.autodiff.SparsityPattern(nnz_indices, size, is_constant, constant_values)

Describes which gradient elements are structurally non-zero.

Enables O(nnz) gradient/Jacobian computation by identifying which variables an expression depends on, without computing values.

1.1.1.1 Attributes

Name	Type	Description
nnz_indices	NDArray[np.intp]	Sorted array of indices (into the variables list) that have non-zero partial derivatives.
size	int	Total number of variables (length of full gradient vector).
is_constant	bool	True if all non-zero gradients are constants (e.g. linear expr).
constant_values	NDArray[np.floating] \| None	If is_constant, the non-zero gradient values at nnz_indices. None if gradients are not all constant.

1.1.2 VectorExpressionPattern

core.autodiff.VectorExpressionPattern()

Classification of patterns in gradient expressions for optimization.

These patterns allow the compiler to generate specialized, high-performance kernels for common vector operations instead of using generic element-wise logic.

1.1.3 VectorGradientPattern

core.autodiff.VectorGradientPattern(
    linear_term,
    constant_term,
    vector,
    linear_type=None,
    linear_scale=0.0,
    linear_diag=None,
)

Represents a gradient of the form ∇f(x) = Ax + b.

Used for vectorized compilation of gradients for: - Linear combinations (A=None, b=c) - Quadratic forms (A=Q+Q’, b=0) - Dot products (A=2I, b=0)

1.2 Functions

Name	Description
analyze_gradient_sparsity	Analyze which gradient elements are structurally non-zero.
analyze_jacobian_sparsity	Analyze the sparsity structure of the full Jacobian matrix.
apply_gradient_rule	Apply the registered gradient rule for an expression type.
compile_hessian	Compile the Hessian for fast evaluation.
compile_jacobian	Compile the Jacobian for fast evaluation.
compile_sparse_jacobian	Compile the Jacobian for fast evaluation, returning sparse output.
compute_hessian	Compute the Hessian matrix of an expression.
compute_jacobian	Compute the Jacobian matrix of expressions with respect to variables.
detect_affine_gradient_pattern	Detect if an expression has a vectorizable gradient pattern: ∇f(x) = Ax + b.
detect_vector_gradient_pattern	Detect the structural pattern of an expression.
gradient	Compute the symbolic gradient of an expression with respect to a variable.
has_gradient_rule	Check if an expression type has a registered gradient rule.
increased_recursion_limit	Temporarily increase Python’s recursion limit.
register_gradient	Decorator to register a gradient rule for an expression type.

1.2.1 analyze_gradient_sparsity

core.autodiff.analyze_gradient_sparsity(expr, variables)

Analyze which gradient elements are structurally non-zero.

Determines which variables the expression depends on by walking the expression tree, then optionally checks if those gradients are constant (for linear expressions).

1.2.1.1 Parameters

Name	Type	Description	Default
expr	Expression	The expression to analyze.	required
variables	list[Variable]	List of variables defining the gradient vector ordering.	required

1.2.1.2 Returns

Name	Type	Description
	SparsityPattern	SparsityPattern describing which gradient elements are non-zero.

1.2.2 analyze_jacobian_sparsity

core.autodiff.analyze_jacobian_sparsity(exprs, variables)

Analyze the sparsity structure of the full Jacobian matrix.

Returns per-row sparsity patterns describing which columns of each Jacobian row are structurally non-zero.

1.2.2.1 Parameters

Name	Type	Description	Default
exprs	list[Expression]	List of m expressions (rows of the Jacobian).	required
variables	list[Variable]	List of n variables (columns of the Jacobian).	required

1.2.2.2 Returns

Name	Type	Description
	list[SparsityPattern]	List of m SparsityPattern objects, one per row.

1.2.3 apply_gradient_rule

core.autodiff.apply_gradient_rule(expr, wrt)

Apply the registered gradient rule for an expression type.

1.2.3.1 Parameters

Name	Type	Description	Default
expr	'Expression'	The expression to differentiate.	required
wrt	'Variable'	The variable to differentiate with respect to.	required

1.2.3.2 Returns

Name	Type	Description
	'Expression'	The gradient expression.

1.2.3.3 Raises

Name	Type	Description
	ValueError	If no gradient rule is registered for this expression type.

1.2.4 compile_hessian

core.autodiff.compile_hessian(expr, variables)

Compile the Hessian for fast evaluation.

1.2.4.1 Parameters

Name	Type	Description	Default
expr	Expression	The expression to differentiate.	required
variables	list[Variable]	List of variables.	required

1.2.4.2 Returns

Name	Type	Description
	Callable[[NDArray[np.floating]], NDArray[np.floating]]	A callable that takes a 1D array and returns the Hessian as a 2D array.

1.2.4.3 Performance

For VectorPowerSum and VectorUnarySum, the Hessian is diagonal, so we use O(n) vectorized computation instead of O(n²).

1.2.5 compile_jacobian

core.autodiff.compile_jacobian(exprs, variables)

Compile the Jacobian for fast evaluation.

1.2.5.1 Parameters

Name	Type	Description	Default
exprs	list[Expression]	List of expressions.	required
variables	list[Variable]	List of variables.	required

1.2.5.2 Returns

Name	Type	Description
	Callable[[NDArray[np.floating]], NDArray[np.floating]]	A callable that takes a 1D array and returns the Jacobian as a 2D array.

1.2.5.3 Performance

Checks for vector patterns (O(1)) per row.
Checks for constant rows (O(1)).
Batches execution per row (m calls instead of m*n).

1.2.6 compile_sparse_jacobian

core.autodiff.compile_sparse_jacobian(exprs, variables, density_threshold=0.5)

Compile the Jacobian for fast evaluation, returning sparse output.

For sparse constraint systems (nnz << mn), returns scipy.sparse.csr_matrix with O(nnz) memory instead of O(mn). Falls back to dense for high-density Jacobians.

1.2.6.1 Parameters

Name	Type	Description	Default
exprs	list[Expression]	List of expressions (rows of the Jacobian).	required
variables	list[Variable]	List of variables (columns of the Jacobian).	required
density_threshold	float	If overall Jacobian density exceeds this, fall back to the dense compile_jacobian (default 0.5).	`0.5`

1.2.6.2 Returns

Name	Type	Description
	Callable[[NDArray[np.floating]], Any]	A callable that takes a 1D array and returns the Jacobian.
	Callable[[NDArray[np.floating]], Any]	Output is scipy.sparse.csr_matrix for sparse problems, or a dense
	Callable[[NDArray[np.floating]], Any]	2D ndarray for high-density problems.

1.2.7 compute_hessian

core.autodiff.compute_hessian(expr, variables)

Compute the Hessian matrix of an expression.

1.2.7.1 Parameters

Name	Type	Description	Default
expr	Expression	The expression to differentiate twice.	required
variables	list[Variable]	List of variables.	required

1.2.7.2 Returns

Name	Type	Description
	list[list[Expression]]	Hessian matrix as H[i][j] = d²(expr)/d(var_i)d(var_j).

1.2.7.3 Note

The Hessian is symmetric, so H[i][j] = H[j][i]. We compute the full matrix but could optimize by exploiting symmetry.

1.2.8 compute_jacobian

core.autodiff.compute_jacobian(exprs, variables)

Compute the Jacobian matrix of expressions with respect to variables.

Uses O(1) vectorized jacobian_row() methods when available for vector expressions (VectorSum, DotProduct, LinearCombination, QuadraticForm). Falls back to individual gradient() calls otherwise.

1.2.8.1 Parameters

Name	Type	Description	Default
exprs	list[Expression]	List of expressions (constraints or objectives).	required
variables	list[Variable]	List of variables to differentiate with respect to.	required

1.2.8.2 Returns

Name	Type	Description
	list[list[Expression]]	Jacobian matrix as J[i][j] = d(expr_i)/d(var_j).

1.2.8.3 Example

x, y = Variable(“x”), Variable(“y”) exprs = [x**2 + y, x*y] J = compute_jacobian(exprs, [x, y]) # J[0][0] = 2*x, J[0][1] = 1 # J[1][0] = y, J[1][1] = x

1.2.9 detect_affine_gradient_pattern

core.autodiff.detect_affine_gradient_pattern(expr)

Detect if an expression has a vectorizable gradient pattern: ∇f(x) = Ax + b.

This enables O(1) compilation of gradients for common patterns like quadratic forms, dot products, and linear combinations, entirely bypassing the potentially large expression tree.

1.2.9.1 Parameters

Name	Type	Description	Default
expr	Expression	The expression to analyze.	required

1.2.9.2 Returns

Name	Type	Description
	VectorGradientPattern \| None	VectorGradientPattern if detected, None otherwise.

1.2.10 detect_vector_gradient_pattern

core.autodiff.detect_vector_gradient_pattern(expr, wrt=None)

Detect the structural pattern of an expression.

Identifies standard forms like sums, dot products, and norms that can be optimized by the compiler.

1.2.11 gradient

core.autodiff.gradient(expr, wrt)

Compute the symbolic gradient of an expression with respect to a variable.

Uses a three-tier approach for optimal performance: 1. Registered gradient rules (O(1) for vector expressions) 2. Cached recursive computation (for shallow trees) 3. Iterative fallback (for deep trees to avoid RecursionError)

1.2.11.1 Parameters

Name	Type	Description	Default
expr	Expression	The expression to differentiate.	required
wrt	Variable	The variable to differentiate with respect to.	required

1.2.11.2 Returns

Name	Type	Description
	Expression	A new Expression representing the derivative.

1.2.11.3 Example

x = Variable(“x”) expr = x**2 + 3x grad = gradient(expr, x) # Returns: 2x + 3

1.2.12 has_gradient_rule

core.autodiff.has_gradient_rule(expr)

Check if an expression type has a registered gradient rule.

1.2.12.1 Parameters

Name	Type	Description	Default
expr	'Expression'	The expression to check.	required

1.2.12.2 Returns

Name	Type	Description
	bool	True if a gradient rule is registered for this expression type.

1.2.13 increased_recursion_limit

core.autodiff.increased_recursion_limit(limit=5000)

Temporarily increase Python’s recursion limit.

This context manager can be used as a workaround for deep expression trees when the automatic iterative/recursive switching isn’t sufficient.

.. warning:: Use with caution - very high limits can cause stack overflow crashes. The iterative gradient implementation is preferred for deep trees.

1.2.13.1 Parameters

Name	Type	Description	Default
limit	int	The temporary recursion limit (default: 5000).	`5000`

1.2.13.2 Returns

Name	Type	Description
	Iterator[None]	A context manager that restores the previous recursion limit when its
	Iterator[None]	`with` block exits.

1.2.13.3 Example

with increased_recursion_limit(5000): … grad = gradient(deep_expr, x)

1.2.14 register_gradient

core.autodiff.register_gradient(expr_type)

Decorator to register a gradient rule for an expression type.

Registered gradient rules are used by the main gradient() function before falling back to recursive tree traversal. This enables O(1) gradient computation for vector expressions.

1.2.14.1 Parameters

Name	Type	Description	Default
expr_type	type	The expression class to register a gradient rule for.	required

1.2.14.2 Returns

Name	Type	Description
	Callable[[GradientFunc], GradientFunc]	A decorator that registers the gradient function.

1.2.14.3 Example

@register_gradient(VectorSum) def gradient_vector_sum(expr: VectorSum, wrt: Variable) -> Expression: # O(1) gradient computation …