Score-Jacobian Information Matrix (SJIM)

• SJIM is a family of score-based second-order matrices that transform score geometry into tractable information operators across sample, latent, or parameter spaces.
• It encompasses forms like Fisher outer products, expected negative Hessians, and pullback matrices, enabling efficient computations in denoising, causal discovery, and inference.
• Practical implementations in SC-VAMP, causal discovery, and derivative-free estimation highlight SJIM’s value despite challenges like score calibration and nonlinear expectation mismatches.

Searching arXiv for the cited papers and topic context. The Score-Jacobian Information Matrix (SJIM) is not a single universally fixed object across the literature. Rather, it denotes a family of score-centered second-order operators that connect gradient fields, Jacobians, Hessians, and information measures. In "Score-Based VAMP with Fisher-Information-Based Onsager Correction" it is defined explicitly as a conditional Fisher information matrix built from score outer products; in "Optimization-Free Topological Sort for Causal Discovery via the Schur Complement of Score Jacobians" it is the expected negative Hessian of the log-density; in "Information Gradient for Directed Acyclic Graphs" and "Score Jacobian Chaining" the term is not explicitly introduced but a natural SJIM arises as a pullback of score power through system or renderer Jacobians; and in likelihood-based inference for state-space models and REML it coincides with the Jacobian of the score, equivalently the Hessian of the log-likelihood up to sign conventions (Wadayama et al., 11 Jan 2026, Wu et al., 28 Apr 2026, Wadayama, 5 Jan 2026, Wang et al., 2022, Doucet et al., 2013, Nemeth et al., 2013, Zhu, 2016).

1. Definitional landscape

The common thread is that SJIM records how score information propagates through a model’s differential structure. What changes across papers is the ambient space: sample space, latent-variable space, parameter space, or rendered scene space. This terminological plurality is explicit in the corpus: some papers define the matrix directly, while others only induce it through their gradient identities or Hessian constructions.

Context	SJIM form	Primary role
SC-VAMP	$$\mathbb{E}[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top]$$	Jacobian-free Onsager correction
Causal discovery	$\mathbb{E}[-\nabla_x^2 \log p(x)]$	Leaf detection and Schur elimination
DAG information gradients	$$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$	Parameter-space curvature from mutual-information scores
Score Jacobian Chaining	$\mathbb{E}[J^\top F_I J]$ or $\mathbb{E}[s_V s_V^\top]$	3D information induced by 2D diffusion scores
Likelihood-based inference	$$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$	Observed information or its negative

In SC-VAMP, the explicit definition is

$$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$

with scalar trace

$$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$

In the causal-discovery formulation, the expected SJIM is instead

$I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$

which the authors identify, under homoscedastic Gaussian additive noise, with the Fisher Information Matrix with respect to location parameters and therefore call the structural information matrix (Wadayama et al., 11 Jan 2026, Wu et al., 28 Apr 2026).

This difference in definition is substantive rather than notational. A common misconception is that SJIM always means a Fisher outer-product matrix. The cited papers show that it may instead mean an expected negative Hessian, a pullback of output score power into parameter space, or the Jacobian of a likelihood score. The unifying feature is the coupling of score geometry with a Jacobian or Hessian operator, not a single canonical formula.

2. Score geometry, Hessians, and Stein-type identities

The score function is the basic primitive:

$s(x)=\nabla_x \log p(x).$

From this starting point, the cited works develop several equivalent or near-equivalent second-order constructions. In SC-VAMP, the score outer product yields conditional Fisher information. In causal discovery, the negative log-density Hessian

$\mathbb{E}[-\nabla_x^2 \log p(x)]$0

is averaged to obtain the structural information matrix. In likelihood-based inference, if $\mathbb{E}[-\nabla_x^2 \log p(x)]$1 and $\mathbb{E}[-\nabla_x^2 \log p(x)]$2, then the score Jacobian is

$\mathbb{E}[-\nabla_x^2 \log p(x)]$3

while the observed information is

$\mathbb{E}[-\nabla_x^2 \log p(x)]$4

(Wadayama et al., 11 Jan 2026, Wu et al., 28 Apr 2026, Doucet et al., 2013).

Stein-type identities furnish the bridge between score norms, divergences, and information traces. SC-VAMP uses

$\mathbb{E}[-\nabla_x^2 \log p(x)]$5

so that for a Tweedie-form estimator

$\mathbb{E}[-\nabla_x^2 \log p(x)]$6

the average divergence becomes

$\mathbb{E}[-\nabla_x^2 \log p(x)]$7

This is the key step by which a Jacobian trace is replaced by the trace of a score-built information matrix, without explicit differentiation of the learned score model (Wadayama et al., 11 Jan 2026).

The causal-discovery paper takes the complementary viewpoint that the expected curvature of the score’s Jacobian encodes topological structure. Since $\mathbb{E}[-\nabla_x^2 \log p(x)]$8 and $\mathbb{E}[-\nabla_x^2 \log p(x)]$9, the authors describe SJIM as capturing an average second-order geometry of $$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$0 across the manifold. This geometry becomes algorithmically useful because Schur complements of block-partitioned SJIMs implement variable elimination under the paper’s additive-noise assumptions (Wu et al., 28 Apr 2026).

These formulations are compatible in Gaussian cases. In linear Gaussian ANMs, the causal-discovery paper states that $$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$1. In SC-VAMP, the Fisher-based trace identity reproduces the divergence term needed by AMP/VAMP. In score-Jacobian chaining, the induced scene-space matrix

$$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$2

is explicitly described as a Fisher-information analogue induced by the renderer and the diffusion score (Wu et al., 28 Apr 2026, Wadayama et al., 11 Jan 2026, Wang et al., 2022).

3. SC-VAMP: SJIM as Fisher-information-based Onsager correction

SC-VAMP introduces the most explicit modern definition of SJIM and makes it operational. The framework uses two SISO modules. The prior module uses the unconditional score

$$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$3

where $$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$4 and $$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$5. The likelihood module uses the conditional score

$$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$6

Both modules construct MMSE estimators through Tweedie’s formula:

$$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$7

and, in the denoiser case,

$$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$8

The same learned score object then determines the Onsager correction through the trace of the SJIM (Wadayama et al., 11 Jan 2026).

The central definition is

$$\mathbb{E}[(D_{\eta_k}Y)^\top \Delta s\,\Delta s^\top (D_{\eta_\ell}Y)]$$9

Its trace is the score-norm statistic

$\mathbb{E}[J^\top F_I J]$0

The Onsager coefficient becomes

$\mathbb{E}[J^\top F_I J]$1

which exactly matches the Jacobian/divergence-based Onsager term in expectation, but avoids explicit Jacobians. The resulting extrinsic recursion is

$\mathbb{E}[J^\top F_I J]$2

and

$\mathbb{E}[J^\top F_I J]$3

The posterior variance per coordinate is

$\mathbb{E}[J^\top F_I J]$4

Algorithmically, the trace of the SJIM is estimated by mini-batch averages of squared score norms rather than by automatic differentiation of Jacobian traces. The paper states that SC-VAMP avoids explicit Jacobians/traces and that its Onsager computation is $\mathbb{E}[J^\top F_I J]$5 per sample, batched, whereas classical Jacobian/divergence computation via AutoDiff is expensive. Learned scores can be obtained by denoising score matching or via reverse Tweedie from a pre-trained denoiser; practical robustness is improved by mini-batch Stein calibration enforcing

$\mathbb{E}[J^\top F_I J]$6

When random orthogonal or unitary mixing is added, the method is intended to restore permutation symmetry, suppress short-range correlations, and place structured channels in the RRI universality class, thereby supporting scalar state evolution even with nonlinear modules (Wadayama et al., 11 Jan 2026).

The information-theoretic interpretation is equally central. The paper states that the entropic CLT acts as a Gaussianizer for mixed errors, and that De Bruijn and I-MMSE identities connect Fisher information, MMSE, and mutual information. In the scalar Gaussian channel $\mathbb{E}[J^\top F_I J]$7 with $\mathbb{E}[J^\top F_I J]$8 and $\mathbb{E}[J^\top F_I J]$9, the SC-VAMP state-evolution fixed point attains

$\mathbb{E}[s_V s_V^\top]$0

Empirically, the paper reports that in a linear RRI system with Bernoulli–Gaussian prior, $\mathbb{E}[s_V s_V^\top]$1, $\mathbb{E}[s_V s_V^\top]$2, $\mathbb{E}[s_V s_V^\top]$3, $\mathbb{E}[s_V s_V^\top]$4, $\mathbb{E}[s_V s_V^\top]$5, and batch $\mathbb{E}[s_V s_V^\top]$6, the SC-VAMP MSE trajectory closely matches state-evolution predictions. For a correlated prior with learned pairwise score, it reports convergence in approximately $\mathbb{E}[s_V s_V^\top]$7 iterations and a final MSE of $\mathbb{E}[s_V s_V^\top]$8 for SC-VAMP versus $\mathbb{E}[s_V s_V^\top]$9 for state evolution, with an approximately $$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$0 gap attributed to score learning and finite-sample effects (Wadayama et al., 11 Jan 2026).

4. Schur-complement SJIM in causal discovery

In causal discovery, SJIM becomes a structural operator on the observational density of an additive noise model

$$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$1

under homoscedastic Gaussian noise $$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$2 and a faithfulness/non-degeneracy assumption requiring

$$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$3

for every true edge $$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$4. The expected SJIM is

$$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$5

The paper’s central theorem states that the diagonal of $$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$6 encodes leaf-node identifiability:

$$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$7

Consequently, a node $$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$8 is a topological sink if and only if $$\nabla_\theta s(\theta)^\top=\nabla_\theta^2 \ell(\theta)$$9, in which case

$$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$0

Thus, acyclicity appears as a geometric signature in the diagonal energies of the SJIM (Wu et al., 28 Apr 2026).

The elimination step is algebraic. For a block partition

$$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$1

the Schur complement is

$$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$2

In linear ANMs, the paper proves exact equivalence between graph marginalization and Schur complementation of the SJIM. If $$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$3 is a leaf, then the SJIM of the marginal $$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$4 is

$$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$5

For Gaussian models, this is also the marginal precision matrix. The proposed Score–Schur Topological Sort therefore extracts a topological order by repeatedly selecting

$$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$6

and updating the active SJIM by a Schur complement, bypassing non-convex acyclicity penalties (Wu et al., 28 Apr 2026).

The nonlinear case is more delicate because expectation and inversion do not commute. The paper states that for a leaf $$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$7,

$$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$8

but also defines an expectation gap

$$\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)\equiv I_{X_{\mathrm{in}}\mid Y} = \mathbb{E}\!\left[s(X_{\mathrm{in}}\mid Y)s(X_{\mathrm{in}}\mid Y)^\top\right],$$9

This gap is localized to the $$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$0 block. The paper proposes Block-SSTS to reduce cumulative error by grouping near-parallel leaves, with a Tikhonov-regularized block Schur update

$$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$1

Empirically, the paper reports exact recovery with zero edge violations up to $$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$2 in linear ANMs using both empirical and population precision matrices. In nonlinear ANMs at $$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$3, SSTS achieves $$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$4 in $$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$5 total, while DAGMA requires approximately $$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$6 with $$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$7. On Sachs data with $$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$8, SSTS achieves $$\mathrm{Tr}\,\mathrm{SJIM}(X_{\mathrm{in}}\mid Y)=J_{X_{\mathrm{in}}\mid Y}.$$9, $I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$0, and total time approximately $I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$1. The same paper also states that post-nonlinear models are a failure mode because the outer nonlinearity couples noise and parents in gradient space, breaking the independence structure used by the SJIM (Wu et al., 28 Apr 2026).

5. Pullback SJIMs in DAG optimization and score chaining

A different line of work uses score-Jacobian constructions in parameter or scene space rather than directly on the sample density. In stochastic DAGs with end-to-end mapping

$I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$2

the mutual-information gradient with respect to parameter block $I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$3 is

$I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$4

Defining

$I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$5

the paper naturally induces a parameter-space matrix

$I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$6

Equivalently, if

$I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$7

then

$I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$8

The paper states that this matrix is positive semidefinite and generalizes Fisher information via the pullback $I=\mathbb{E}_{p(x)}[-\nabla_x^2 \log p(x)],$9. In the single-stage additive Gaussian case where the conditional term vanishes in expectation, $s(x)=\nabla_x \log p(x).$0 reduces to $s(x)=\nabla_x \log p(x).$1 and the matrix becomes the pullback of the marginal Fisher information of $s(x)=\nabla_x \log p(x).$2 (Wadayama, 5 Jan 2026).

This pullback formulation is computationally aligned with automatic differentiation. The same paper emphasizes that the information gradient requires only vector-Jacobian products, implemented through the scalar VJP loss

$s(x)=\nabla_x \log p(x).$3

for which

$s(x)=\nabla_x \log p(x).$4

A batch SJIM estimator is then

$s(x)=\nabla_x \log p(x).$5

The paper reports that DSM-trained and Stein-calibrated scores accurately reproduce analytic gradients in a linear multipath DAG and coincide with finite-difference gradients in the nonlinear scalar channel $s(x)=\nabla_x \log p(x).$6 (Wadayama, 5 Jan 2026).

Score Jacobian Chaining applies an analogous pullback principle to generative 3D optimization. If $s(x)=\nabla_x \log p(x).$7 is a differentiable renderer with Jacobian

$s(x)=\nabla_x \log p(x).$8

and if a pretrained diffusion model provides a 2D score

$s(x)=\nabla_x \log p(x).$9

then the chained 3D score is

$\mathbb{E}[-\nabla_x^2 \log p(x)]$00

and the aggregated score is

$\mathbb{E}[-\nabla_x^2 \log p(x)]$01

The induced SJIM is

$\mathbb{E}[-\nabla_x^2 \log p(x)]$02

where $\mathbb{E}[-\nabla_x^2 \log p(x)]$03 is the 2D score second moment. The paper describes this as a Fisher-information analogue induced by the renderer and the diffusion score, and notes that its nullspace corresponds to 3D perturbations invisible to the current views and renderer (Wang et al., 2022).

A technical complication is distribution mismatch: feeding a clean render directly into a denoiser trained on Gaussian-corrupted inputs is out-of-distribution. The paper therefore proposes Perturb-and-Average Scoring,

$\mathbb{E}[-\nabla_x^2 \log p(x)]$04

and proves that it computes the gradient of a lower bound to $\mathbb{E}[-\nabla_x^2 \log p(x)]$05. The corresponding chained score and SJIM use this PAAS-adjusted image score (Wang et al., 2022).

These pullback constructions suggest a broader interpretation: SJIM can be viewed as the second moment of a pathwise score signal after transport through the system Jacobian. This interpretation is explicit in the DAG information-gradient paper and implicit in score-Jacobian chaining. It places SJIM at the intersection of Fisher geometry, backpropagation, and score-based modeling (Wadayama, 5 Jan 2026, Wang et al., 2022).

6. Derivative-free estimation and large-scale computation

In classical likelihood-based inference, SJIM often coincides with the Hessian of the log-likelihood or the observed information up to sign. The derivative-free estimation paper formulates this explicitly. For $\mathbb{E}[-\nabla_x^2 \log p(x)]$06,

$\mathbb{E}[-\nabla_x^2 \log p(x)]$07

Using an artificial Gaussian prior $\mathbb{E}[-\nabla_x^2 \log p(x)]$08 and the artificial posterior

$\mathbb{E}[-\nabla_x^2 \log p(x)]$09

Stein’s lemma yields moment identities from which derivative-free estimators follow:

$\mathbb{E}[-\nabla_x^2 \log p(x)]$10

The paper states that the score estimator has bias $\mathbb{E}[-\nabla_x^2 \log p(x)]$11 and the observed-information estimator also has bias $\mathbb{E}[-\nabla_x^2 \log p(x)]$12, while Monte Carlo implementations based on importance sampling or SMC are robust to noisy likelihood estimators (Doucet et al., 2013).

For general state-space models, particle methods estimate the score and observed information through Fisher’s and Louis’ identities. In that literature the matrix estimated via Louis’ identity is the observed information, and the details explicitly note that “Score-Jacobian” typically refers to the Jacobian of the score map, i.e.

$\mathbb{E}[-\nabla_x^2 \log p(x)]$13

The proposed particle method combines KDE shrinkage and Rao–Blackwellisation to achieve $\mathbb{E}[-\nabla_x^2 \log p(x)]$14 cost per time step, with score estimator

$\mathbb{E}[-\nabla_x^2 \log p(x)]$15

and observed-information estimator

$\mathbb{E}[-\nabla_x^2 \log p(x)]$16

The paper reports that the estimator is robust to the bandwidth choice and that for any fixed $\mathbb{E}[-\nabla_x^2 \log p(x)]$17, at the true parameter the expected approximate score is zero (Nemeth et al., 2013).

In REML estimation for linear mixed models, the score-Jacobian again appears as the Hessian of the restricted log-likelihood, with the observed information defined as its negative:

$\mathbb{E}[-\nabla_x^2 \log p(x)]$18

The computational obstacle is that exact Jacobian entries involve trace terms such as $\mathbb{E}[-\nabla_x^2 \log p(x)]$19 and $\mathbb{E}[-\nabla_x^2 \log p(x)]$20, which are expensive in large-scale problems. The paper therefore uses the average information matrix $\mathbb{E}[-\nabla_x^2 \log p(x)]$21, with the splitting

$\mathbb{E}[-\nabla_x^2 \log p(x)]$22

and each $\mathbb{E}[-\nabla_x^2 \log p(x)]$23 entry computed with four matrix-vector multiplications. This replaces costly matrix-matrix trace evaluations by sparse linear solves and structured multiplications, enabling AI-REML for large GWAS-scale models (Zhu, 2016).

Across these settings, the practical message is stable: exact score Jacobians are often computationally prohibitive, and workable SJIM surrogates arise from moment identities, score outer products, Rao–Blackwellisation, Schur complements, or average-information splittings. The specific approximation differs by domain, but the organizing principle remains the same: transform second-order score information into a form compatible with the underlying numerical architecture (Doucet et al., 2013, Nemeth et al., 2013, Zhu, 2016).

7. Assumptions, limitations, and unresolved issues

The first limitation is definitional. SJIM is context-dependent. In some papers it is a Fisher outer product; in others an expected negative Hessian; in still others a score pullback through a system Jacobian or simply the score Jacobian itself. Treating these as interchangeable can obscure sign conventions, conditioning properties, and the space in which the matrix acts. This suggests that any use of the term should specify whether the object lives in sample space, latent space, parameter space, or observation space (Wadayama et al., 11 Jan 2026, Wu et al., 28 Apr 2026, Wadayama, 5 Jan 2026, Wang et al., 2022).

The second limitation is model validity. SC-VAMP requires decoupling, statistical symmetry, and suppression of short-range correlations; the paper states that performance may degrade when decoupling fails, learned scores are biased or poorly calibrated, $\mathbb{E}[-\nabla_x^2 \log p(x)]$24 becomes very small, or priors are non-differentiable such as mixtures with Dirac masses. Its entropic-CLT-based state-evolution argument also leaves open a rigorous characterization of double-sided random mixing for general nonlinear modules (Wadayama et al., 11 Jan 2026).

The third limitation is nonlinear expectation mismatch in causal discovery. The Schur-complement theory is exact in linear ANMs, but in nonlinear systems the gap

$\mathbb{E}[-\nabla_x^2 \log p(x)]$25

accumulates across eliminations. Block-SSTS mitigates this by compressing extraction depth, but does not nullify it. The same paper reports that post-nonlinear models are a failure mode, while multiplicative noise models are surprisingly robust empirically. It also emphasizes that, once non-convex optimization is bypassed, structural fidelity is bounded by the finite-sample estimation variance of the global score geometry (Wu et al., 28 Apr 2026).

The fourth limitation is score quality and drift. The DAG information-gradient framework notes that stale scores bias gradients as $\mathbb{E}[-\nabla_x^2 \log p(x)]$26 and $\mathbb{E}[-\nabla_x^2 \log p(x)]$27 drift with changing parameters, and recommends iterative retraining or conditional score networks. It also remarks that DSM can be unstable in low-noise regimes and that velocity-field alternatives such as flow matching or rectified flow may be preferable. In score-Jacobian chaining, failure to correct distribution mismatch yields misleading gradients; the PAAS construction is introduced precisely because naive denoiser calls on clean renders are out-of-distribution (Wadayama, 5 Jan 2026, Wang et al., 2022).

Finally, large-scale estimation remains numerically delicate. In derivative-free Stein-based estimation, small perturbation covariances reduce bias but worsen weight degeneracy. In particle observed-information estimation, ill-conditioned Hessian surrogates require ridge stabilization or damping. In AI-REML, the average-information approximation is most faithful when the variance-component parameterization is linear and can degrade when second-derivative terms are large. These issues do not invalidate SJIM-based methods, but they do show that the matrix is only as reliable as the score estimates, structural assumptions, and numerical approximations that produce it (Doucet et al., 2013, Nemeth et al., 2013, Zhu, 2016).

Taken together, the literature presents SJIM as a versatile but non-uniform concept: a score-built information operator that can replace explicit Jacobians in message passing, reveal causal hierarchy through Schur complements, pull mutual-information score power back into parameter space, or stand in for observed information in derivative-free inference. Its importance lies less in a single formula than in the recurring idea that score geometry can be converted into tractable second-order structure.