Implement pointobservation

[finsberg]

Motivation

Inverse problems are usually posed against data that lives at points, not on the computational mesh: sensor readings, well measurements, or the voxel centres of an image. Today dolfinx-adjoint can only build functionals from assemble_scalar, i.e. integrals over the mesh, so matching point data means first interpolating the measurements onto a finite element field. That is unfortunate for two reasons:

1. It invents data where nothing was measured. Interpolating a few hundred sensor values onto a fine mesh manufactures exactly the information the inverse problem is supposed to extract.
2. It changes the problem. The resulting functional is a mesh integral of an interpolant, not the least-squares misfit of the measurements, so the gradient is the gradient of a different objective.

This PR adds the pointwise observation operator and a differentiable least-squares misfit built on top of it, so that the parameter-to-observable map $m \mapsto B,u(m)$ of PDE-constrained optimization and Bayesian inversion can be written down directly and taped.

The math

The observation operator

A finite element function is $u(x) = \sum_j u_j \phi_j(x)$, with $\phi_j$ the basis of $V$ and $u_j$ the degrees of freedom. Evaluating it at a point is therefore linear in the degrees of freedom: it is just the basis functions sampled at that point. Collecting $n_d$ observation points $x_1,\dots,x_{n_d}$ into a matrix gives the discrete observation operator

$$ d = B u, \qquad B \in \mathbb{R}^{n_d \times N}, \qquad B_{ij} = \phi_j(x_i), $$

where $N = \dim V$. Row $i$ of $B$ evaluates $u$ at $x_i$. The matrix is extremely sparse: $\phi_j(x_i)$ is nonzero only for the handful of basis functions supported on the cell containing $x_i$.

For a vector-valued space with block size $b$ the same construction applies per component; row $ib + c$ observes component $c$ at point $x_i$.

The misfit

Given measured data $d \in \mathbb{R}^{n_d}$ we minimize

$$ J(u) = \frac{1}{2\sigma^2},\lVert W (B u - d) \rVert^2, $$

with $\sigma^2$ the variance of the additive Gaussian observation noise (use $\sigma^2 = 1$ for a purely deterministic problem) and $W$ an optional diagonal matrix of per-observation weights. A $0/1$ weight masks individual sensors out; a general $W$ lets you use $\sigma^{-2} W^2 = \Gamma_{\text{noise}}^{-1}$ for a diagonal noise covariance, which is what makes $J$ the negative log-likelihood of the data in the Bayesian setting.

Derivatives

$J$ is a quadratic function of the state, so no linearization point has to be stored and every derivative is available in closed form. Writing $r = Bu - d$:

$$ \frac{\partial J}{\partial u} \delta u = \frac{1}{\sigma^2}, (W r)^{!\top} W B, \delta u \quad\Longrightarrow\quad \nabla_u J = \frac{1}{\sigma^2}, B^{\top} W^2 r , $$

$$ \nabla_u^2 J = \frac{1}{\sigma^2}, B^{\top} W^2 B . $$

Three things fall out of this:

• The adjoint (reverse) action is $\bar u = \bar J,\sigma^{-2} B^{\top} W^2 r$: apply $B$, subtract the data, weight, and apply $B^{\top}$. Since $B$ evaluates, $B^{\top}$ scatters: it deposits $\phi_j(x_i)$-weighted contributions from each observation point back onto the degrees of freedom of the cell that contains it. This is a discrete sum of Dirac deltas, which is precisely the right-hand side one gets by hand when deriving the adjoint equation for point data. - The tangent-linear (forward) action of $J$ is the scalar $\sigma^{-2}(Wr)^{!\top} W B \hat u$.
• The Hessian is exact and equals the Gauss–Newton operator — there is no second-order term to drop, because $B$ is linear. The Hessian action of the composed problem is still only Gauss–Newton in $u$; the PDE contributes its own curvature through the rest of the tape as usual.

Because the Hessian never depends on $u$, the block stores nothing but $d$, $\sigma^2$, and $W$.

Implementation

PointObservation (src/dolfinx_adjoint/observation.py)

Builds $B$ once and exposes its action:

B = dolfinx_adjoint.PointObservation(V, points)   # points: (num_points, gdim)
values = B.evaluate(u)                            # u at every point, on every rank

The matrix itself is assembled by fenicsx_ii. The key observation is that we do not have to write a point-evaluation kernel at all: a DG-0 space on a point mesh (one cell per observation point, dolfinx.mesh.create_point_mesh) has exactly one degree of freedom per point, so the interpolation matrix from $V$ onto that space is $B$. fenicsx_ii.create_interpolation_matrix assembles it as a distributed PETSc.Mat, including all cross-process communication, and its transpose gives $B^{\top}$ for free.re is a small _PointCloudTrace reduction operator that hands fenicsx_ii the coordinates: its PointwiseTrace is written for 1D line meshes and gets the physical coordinates by compiling a SpatialCoordinate expression, which FFCx cannot do on a point cell — on a point mesh each cell is a geometry node, so the coordinates are read straight off the geometry.

This adds fenicsx-ii as a dependency.

The misfit.

J = dolfinx_adjoint.point_observation_misfit(u, B, data, noise_variance=sigma2, weights=W)

returns a pyadjoint.AdjFloat and records a PointObservationBlock on the tape, so it composes with LinearProblem/NonlinearProblem and adds to assemble_scalar terms (e.g. Tikhonov regularization) exactly like any other functional. The block implements recompute, evaluate_adj, evaluate_tlm and evaluate_hessian from the closed forms above. data and weights are copied on construction, so a caller reusing one buffer across a time-stepping loop cannot retroactively change a block the tape is still holding. In the Hessian, the second-order seed and the curvature term are summed in row space before $B^{\top}$ is applied, so the communication round-trip runs once per Hessian action rather than twice.

Parallel semantics

This is where most of the care went.

• points must be replicated on every rank — they typically come from a file or an instrument layout every rank can read. This is checked with a length reduction plus a position-weighted checksum, so equally many but different (or merely permuted) points are rejected rather than silently corrupting the operator.
• Each point is assigned to exactly one owner: only owned cells are searched, and the remaining ties on shared facets and vertices are broken by lowest rank.
• Points no rank can locate are reported in B.found / B.owner and excluded from the operator rather than becoming zero rows. A zero row silently contributes $d_i^2$ to the misfit and biases the result.
• The default bounding-box padding is computed from the global bounding box, so the same points are located regardless of how the mesh happens to be partitioned.
• Validation that can fail on one rank only (bad points shape, mixed elements, elements needing DOF transformations) is checked collectively before anyone raises. Raising locally would let one rank leave the constructor while the others block forever in the next collective call.

apply/apply_transpose work in the distributed row layout; gather/restrict/evaluate convert to and from a replicated global array, with a fill value (default nan) for points outside the mesh.

Testing

37 tests in tests/test_observation.py, run in serial and in parallel:

• Correctness of $B$: exact on polynomials up to the element degree, agrees with a manual dolfinx point evaluation, rows sum to one (partition of unity), component-fastest ordering for vector spaces, simplices/quads/hexes including distorted ones, and the generic pull-back path.
• Parallel behaviour: unique ownership agreed on by all ranks, ambiguous points on shared entities owned exactly once, partition-independent padding, ranks owning no points, an empty point set, and the deadlock-free error paths.
• Adjointness: $\langle Bu, v\rangle = \langle u, B^{\top}v\rangle$ to machine precision.
• Derivatives: gradient against the closed form, weights applied twice, the Hessian equal to $\sigma^{-2}B^{\top}W^2B$, and Taylor tests both for the bare misfit and through a full PDE solve.
• End-to-end: recovery of a known source from point data.

Demo

demos/point_observations.py (added to the toc) solves the mother problem against noisy sensor data on a jittered $12\times 12$ grid: build $B$, generate synthetic obse+ Tikhonov with L-BFGS-B, verify the gradient with a Taylor test, and compare the recovered source against the truth. The sensor residual settles at the noise floor rather than below it, which is the outcome to hope for — driving it to zero would mean fitting the noise.

AI assistance

I used Claude Code (Claude Sonnet 5 and Claude Opus 5) to help implement, test, and iterate on this feature, and to draft this PR description. I reviewed, tested, and take full responsibility for the final contribution.