Task Jacobian in Applied Mathematics

• The task Jacobian is a differential map that represents the derivative of a task map with respect to its input variables, used across robotics, optimization, dynamical systems, numerical simulation, machine learning, and other fields. Its applications include mapping kinematic quantities in robotics and describing sensitivity in machine learning models.
• In dynamical systems, the Jacobian describes the local behavior of a system perturbed from a given state, which can be applied to stability analysis and subsystem control further supported by concepts such as the reachability Gramian, providing insights into control properties over time. JacobianODE, meantime, relies specifically on time-series data to estimate state Jacobians without learning a vector field.
• Methodologies for efficient computation of Jacobians, such as matrix-free differentiation and scheduled Jacobian chaining, help streamline complex computations and optimizations, whereas the Jacobian conjecture in algebraic geometry explores the conditions under which polynomials are automorphisms.

A task Jacobian is the derivative of a task map with respect to its input variables. For a differentiable mapping $\boldsymbol{y}=\boldsymbol{f}(\boldsymbol{x})$, it is

$$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$

The term is used across robotics, optimization, dynamical systems, numerical simulation, machine learning, finite-element methods, computational physics, and inverse problems. Its interpretation depends on the task variables: it may map joint velocities to end-effector velocities, parameters to multiple objective values, states to vector-field velocities, inputs to neural-network logits, or generalized coordinates and their derivatives to kinematic and dynamic quantities. Although these objects share the same differential definition, their dimensions, sparsity, conditioning, computational representation, and role in downstream algorithms can differ substantially.

1. Mathematical definition and differential interpretation

Let $\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ be differentiable. The task Jacobian is the matrix

$$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$

with entries

$$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$

Its rows are output sensitivities with respect to all input coordinates, whereas its columns describe the output variation induced by perturbing individual input coordinates. The first-order approximation

$$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$

is the local linearization underlying task-space control, Newton methods, sensitivity analysis, stability analysis, and gradient-based optimization.

The dimensions of the task Jacobian determine the local structure of the mapping. If $m<n$, the task is underdetermined relative to the input dimension and the Jacobian has a nontrivial input-space null space. If $m=n$, it is locally square. If $m>n$, the task contains more output constraints than input degrees of freedom and is generally treated in a least-squares or constrained sense. For a sequence of differentiable tasks,

$$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$

the total Jacobian is the ordered product

$$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$0

The order is significant, although the parenthesization can be changed without changing the exact product. This distinction motivates Jacobian matrix-chain optimization and scheduled Jacobian chaining (Naumann, 2024, Märtens et al., 9 May 2025).

A task Jacobian is local. Its entries describe infinitesimal or sufficiently small perturbations around a particular state, trajectory, activation pattern, or parameter configuration. A small derivative does not by itself establish global independence, global robustness, or global causal irrelevance. Conversely, an accurate prediction map does not necessarily imply an accurate derivative map: a sequence of functions can converge pointwise while their derivatives fail to converge. This issue is explicit in data-driven dynamical-system Jacobian estimation (Eisen et al., 2 Jul 2025).

The term should also be distinguished from related derivatives. An input-output Jacobian differentiates model outputs with respect to inputs; a parameter-function Jacobian differentiates outputs with respect to parameters; a loss gradient is a vector rather than the Jacobian of an entire task map; and a training Jacobian differentiates the final trained parameters with respect to the initialization parameters,

$$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$1

rather than differentiating predictions with respect to inputs (Belrose et al., 2024).

2. Task Jacobians in robotics and motion control

In robot kinematics, the task map is commonly the forward-kinematics function

$$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$2

where $$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$3 is the joint-coordinate vector and $$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$4 is an end-effector or task-space representation. Differentiation gives

$$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$5

The columns of $$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$6 are the instantaneous task-space velocities produced by unit joint velocities. In visual servoing, $$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$7 can contain image-plane feature coordinates, so the Jacobian includes the sensor and projection mapping and is often called a visual-motor or image Jacobian. It is not a dynamic Jacobian: inertial, Coriolis, gravitational, and torque relationships are separate objects (Przystupa et al., 2021).

Differential inverse kinematics commonly uses the Moore–Penrose pseudoinverse:

$$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$8

The first term addresses the task error. The second lies in the input-space null space and can represent redundant motion. Under the full-row-rank assumption,

$$\boldsymbol{J}=\frac{\partial \boldsymbol{y}}{\partial \boldsymbol{x}}, \qquad \delta\boldsymbol{y}=\boldsymbol{J}\,\delta\boldsymbol{x}.$$9

Near rank deficiency, this inverse can become unstable and produce large joint commands; the cited visual-servoing study analyzes conditioning but does not introduce a damped least-squares inverse (Przystupa et al., 2021).

For dual-arm manipulation, the task may be decomposed into absolute and relative motion. With end-effector twists $\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$0,

$\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$1

The conventional relative Jacobian is

$\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$2

where $\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$3 is the block-diagonal Jacobian of the two manipulators and $\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$4 maps end-effector twists to object-frame twists. The relative task has an absolute-motion redundancy: many pairs of end-effector velocities produce the same relative velocity. This redundancy can support collision avoidance, joint-limit avoidance, posture optimization, or reduced joint velocity.

An extended relative Jacobian introduces a user-selected asymmetry parameter $\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$5:

$\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$6

with

$\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$7

The corresponding task-space resolution is

$\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$8

Thus $\boldsymbol{f}:\mathbb{R}^{n}\to\mathbb{R}^{m}$9 gives symmetric execution, while $$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$0 assigns nearly all relative motion to arm 2 and $$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$1 assigns nearly all of it to arm 1. The proposed construction places asymmetry in the relative task rather than redefining the physical meaning of absolute motion (Almeida et al., 2019).

For higher-order motion optimization, the input may include

$$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$2

and the output may be a pose, velocity, acceleration, jerk, momentum, force, torque, or torque rate. The relevant Jacobian is then a derivative with respect to generalized coordinates and their higher-order derivatives. A conventional velocity Jacobian may reproduce the correct instantaneous velocity while failing to equal the mathematical derivative of a higher-order task map under arbitrary perturbations. Structured representations based on comprehensive motion transformation matrices preserve the multi-link tree structure and avoid repeatedly expanding higher-order chain and product rules (Ishigaki et al., 15 May 2026).

3. Objective, input-output, and training Jacobians in machine learning

For a classifier with logits

$$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$3

the input-output task Jacobian is

$$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$4

It describes how each output logit responds to each input coordinate. A first-order perturbation satisfies

$$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$5

Jacobian regularization penalizes the squared Frobenius norm,

$$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$6

through an objective of the form

$$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$7

The Frobenius norm controls aggregate sensitivity rather than the worst-case sensitivity, which is governed by the spectral norm. Since

$$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$8

reducing the Frobenius norm provides an upper-bound-type control on local Euclidean amplification but does not directly minimize the most sensitive direction. In classification, smaller derivatives of logit differences can increase the local linearized distance to a decision boundary. Empirically, Jacobian regularization reduced measured input-output Jacobians and improved robustness to random and adversarial perturbations in the reported experiments, with a robustness–accuracy trade-off (Hoffman et al., 2019).

Explicitly constructing a $$\boldsymbol{J}(\boldsymbol{x}) = \frac{\partial \boldsymbol{f}(\boldsymbol{x})}{\partial \boldsymbol{x}} \in\mathbb{R}^{m\times n},$$9 Jacobian is expensive. A stochastic output projection uses a random unit vector $$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$0 and the identity

$$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$1

The resulting estimator uses vector-Jacobian products and requires differentiating a gradient-dependent penalty with respect to model parameters, commonly called double backpropagation. The cited method reports approximately constant-factor rather than class-count-proportional overhead, although the procedure remains more expensive than ordinary training (Hoffman et al., 2019).

In multi-objective optimization, the task Jacobian instead contains one objective gradient per row. For

$$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$2

$$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$3

A Jacobian-descent update has the form

$$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$4

where the aggregator maps the complete objective Jacobian to a parameter-space direction. The proposed $$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$5 aggregator projects each task gradient into the dual cone

$$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$6

and averages the projected gradients. It is designed to be non-conflicting, linear under positive row scaling, and weighted in the task-gradient span (Quinton et al., 2024).

The associated Gramian,

$$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$7

contains pairwise task-gradient inner products. Positive off-diagonal entries indicate alignment, negative entries indicate conflict, and zero entries indicate orthogonality. When the number of objectives is much smaller than the number of parameters, Gramian-based computation can avoid storing the full $$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$8 Jacobian, although it still requires $$J_{ij}(\boldsymbol{x}) = \frac{\partial f_i(\boldsymbol{x})}{\partial x_j}.$$9 storage and can be expensive for very large objective sets (Quinton et al., 2024).

A different object is the training Jacobian. Fixing data, labels, optimizer, learning-rate schedule, minibatch sequence, and training duration defines a map

$$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$0

The training Jacobian is

$$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$1

Its singular values characterize local amplification, preservation, or attenuation of initialization perturbations. The reported experiments identify chaotic, bulk, and stable spectral regions; bulk directions have singular values near one and nearly aligned left and right singular vectors, indicating approximate transport of initialization perturbations through training. This quantity is neither the input-output Jacobian nor the Hessian (Belrose et al., 2024).

For transformer interpretability, Jacobian Sparse Autoencoders apply the same differential concept to the activation-to-activation map of an individual MLP in learned sparse feature coordinates. If $$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$2 and $$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$3 are sparse input and output SAE features, the relevant Jacobian is

$$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$4

Its entries are interpreted as weighted feature-to-feature computational edges. TopK sparsification reduces the relevant Jacobian to an active $$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$5 block, and an $$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$6 penalty on this block encourages computational sparsity in addition to activation sparsity. Efficient chain-rule evaluation avoids materializing the full latent Jacobian (Farnik et al., 25 Feb 2025).

4. Dynamical systems, subsystem control, and inverse problems

For a nonlinear continuous-time system,

$$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$7

the state Jacobian is

$$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$8

A perturbation $$\boldsymbol{f}(\boldsymbol{x}+\delta\boldsymbol{x}) \approx \boldsymbol{f}(\boldsymbol{x}) + \boldsymbol{J}(\boldsymbol{x})\delta\boldsymbol{x}$$9 evolves according to the tangent dynamics

$m<n$0

If the state is partitioned into subsystems $m<n$1 and $m<n$2, the cross-block $m<n$3 describes the instantaneous influence of perturbations in $m<n$4 on the velocity of $m<n$5. Treating $m<n$6 as a local control input gives

$m<n$7

where

$m<n$8

The reachability Gramian

$m<n$9

combines cross-subsystem coupling with intrinsic target dynamics. Its trace measures average reachability, while its minimum eigenvalue measures the hardest local control direction. Thus a large instantaneous cross-block does not by itself imply large long-horizon control: the state-transition matrix can amplify, rotate, or suppress the perturbation (Eisen et al., 2 Jul 2025).

JacobianODE directly estimates state Jacobians from time-series data rather than learning a vector field and differentiating it. Its loss combines trajectory reconstruction with a loop-closure penalty. Since each row of a true Jacobian is a gradient field, its line integral around a closed path should vanish:

$m=n$0

This integrability constraint distinguishes a dynamically adequate but geometrically incorrect derivative estimate from a valid Jacobian. The method was evaluated on the Van der Pol oscillator, Lorenz systems, and Lorenz96 systems, where Jacobian accuracy differed substantially from trajectory-prediction accuracy (Eisen et al., 2 Jul 2025).

In elliptic inverse problems, a generalized task Jacobian is formed from several scalar solutions. For solutions $m=n$1,

$m=n$2

A rank-$m=n$3 condition ensures that some $m=n$4 columns form a nonzero determinant. For piecewise-regular coefficients, raw gradients may jump across interfaces. An interface-adapted quantity replaces the normal gradient by the continuous normal flux while retaining tangential derivatives:

$m=n$5

This quantity is continuous under the stated interface conditions and allows non-vanishing determinant constraints to be propagated across discontinuous coefficients. The resulting open-and-dense families of elliptic solutions support internal-data reconstruction in multi-wave inverse problems (Capdeboscq et al., 2023).

The classical Jacobian determinant also appears in the Jacobian conjecture. For a polynomial map

$m=n$6

with Jacobian matrix

$m=n$7

the conjecture asks whether

$m=n$8

implies that $m=n$9 is a polynomial automorphism. Compactification, blow-ups, exceptional curves, pole orders, and Stein factorization provide a geometric framework for studying behavior at infinity in dimension two. The cited geometric method obtains structural constraints and conditional automorphism criteria but does not prove the full conjecture (Sun, 2014). Other proposed approaches use reductions to Drużkowski maps, optimization arguments, Hadamard’s global inverse theorem, or perturbative tree expansions; the supplied analyses emphasize that the central properness or finite-tree steps in such claimed proofs remain unestablished (Liu, 2020, Magnen, 2023).

5. Numerical Jacobian construction and automatic differentiation

In numerical simulation, task Jacobians often arise from discretized residuals. A finite-element discretization produces

$m>n$0

Newton–Raphson solves

$m>n$1

Because finite-element basis functions are local, the global Jacobian is sparse. Differentiating the element residual instead of the assembled global residual yields a small dense local Jacobian $m>n$2, which is assembled as

$m>n$3

Source-transformation automatic differentiation with TAPENADE can generate the local derivatives. Forward vector mode seeded with the identity produces the complete local Jacobian, after which the existing sparse assembly infrastructure constructs the global matrix. The approach avoids requiring the AD system to discover the mesh-dependent global sparsity pattern, although sensitivity arrays and generated code can remain expensive (Zwicke et al., 2017).

Chemical-kinetics task Jacobians have a different structural issue. For a constant-pressure thermochemical ODE with state

$m>n$4

pyJac constructs the full Jacobian

$m>n$5

This includes temperature-energy coupling, density derivatives, heat-capacity derivatives, species production derivatives, and the dependent-species constraint. Under constant pressure, $m>n$6, making the mass-fraction formulation generally dense even when individual reactions involve few species. pyJac uses analytical reaction derivatives, common-subexpression reuse, reaction-level accumulation, and CPU/GPU code generation rather than finite-difference source evaluations (Niemeyer et al., 2016).

Entropy-stable summation-by-parts schemes exploit a different form of structure. If the residual is

$m>n$7

where $m>n$8 contains symmetric two-point fluxes, then for skew-symmetric $m>n$9,

$$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$0

Only derivatives of the small two-point flux are required. Automatic differentiation can be applied to that local flux function rather than to the full spatial residual, preserving SBP sparsity and reducing the cost relative to direct AD or finite differences (Chan et al., 2020).

Power-flow Jacobians exploit sparse network structure in a related way. With complex bus voltage $$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$1 and admittance matrix $$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$2,

$$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$3

The Newton Jacobian is partitioned into derivatives of active and reactive power with respect to voltage angles and magnitudes:

$$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$4

The derivative matrices share the admittance matrix’s CRS sparsity pattern. Direct row-wise loops can therefore construct the final CRS Jacobian without repeatedly creating sparse diagonal matrices, multiplying sparse matrices, or rebuilding row and column metadata (Schäfer et al., 2018).

For a modular program, the complete task Jacobian is a chain of elemental Jacobians that may not be explicitly available. Tangent routines compute

$$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$5

whereas adjoint routines compute

$$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$6

Matrix-free Jacobian chaining chooses between tangent propagation, adjoint propagation, explicit intermediate matrices, and dense products. Dynamic programming minimizes arithmetic cost while incorporating graph sizes and tape memory. The resulting recurrence differs from classical matrix-chain multiplication because the factors are implicit and available only through tangent and adjoint operations (Naumann, 2024).

Scheduled Jacobian chaining extends this problem to parallel execution. A parenthesization induces precedence constraints among elemental accumulations, tangent eliminations, adjoint eliminations, and explicit multiplications. The resulting scheduling problem is modeled as an in-tree precedence problem on identical machines. A dynamic-programming heuristic optimizes both bracketing and machine allocation, while branch-and-bound supplies a global optimum for small instances (Märtens et al., 9 May 2025).

6. Computational properties, limitations, and terminology

The practical value of a task Jacobian depends on four interacting properties: accuracy, sparsity, conditioning, and the required mode of access.

Accuracy concerns whether the derivative corresponds to the intended mathematical task. Analytical differentiation, structured formulas, and automatic differentiation remove finite-difference step-size selection and subtraction cancellation, but automatic differentiation differentiates the implemented program. Conditional branches, numerical approximations, matrix-function algorithms, stabilization terms, and singular formulas may produce derivatives that are mathematically inappropriate at limiting points. The finite-element TAPENADE study therefore treats AD as requiring inspection of branches, special functions, initialization, and source-code compatibility (Zwicke et al., 2017).

Sparsity may arise from locality, network topology, reaction structure, SBP operators, active feature sets, or a chosen coordinate representation. It is not invariant under arbitrary changes of basis. A constant-pressure mass-fraction Jacobian can be dense because of density coupling, while a concentration-based formulation can expose sparsity. In finite elements, global sparsity follows from basis-function support even when element Jacobians are dense. In JSAEs, computational sparsity is explicitly learned in a selected SAE basis rather than being a basis-independent property (Niemeyer et al., 2016, Farnik et al., 25 Feb 2025).

Conditioning determines how reliably a Jacobian can be inverted or used in a solver. In robotics, poorly conditioned Jacobians can generate large joint commands. In power flow, the Jacobian enters the Newton linear solve. In multi-objective optimization, the Gramian identifies aligned, orthogonal, and conflicting task gradients. In dynamical systems, reachability-Gramian eigenvalues characterize the ease of controlling different directions. A learned Jacobian can be biased yet better conditioned than the physical Jacobian and consequently yield better closed-loop behavior (Przystupa et al., 2021, Eisen et al., 2 Jul 2025).

Access mode distinguishes explicit Jacobians from Jacobian actions. A full matrix may be necessary for eigendecomposition, preconditioning, diagnostics, or matrix assembly. Many algorithms need only a Jacobian-vector product $$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$7 or vector-Jacobian product $$\boldsymbol{F} = \boldsymbol{F}_q\circ\cdots\circ\boldsymbol{F}_1,$$8. Forward-mode AD is generally favorable for a small number of input directions; reverse-mode AD is generally favorable for a small number of output directions, at the cost of tape or intermediate-state storage. Matrix-free chaining formalizes this trade-off, while Gramian propagation replaces the full parameter Jacobian by pairwise task-gradient information when the objective dimension is manageable (Naumann, 2024, Quinton et al., 2024).

Several common misconceptions follow from these distinctions. A task Jacobian is not necessarily a Hessian, a determinant, a gradient, a causal graph, or a global sensitivity certificate. A nonzero cross-subsystem block indicates instantaneous local influence, not necessarily long-term control or biological causality. A small input-output Jacobian controls local sensitivity but does not establish global robustness. A sparse Jacobian in learned coordinates suggests a sparse local computational description but does not by itself prove a complete causal circuit. Finally, an analytic inverse supplied by a global inverse theorem is not automatically polynomial, and a proposed Jacobian-conjecture proof must separately establish the algebraic conclusion.

Across these domains, the unifying role of the task Jacobian is the representation of first-order task sensitivity. Its form is always a derivative, but its scientific meaning is determined by the task variables, the operating point, the representation, the numerical method used to construct it, and the downstream operation—control, optimization, inference, stability analysis, or interpretation—for which it is employed.