arXiv CSopen access

Perturbative Contrastive Physical Learning

Perturbative Contrastive Physical Learning Kyungeun Kim,1, ∗ Amanuel Anteneh,2 Israel Klich,2, 3 Olivier Pfister,2, 4, † and J. M. Schwarz5, 6, ‡ 1

arXiv:2606.09756v1 [cs.LG] 8 Jun 2026

Department of Mathematics, University of British Columbia, Vancouver, BC Canada 2 Department of Physics, University of Virginia, 382 McCormick Rd, Charlottesville, VA 22903, USA 3 Max Planck Institute for the Physics of Complex Systems, 01187 Dresden, Germany 4 Charles L. Brown Department of Electrical and Computer Engineering, University of Virginia, 351 McCormick Road, Charlottesville, VA 22903, USA 5 Department of Physics, Syracuse University, Syracuse, NY 13244, USA 6 Indian Creek Farm, Ithaca, NY 14850, USA (Dated: June 9, 2026)

Responses to perturbations are key to understanding physical systems. The ability to contrast such responses by comparing how a system reacts under slightly different conditions provides a mechanism for learning. Here, we introduce Perturbative Contrastive Physical Learning (PCPL), a general framework in which learning emerges from measurable contrasts between physical states produced by controlled changes to inputs, boundary conditions, parameters, or interpreter functions. PCPL unifies and extends prior approaches: Equilibrium Propagation is rooted in contrasts between free and nudged equilibria in energy-based systems, while Frequency Propagation corresponds to contrasts extracted from sinusoidally driven, frequency-demodulated responses. We show that contrast-driven updates can reflect either local sensitivities or global inverse-problem structure, yet do not require centralized gradient computation. Instead, effective learning geometry emerges implicitly from the system’s own physical response, allowing learning behavior to arise without an external processor or explicit backpropagation. We demonstrate PCPL in two platforms: (i) spring networks that update bond stiffness using measured displacements and forces, and (ii) continuousvariable photonic circuits trained via x quadrature measurements and finite-difference estimates of the Jacobian. Both platforms successfully learn classification tasks. We further show that a continuous-variable photonic circuit can be trained to implement analog multiplication, illustrating a step toward more autonomous physical learning systems.

I.

INTRODUCTION

Learning requires the ability to modify behavior based on experience. In biology, learning is often associated with brains and centralized neural processing [1, 2]. Yet organisms without nervous systems, such as slime molds, are also capable of learning [3]. Moreover, nonliving brain-inspired systems, such as convolutional neural networks, also learn [4]. These observations suggest that learning may emerge from general principles that can be realized directly in physical processes. With this viewpoint, learning arises from measurable changes in system behavior produced by controlled perturbations of inputs, boundary conditions, or internal parameters [5–8]. Physical learning, therefore, represents a fundamental shift from brain-inspired learning at the software level to learning that emerges directly from physical processes, with the hardware effectively becoming the software [9]. A key physical principle underlying such learning is the response of a system to perturbations. Phys-

∗ [email protected] † [email protected] ‡ [email protected]

ical systems reveal their structure through how observable quantities change when inputs, boundary conditions, or internal parameters are varied. The ability to contrast these responses naturally provides a mechanism for learning. If a system can be perturbed in a controlled way and the resulting changes in its measurable outputs can be detected, those contrasts can be used to guide updates to modifiable parameters. Learning then becomes a process of probing, measuring, and remodeling a physical system based on how it responds. Here, we introduce Perturbative Contrastive Physical Learning (PCPL), a general framework that formalizes this principle. In PCPL, parameter updates are driven by measurable contrasts between physical states generated by controlled perturbations. These perturbations may involve changes to inputs, boundary conditions, internal parameters, or readout mappings, but in all cases learning is derived from the difference between two nearby physical responses. The perspective that physical learning can emerge from contrasts between nearby physical states was crystallized by Equilibrium Propagation (EP), which showed that measurable differences between free and weakly nudged equilibria can generate learning signals equivalent to gradient-based optimization [5, 10]. PCPL builds on this insight by

2 treating contrast itself as the central primitive. And rather than requiring an explicit energy function, equilibrium convergence, or symbolic gradient computation, PCPL treats learning as an experimentally accessible operation: apply a perturbation, measure the contrast in observables, and translate that contrast into a parameter update. This perspective reveals that several existing physical learning schemes can be understood as special cases of a broader contrast principle. In addition to EP [5, 10], there is Coupled Learning with its experimental realizations [11–15] and Multimechanism learning with the contrast encoded in different modalities [16]. Frequency Propagation extracts contrasts from sinusoidally driven, frequencydemodulated responses [17]. In addition, explicit temporal contrasts between free and nudged states using a non-equilibrium memory kernel also make learning possible [18]. PCPL generalizes this idea beyond equilibrium settings and energy-based formulations, allowing contrast-driven learning to be implemented in systems with complex dynamics, multiple steady states, or black-box measurement and interpretation layers. By elevating the contrast itself to the central primitive, PCPL provides a unifying framework for learning across physical substrates. These ideas are not limited to classical platforms. Recent work has shown that contrast-based learning rules can also be implemented in quantum systems, where observables replace classical state variables and measurement outcomes provide the signals needed for adaptation. In particular, quantum generalizations of Equilibrium Propagation demonstrate that contrasts between weakly perturbed quantum steady states can be used to extract local update signals through measurable expectation values and correlations [19, 20]. From the PCPL perspective, such schemes represent a specific instance of a broader contrast principle: learning is driven by observable differences between nearby quantum states. Continuous-variable photonic platforms offer a natural testbed for this view as Gaussian operations, which include quantum squeezing, singlephoton gain/loss, interferometric couplings, and homodyne detection, enable controlled perturbations and precise readout of system responses—including at the quantum computing level [21–25]. In this way, quantum physical learning appears not as a fundamentally different paradigm, but as a direct extension of contrast-driven adaptation into regimes where information is encoded in quantum states and accessed through measurement. While such quantum realizations are an important and rapidly developing direction, the present work focuses on classical mechanical and optical systems, leaving a detailed

exploration of quantum platforms to future studies. The organization of the manuscript is as follows. First, we introduce the overall framework, then we apply the framework to a classical spring network followed by several photonic examples. We conclude with a discussion containing a summary and both experimental and theoretical implications.

II. PERTURBATIVE CONTRASTIVE PHYSICAL LEARNING: A TWO MODE FRAMEWORK

We begin with a dataset D of paired observations (x, z), where x denotes the set of inputs and z denotes the corresponding set of outputs. In physical learning systems, z represents a physical readout, such as a desired displacement or stress pattern. A batch refers to a finite subset of such pairs sampled from D. The mapping from inputs x to measurable outputs z may involve an initial physical encoding, or preprocessing, stage y = F (x) followed by a parameterized physical response G(y, θ) = z. Here, G denotes the observable response of the physical system under modifiable parameters θ to enable learning. Throughout this work, learning acts on the response map G, while F is treated as fixed. This decomposition allows fixed physical architectures to solve more complex tasks through structured input representations rather than architectural expansion. We consider two related learning modes that update parameters θ of G, the response of a physical system, from weak nudging perturbations. In both cases, a batch is sampled, features are computed, a parameter update ∆θ is applied via θ ′ = θ + ∆θ. Note that throughout the paper we denote by δθ a small imposed probe perturbation used to generate a measurable contrast, and by ∆θ the parameter update inferred from that contrast. In simple implementations these may be proportional, but they play distinct conceptual roles. To determine ∆θ, the resulting response G(y, θ ′ ) is used to form a contrast measurement functional M . This functional is task- and implementation-dependent. We will consider learning rules based on measurable contrasts between nearby physical responses via observables z1 and z2 such that ∆θ ∼ M (z1 , z2 ).

(1)

In other words, M maps observable contrasts into parameter updates. Specific forms of M will be discussed below. We note that it need not represent an analytically computed gradient; rather, it encodes how physically measurable contrasts can be used for learning.

3

FIG. 1. Schematic for PCPL with Mode A (implicit target) and Mode B (explicit target).

The distinction between the two modes is the reference against which the nudged response is compared: Mode A contrasts the nudged response against the system’s own free response, while Mode B contrasts it against an externally specified target, often realizable as a physical constraint. Let us now be more specific about the two modes.

A.

Mode A: Self-referenced contrast: Linear response learning

In Mode A, learning is driven by measurable contrasts between free and weakly perturbed, or nudged, physical responses, enabling parameter updates through local response probing. The system is first allowed to relax in a free phase, producing zf = G(y, θ),

(2)

followed by a nudged phase with θ ′ = θ+δθ, yielding zn = G(y, θ + δθ).

(3)

For weak nudges, a first-order expansion gives the local parameter-response map G(y, θ+δθ) ≈ G(y, θ) + Jθ (y, θ) δθ,

Jθ ≡

∂G . ∂θ (4)

Thus, the contrast between phases is ∆z ≡ zn − zf ≈ Jθ δθ,

(5)

which is directly accessible as a physical measurement.

The Mode A learning signal depends on the phase contrast relative to target z∗ via M (zn , zf )|z∗ and is naturally interpretable as a Jacobian-driven update. If M extracts a correlation between internal degrees of freedom and the observed change ∆z, then the update inherits the structure of linear response theory: one probes the system, measures the response, and uses that response to adapt parameters. This self-referential baseline subtraction also reduces bias from offsets or slow drifts in the operating point, making Mode A especially natural for physical systems where free relaxation, such as energy minimization, is meaningful. Mode A does not require the target to appear directly in the perturbation step. Instead, the target enters implicitly through the choice of M . Below is the pseudo code for learning using Mode A. Perturbations can be implemented using either one-sided or symmetric nudging. In onesided nudging, parameters are shifted as θ → θ + δθ, whereas symmetric nudging applies paired perturbations θ ± δθ, enabling a central-difference estimate of the response that reduces bias and improves accuracy at the cost of additional measurements.

4 Algorithm 1 Mode A: simple nudging Require: Dataset D; feature map F ; system response G; contrast measurement functional M ; initial parameters θ0 ; learning rate η; perturbation rule for δθ; iterations T

Defining the instantaneous output error e ≡ z∗ − z, the linearized goal e ≈ Jθ ∆θ is posed as a leastsquares problem, 2

∆θ = arg min ∥e − Jθ ∆θ∥2 , ∆θ

1: θ ← θ0 2: for i = 0 to T − 1 do 3: (x, z) ← SampleBatch(D) 4: y ← F (x) 5: Free phase: zf ← G(y, θ) 6: Nudge: draw small δθ and set θ ′ ← θ + δθ 7: Nudged phase: zn ← G(y, θ ′ ) 8: Contrast-to-update: corr ← M (zn , zf )|z∗ 9: Parameter update: θ ← θ + η corr 10: end for 11: return θ

whose solution is the Moore–Penrose pseudo-inverse, ∆θ ≈ Jθ+ e,

∆θ ≈

JθT Jθ + λI

−1

JθT e.

(10)

Consequently, Mode B does not merely probe local sensitivity (as in Mode A), but instead performs a local, metric-aware correction step in parameter space, where the effective metric is given by Jθ⊤ Jθ . Learning therefore follows the natural geometry induced by the system’s input–output map rather than the Euclidean geometry of parameter space. Mode B performs target-referenced learning by translating output errors into parameter updates through the local response geometry of the system. Depending on the implementation, these updates may take the form of gradient-descent, Gauss–Newton, pseudoinverse, or other inverse-geometry corrections. However, in the implementations presented here, the Jacobian is measured through controlled perturbations and its pseudoinverse (or a regularized approximation) is computed externally. The key advantage of this metric-aware update is robustness: because (Jθ⊤ Jθ + λI)−1 automatically normalizes the update direction according to the local response geometry, Mode B remains stable across a wide range of learning rates and is insensitive to ill-conditioning of the input space in comparison with gradient descent, for example. Please see Appendix A for more details. Here is the pseudocode for Mode B. Algorithm 2 Mode B: direct inverse Require: Dataset D; feature map F ; system response G; initial parameters θ0 ; learning rate η; iterations T

Mode B: Target-referenced contrast: Anchored Learning

Unlike Mode A, which compares nearby responses around the system’s free state, Mode B compares the system’s response to an externally specified target or imposed constraint. Assume we want to modify the weights θ so that the response to a given y is z ∗ . Linearizing the actual response about θ gives G(y, θ + ∆θ) ≈ G(y, θ) + Jθ ∆θ, z∗ ≈ z + Jθ ∆θ.

(9)

or with Tikhonov regularization,

Mode A provides a broader perturbative contrast framework that generalizes several existing physical learning algorithms. In Equilibrium Propagation (EP) [5], learning arises from contrasts between free and weakly nudged equilibrium states in energy-based systems. Coupled Learning similarly extracts update signals from measurable differences between nearby physical responses with experiments implementing two copies of the system to determine differences[11], while Multi-mechanism learning eliminates the need for two copies by using two types of non-interfering responses, i.e. mechanical and electrical should they be decoupled [16]. Frequency Propagation considers different frequency domains [17]. Mode A unifies these approaches by treating learning as arising from measurable contrasts generated by controlled perturbations, independent of whether the underlying system is equilibrium-based, dynamical, classical, or quantum. In this view, the central object is not an explicit gradient or energy function, but the experimentally accessible response difference that acts as a physically measurable probe of the local response geometry of the system.

B.

(8)

(6) (7)

1: θ ← θ0 2: for i = 0 to T − 1 do 3: (x, z) ← SampleBatch(D) 4: y ← F (x) 5: Reference output: z ← G(y, θ) 6: Measure J: vary θ and evaluate Jθ by measur-

ing system response z Target mismatch: e ← z∗ − z Parameter update: θ ← θ + η Jθ† e (e.g. pseudoinverse or other invertible forms) 9: end for 10: return θ 7: 8:

5 While Algorithm 2 presents Mode B in terms of an explicit local inverse, an alternative implementation can be constructed using perturbative contrasts and a measurement functional. In this approach, parameter updates are inferred from correlations between nudged responses and target outputs, allowing the effective inverse mapping to emerge from physically measurable quantities without requiring explicit evaluation of the Jacobian or its pseudoinverse. C.

Two Learning Geometries: Sensitivity Probing vs. Local Inversion

Mode A is most directly physical when the free response zf is a well-defined operating point (or steady state) and nudging corresponds to a small, controlled perturbation: the learning signal is then a bona fide linear response measurement, Eq. (5). Mode B becomes comparably physical when the target z∗ is not a symbolic label but is embodied as a constraint—a boundary condition, applied force pattern, clamped voltage, imposed concentration field, or geometric confinement. In that regime, the “target” is simply a physical condition the system is required to satisfy, and learning corresponds to remodeling internal parameters so that the constrained configuration becomes achievable (or stable) with progressively weaker external enforcement. A useful practical hybrid is therefore: employ Mode B to imprint or instruct desired behaviors under constraints, then use Mode A to consolidate them by aligning the free dynamics with the nudged behavior via contrastive, Jacobian-driven refinement. Regarding the measurement functional M , it abstracts the physically accessible measurement used to form a learning signal from pairs of system responses. In Mode A, M (zf , zn ) extracts information proportional to the phase contrast ∆G ≡ ∆z = zn − zf ,

(11)

which, in the weak-perturbation limit, approximates the Jacobian action Jθ δθ (Eq. 5). In Mode B, M (zn , z∗ ) encodes a local discrepancy between the nudged response (obtained via random perturbation or a physics-guided mechanism) and the target, which can be interpreted as an error signal driving a local inverse problem (Eq. 8). While Mode B superficially resembles pseudoinverse-based learning rules used in machine learning, the underlying perspective is different. Conventional inverse methods begin with an explicit computational model and solve an inverse problem within that model. In contrast, Mode B infers the local learning geometry from

physically measurable perturbations and response contrasts. In the proof-of-principle implementations presented here, this geometry is represented through an experimentally estimated Jacobian and an externally computed pseudoinverse update. As a result, Mode B can be interpreted as a physically realized pseudoinverse update, where the geometry of the parameter space is encoded in the system’s response to perturbations and read out through contrasts, rather than computed analytically, thereby realizing numerical optimization through physical response measurements. This distinction is crucial: the pseudoinverse update is an abstract linear-algebraic procedure, whereas Mode B implements metric-aware learning directly through physically measurable contrasts. Finally, Mode B is conceptually related to target propagation [26]; however, target propagation constructs desired hidden-state targets using approximate inverse mappings, whereas Mode B directly solves a local inverse problem in parameter space using the measured response geometry of the physical system.

III.

PHYSICAL REALIZATIONS OF PCPL IN MECHANICAL AND PHOTONIC SYSTEMS

We now present two concrete realizations of PCPL: a classical mechanical network and a continuous-variable photonic circuit. These examples are not intended as domain-specific optimizations, but as conceptual demonstrations of how the same perturbation–contrast principle gives rise to learning across different physical substrates. In particular, we emphasize how both systems naturally implement Mode B (target-referenced) learning through pseudoinverse or Gauss–Newton–like geometry, while remaining compatible with Mode A–style linear-response probing. See Appendix B for more details on PCPL.

A.

Spring Networks as Contrastive Learners

The Network. Consider a spring network with modifiable spring constants and equilibrium spring lengths. An illustrative example is shown in Fig. 2. Here, the blue nodes are the input ports, where datadependent forces are applied, while the green nodes are fixed anchors that impose the boundary conditions. The rest of the network relaxes in response to these forces, producing a displacement field that serves as the physical response. The learning objective is to adjust mechanical parameters so that the displacement field describes a desired output given

6 an input. The network will be trained to perform classification. We begin with the energy functional: E= =

1X 2 ij

gives  2 −1 ∂ ∂u∗ ∂ E =− (∇u E) , ∂kij ∂u2 ∂kij

(0)

kij (lij − lij )2

X 1X 1X 2 Aij lij + Bij lij + Cij 2 ij 2 ij ij (0)

(12)

(0)2

where Aij = kij , Bij = −kij lij , Cij = kij lij . Given this energy functional, two training issues (0) arise. First, to train rest length lij and stiffness kij , we need to update A, B. Since kij contributes to both A and B, this creates an identifiability prob(0) lem: multiple combinations of kij and lij can produce the same energy, making the optimization illconditioned and prone to overfitting.

which involves the inverse stiffness/Hessian matrix. To remain consistent with the PCPL philosophy and avoid this global adjoint calculation, we instead estimate parameter sensitivities from measured perturbative responses and construct updates from the resulting Jacobian. For the proof-of-principle example considered here, we restrict attention to a small linear spring network, where this Jacobian can be computed efficiently from the output displacements. Given the first identifiability issue, we introduce coordinate vectors ui : setting the rest lengths as (0) fixed to lij (ui , uj ) = 1 and add a gravitational force term. Thus, the energy equation becomes: E=

X 1X kij (lij (ui , uj ) − 1)2 + g mk hk (uk ) 2 ij k

1X kij (lij (ui , uj )2 − 2lij (ui , uj ) + 1) = 2 ij X +g mk hk (uk )

(13)

k

FIG. 2. Spring Network Schematic. The green nodes are fixed and the 4 blue nodes denote where the input forces containing the Iris dataset information are applied.

The second issue concerns gradient evaluation. In a mechanical network, the equilibrium node positions depend implicitly on all spring parameters through force balance, ∇u E(u, k) = 0. Thus, although a partial derivative such as 2 1 ∂E (0) = lij − lij ∂kij 2 appears local, changing kij shifts the equilibrium configuration u∗ (k). More generally, for an objective L(u∗ (k), k), dL ∂L ∂L ∂u∗ = + · , dkij ∂kij ∂u ∂kij so exact gradient computation requires the global sensitivity of the relaxed configuration to parameter changes. Linearizing the force-balance condition

where hk (uk ) represents the vertical height of node k. This approach enables us to separate the quadratic and linear terms, allowing us to train both spring constants kij and masses mk simultaneously using the Jacobian pseudo-inverse method. To address the complexity of gradient computation, we restrict to a small linear spring network (Fig. 2), in which parameter updates are computed via the Jacobian of the output displacements with respect to the spring constants. Iris dataset features as force inputs.We now demonstrate how a mechanical system equipped with a quadratic readout can implement targetreferenced contrastive learning using a pseudoinverse update rule. While the classification task itself is standard, the purpose of this example is not to optimize accuracy but to show how a physical system can internalize a nonlinear decision boundary through Mode B geometry. We use a 6 × 6 triangular lattice as a linear spring network. The Iris features are encoded as dimensionless input vectors f̃ ∈ R4 , which act as abstract forcing patterns rather than physical loads. Instead of interpreting these inputs as displacements or forces with physical units, we treat them as coordinates in a normalized input space and define a dimensionless objective function that couples the network’s mechanical response to a learned quadratic form. Specifically, the spring network produces a physical elastic energy Espring ,

7 which we nondimensionalize by a reference energy scale E0 , Ẽspring =

Espring . E0

(14)

We then define a dimensionless quadratic functional F (f̃ ) = f̃ ⊤ W f̃ + b⊤ f̃ + c,

(15)

where W , b, and c are dimensionless parameters learned during training. The total dimensionless readout is F̃total = F (f̃ ) + αẼspring ,

(16)

where α is a dimensionless coupling coefficient that controls the relative contribution of the mechanical response. The two-dimensional feature space is (F̃total , d˜max ), where d˜max is the normalized maximum displacement, is used for classification. Training using PCPL Mode B. Training is performed using a pseudoinverse learning rule. In particular, the Jacobian takes the form   J = vec(f̃ f̃ ⊤ ), f̃ , 1 . (17) To then determine the update ∆θ, we use Tikhonov regularization in Equation 10. Finally, the parameters are updated using a learning rate η as θnew = θ + η ∆θ,

(18)

which is used to update the parameters W , b, and c. In each epoch, the spring network is solved using the current input vector, producing Ẽspring , while the quadratic functional is evaluated using the same input. Classification thresholding. We evaluated three decision strategies for classification, which differ in how the learned mechanical responses are mapped to discrete class labels: • Case 1: Scalar level matching in Ftotal . Each sample is assigned to the class whose target value tc is closest to the observed scalar readout Ftotal , i.e., arg min |Ftotal − tc | . c

Here, the tc are class-specific reference values learned during training. Using this nearestscalar readout rule, the system achieved 90.67% test accuracy after 500 epochs with learning rate η = 0.001. • Case 2: Adaptive boundary classification in scalar–displacement space. To

map mechanical responses to discrete behavioral classes, we also employ an adaptive twothreshold classifier derived directly from the distribution of training data. For each input force pattern, the model computes a combined scalar response Etot = EFEM +f ⊤ W f +b·f +c, where EFEM is the physical strain energy from the lattice simulation and the quadratic form represents the learned energetic contribution. In parallel, we extract the maximum nodal displacement Dmax from the mechanical solution. Class boundaries are then determined adaptively from training-set extrema rather than fixed parameters. Specifically, an energy threshold
E ∗ = 12 min E0 + max(E1 , E2 ) separates class 0 from the union of classes 1 and 2. Within the low-energy regime Etot ≤ E ∗ , a second threshold
D∗ = 21 max D1 + min D2 distinguishes class 1 from class 2. The resulting decision rule is therefore hierarchical: energy first partitions responses into “highenergy” versus “low-energy” behaviors, after which displacement resolves the low-energy cases into distinct mechanical states. Because both thresholds are recomputed from the current training distribution, the decision boundary adapts automatically to shifts in system stiffness, geometry, or learned parameters. • Case 3: Mechanically grounded scalar classification. In this strategy, the learned quadratic readout is removed, and classification is based solely on a scalar summary of the spring network’s mechanical response. For each input force pattern, the network is solved to mechanical equilibrium, producing a displacement field u. From this state, we compute a dimensionless scalar measure of deformation, Smech = 12 u⊤ Ku, which reflects the global elastic response of the lattice rather than a direct algebraic function of the inputs. Because the equilibrium displacements arise from a coupled force-balance problem, the mapping from input features to Smech is geometrically nonlinear, even though the energy expression itself is quadratic in u. Class boundaries are determined adaptively from gaps in the training distributions of Smech . Specifically, two scalar thresholds are defined as

b1 = 21 max S0 +min S1 , b2 = 12 max S1 +min S2 ,

8 where Si denotes the set of mechanical scalar responses for class i. The resulting decision rule is therefore a two-level threshold classifier operating in a one-dimensional, mechanically generated feature space. Because the thresholds are recomputed from the current training distribution, the classifier adapts automatically to changes in stiffness, geometry, and learned parameters while remaining grounded in physically measurable quantities. Performance. Figure 3 shows the decision boundary of Case 2 and Case 3, which both yielded higher accuracy than Case 1. The difference in classification accuracy across Cases 1–3 reflects the geometric expressivity of the feature spaces rather than merely the number of thresholds. Case 1 relies on a single scalar readout with nearest-level matching, which restricts the decision rule to affine separations and yields lower accuracy. In contrast, Cases 2 and 3 both employ two-threshold decision rules, but operate in nonlinear feature spaces generated by the physical system. Case 2 uses a two-dimensional feature space (Etot , Dmax ), enabling separations that exploit independent mechanical response modes and thereby achieving the highest accuracy. Case 3, while also nonlinear, compresses all class information into a single mechanically generated scalar, limiting separability when class clusters overlap along that axis and resulting in slightly reduced performance relative to Case 2. Interestingly, using a spring network physical learning system with linear classification, prior work reports an average classification accuracy of 96% on the Iris dataset using a 70/30 train–test split [27]. Under slightly more challenging splits, such as 50/50, this performance may be expected to decrease further due to both the reduced amount of training data and the restricted representational capacity of linear models, which can only implement affine input–output mappings of the form Ay + b = z. Such mappings are unable to represent curved or topologically nontrivial decision boundaries. This limitation is reflected in our Case 1, which relies on a scalar-level matching rule in a one-dimensional readout space and therefore inherits the expressivity constraints of linear classifiers. As expected, geometrically nonlinear decision rules can achieve substantially higher performance on this dataset. Indeed, a wide range of nonlinear models report test accuracies approaching ∼ 98% [28], which appears to represent a practical generalization ceiling imposed by intrinsic class overlap. Our Case 2, which constructs adaptive decision boundaries in the two-dimensional (Ftotal , D) feature space, falls into this category: by allowing curved and data-adaptive separations, it achieves

FIG. 3. Decision boundaries of case 2 and case 3 for a classical spring network in a gravitational field: Left - Hybrid system classification using adaptive thresholds in Ftotal vs. displacement space achieves 100% test accuracy. Right - Pure mechanical system classification using Etotal (linear strain energy) with gap-based energy thresholds achieves 96% accuracy, where both stiffness matrix K and masses b are trained via Jacobian pseudoinverse with parameter perturbation.

100% test accuracy for particular train–test partitions. Importantly, such perfect separation should be interpreted as a finite-sample geometric effect rather than evidence that the intrinsic limit has been exceeded. For specific splits, the learned boundaries may align favorably with the data geometry, but this behavior is not expected to persist under repeated randomized splits or cross-validation. Finally, Case 3 introduces an additional layer of nonlinearity through the mechanically-mediated mapping between inputs and scalar readouts, allowing the classifier to exploit structure induced by the spring network itself. While still subject to

9 the same intrinsic data overlap constraints, this mechanically grounded nonlinearity provides a distinct route to enhanced expressivity that differs fundamentally from purely algebraic nonlinear classifiers. Together, these results highlight the tradeoff between computational simplicity and geometric expressivity, and clarify the distinction between splitspecific performance and intrinsic generalization limits.

B.

Continuous Variable Photonic Circuits as Contrastive Learners

We now show that the same PCPL framework applies to continuous-variable (CV) photonic systems. For convenience, we employ the standard Gaussian-state formalism of continuous-variable optics and represent optical modes using bra-ket notation. Because the examples considered here remain entirely within the Gaussian sector, they also admit an equivalent description in terms of classical optical field amplitudes. The learning logic, however, remains unchanged: perturb, contrast, and update. Let us build continuous variable photonic circuits out of the following components: displacement gates D(α), phase shift, a.k.a. phase-space rotation, gates R(ϕi ), beam splitter gates, BS(θ, ϕ) to create two-mode couplings, and squeezing gates. The circuit will contain modifiable parameters so that it can learn to classify objects, for example. Using our learning notation, F prepares photonic states (F : x 7→ |ψ(x)⟩), G applies parameterized operations: G : |ψ(x)⟩ , θ 7→ |ψ(x; θ)⟩. As the parameter contrast in Mode A [29] is given by ∆θ = M (ψfree , ψnudged ) depending on the readout, the update becomes one of the following,

1.

Gaussian CV photonic (using linear optical gates) circuit

This example demonstrates how Mode B PCPL can be implemented in a CV photonic circuit to perform Iris classification. Classical input features are encoded into Gaussian optical states, and learning proceeds through measurement contrasts that drive pseudoinverse (Gauss–Newton–like) parameter updates. Circuit architecture and component operations. The input feature vector x ∈ R4 is encoded into a four-mode photonic system using displacement gates D(xi ) applied to the vacuum state. These displacements shift the expectation values of the position quadratures, embedding classical data into photonic phase space. The displaced modes then pass through parameterized Gaussian operations consisting of rotation gates R(ψi ) and a network of beam splitter gates BSij (θk , ϕk ) that mix quadratures and couple modes. Finally, a set of beam splitters concentrates the information into a single mode, and the expectation value of its position quadrature ⟨x̂0 ⟩ is measured via x-quadrature mean measurements. See Fig. 4. Let us more specifically analyze how the sequence of Gaussian operations transforms ⟨x̂0 ⟩. 1. Initial State: The vacuum state has ⟨x̂i ⟩ = ⟨p̂i ⟩ = 0, or equivalently αi ≡ ⟨âi ⟩ = 0, for all modes i ∈ {0, 1, 2, 3}. 2. Displacement Operation: The amplitude displacement gate D (x ∈ R) = i i √
exp −i 2 xi Pi , applied to each mode i, shifts the initial amplitude αi : αi → αi + xi

(25)

so that after displacement, αi becomes 0 + xi . i ∆θphase = αi · (⟨X̂i ⟩nudged − ⟨X̂i ⟩free ),

i ∆θsq = (Var(X̂i )nudged − Var(X̂i )free ),

(19) (20)

ij ∆θbs = (Cov(X̂i , X̂j )nudged − Cov(X̂i , X̂j )free ). (21)

Similarly, a Mode B (target-referenced) contrast can be written generally as ∆θ = M (ψtarget , ψnudged ) with the following possible readouts 



i ∆θphase = αi ⟨X̂i ⟩nudged − ⟨X̂i ⟩target , (22)   i ∆θsq = Var(X̂i )nudged − Var(X̂i )target , (23)   ij ∆θbs = Cov(X̂i , X̂j )nudged − Cov(X̂i , X̂j )target .

(24)

3. Rotation Operation: The rotation gate Ri (ψi ) = exp(iψi n̂i ) acts on the complex amplitude as a phase shift: αi → eiψi αi = eiψi xi .

(26)

Collectively across all four modes, this defines the diagonal unitary UR = diag(eiψ0 , eiψ1 , eiψ2 , eiψ3 ) ∈ U (4)

(27)

so that α ⃗ → UR x. The rotation angles {ψi }3i=0 therefore control the complex phase of each mode’s contribution to the output. 4. Beam Splitter Network: The beam splitter BSij (θk , ϕk ) acts as a U (2) rotation on the pair

10 |0⟩0

D(x1 )

x0

R(ψ1 ) UL1 : BS01 (θ1 , ϕ1 )

|0⟩1

D(x2 )

Uf : BS01 π4 , 0

UL2 : BS02 (θ3 , ϕ3 )

R(ψ2 ) Uf : BS03 π4 , 0

UL2 : BS13 (θ4 , ϕ4 ) |0⟩2

D(x3 )

R(ψ3 )

|0⟩3

D(x4 )

R(ψ4 )

UL1 : BS23 (θ2 , ϕ2 )

FIG. 4. Four-mode Gaussian photonic circuit for linear classification. Each mode is initialized in vacuum |0⟩, displaced by D(xi ), and rotated by R(ϕi ), giving UR = diag(eiϕ1 , . . . , eiϕ4 ). This can be implemented by phase modulating lasers of tunable power. The first beam splitter layer UL1 = BS01 (θ1 , ϕ1 ) ⊕ BS23 (θ2 , ϕ2 ) and second layer UL2 = BS02 (θ3 , ϕ3 ) ⊕ BS13 (θ4 , ϕ4 ) mix the modes using 8 trainable parameters. The fixed concentration layer Uf = BS √ 0) concentrates information into mode 0. x-quadrature mean measurements yields P03 (π/4, 0) · BS01 (π/4, ⟨x̂0 ⟩ = i Wi xi , where Wi = 2 Re[(UBS )0i eiϕi ] and UBS = Uf UL2 UL1 .



(αi , αj ):     αi cos θk −eiϕk sin θk αi → . (28) αj αj e−iϕk sin θk cos θk

Its extension to the full four-mode space is:   gij (θk , ϕk ) = BSij (θk , ϕk ) 0 ∈ U (4), BS 0 I2 (i,j,k,l) (29) where the subscript (i, j, k, l) indicates the row/column ordering with {k, l} = {0, 1, 2, 3}\ {i, j} being the unchanged modes. The layer unitaries are then products of these extended operators: g01 (θ1 , ϕ1 ) · BS g23 (θ2 , ϕ2 ) ∈ U (4) UL1 = BS

(30)

g02 (θ3 , ϕ3 ) · BS g13 (θ4 , ϕ4 ) ∈ U (4) UL2 = BS

g01 π , 0 ∈ U (4) g03 π , 0 · BS Uf = BS 4 4

(31) (32)

The full circuit applies the following sequence of operations to the initial state α ⃗ = x (see Fig. 4): α ⃗ out = Uf UL2 UL1 UR x ≡ UBS UR x,

(33)

where UBS ≡ Uf UL2 UL1 ∈ U (4). Reading off the zeroth component, with e0 = (1, 0, 0, 0)⊤ : 3 X

α0out = e⊤ U U x = UBS 0i eiϕi xi . BS R 0

(34)

i=0

The x-quadrature measurement of mode 0 then yields: ⟨x̂0 ⟩ = ≡

√

3
√ X 
 2 Re α0out = 2 Re UBS 0i eiϕi xi i=0

3 X i=0

Wi xi ,

(35)

where the effective classifier weights are: √ 
 Wi = 2 Re UBS 0i eiϕi .

(36)

The 12 trainable parameters {ϕi }3i=0 ∪ {θk , ϕk }4k=1 jointly determine the effective linear classifier weights {Wi } through the composite map UBS UR . Consequently, the circuit output can be viewed as an inner product: ⟨x̂0 ⟩ = W · x =

3 X

Wi xi

(37)

i=0

where W = (W0 , W1 , W2 , W3 )⊤ ∈ R4 is the effective weight vector and x ∈ R4 is the input feature vector. Despite operating on photonic states, the measured output is thus a linear function of the classical inputs, with the photonic circuit parameterizing the weight vector W through the unitary degrees of freedom. Nonlinearity via Jacobian pseudoinverse updates. Although the instantaneous circuit output is linear in the inputs, ⟨x̂0 ⟩ = W(θ) · x, the parameter update rule introduces an effective nonlinearity. The Jacobian matrix J ∈ RN ×12 , stacked over a batch of N samples, has entries: Jki =

∂⟨x̂0 ⟩(k) ∂W(θ) (k) = · x ≡ gi (θ) · x(k) , ∂θi ∂θi (38)

where gi (θ) = ∂W/∂θi is generally nonlinear in θ through the trigonometric dependence of UBS . Moreover, there are other sources of nonlinearity. First, if the true input-output relationship is nonlinear, the stacked system J ∆θ = e is overdetermined and inconsistent — no single linear weight vector W can satisfy all N equations simultaneously. The pseudoinverse solution then finds a least-squares compromise across the batch, whose effective behavior depends nonlinearly on the data distribution

11 {(x(k) , e(k) )}. Second, given that the update rule for Mode B contains J ⊤ J, the entries ⊤

(J J)ij =

N X k=1

Using this notation, the pseudoinverse update becomes: ∆θ = (J T J + λI)−1 J T e

gi · x

(k)

gj · x

(k)

(39)

T

= ((XP ) (XP ) + λI) T

T

= (P X XP + λI)

−1

T

(XP ) e

(43)

T

P X e

(44)

T

⊤

(k)

are quadratic in x . The matrix inverse (J J + λI)−1 then promotes this to a rational function of {x(k) }, analogous to a polynomial kernel K(x, x′ ) = (x·x′ )2 . Therefore, although Jki is linear in x(k) , the multiplication of two such linear terms in J ⊤ J that breaks linearity of the update rule, while the minibatch overdetermination additionally encodes the nonlinear structure of the data distribution. In sum, the learned input-output map remains linear in the input features once training is complete. The nonlinearity enters the learning dynamics rather than the inference map. Mini-batch implementation. Training proceeds using mini-batches of input samples, each encoded into a separate preparation of the photonic circuit. For a batch of size N , the Jacobian becomes an N × q matrix, and the pseudoinverse update integrates response information across multiple inputs simultaneously. This improves numerical conditioning of the update and reduces the influence of measurement noise, while remaining compatible with sequential experimental execution in hardware implementations. In our implementation, a mini-batch of size 5 corresponds to evaluating five independent four-mode Gaussian circuits. Each input x ∈ R4 is encoded via displacement gates across four modes. Denoting (+) zkj = ⟨x̂0 ⟩k (θj + ε) and zkj = ⟨x̂0 ⟩k (θj ) as the perturbed and unperturbed x-quadrature mean outputs for the k-th sample and j-th parameter, the Jacobian entries can be directly expressed using the bilinear form: (+)

Jkj =

(42) −1

(+)

M (zkj , zkj ) zkj − zkj ∂⟨x̂0 ⟩k ≃ = ∂θj (θj + ε) − θj ε =

3 X i=0

(40)

Pij · xki

where xki is the i-th feature of the k-th sample. This makes clear that P emerges implicitly from the physical measurements without ever being explicitly constructed. For a batch of N samples, we can write the full Jacobian matrix as: J = XP

(41)

where X ∈ RN ×4 is the input data matrix and P ∈ R4×q contains the fixed coefficients Pij .

where e = z∗ − z is the error vector between the target and x-quadrature mean outputs, as defined previously. During training, the parameter vector θ traces a trajectory through the parameter space Rq . This trajectory is guided by: θt+1 = θt + η∆θt .

(45)

The learning rate η controls the step size along this trajectory. Classification from continuous measurements. Because the circuit produces a continuous scalar output ⟨x̂0 ⟩, classification is performed by mapping this value into discrete labels using two adaptive thresholds t1 < t2 . At each training step, the thresholds are updated as the midpoints between the mean outputs of adjacent classes, tk ← (1 − γ) tk + γ

µk + µk+1 , 2

k = 1, 2,

where µk denotes the mean measured output for class k within the current mini-batch and γ > 0 is a damping factor. The decision rule is then   0, ⟨x̂0 ⟩ < t1 , class = 1, t1 ≤ ⟨x̂0 ⟩ < t2 ,  2, ⟨x̂ ⟩ ≥ t . 0 2 This adaptive scheme ensures the thresholds track shifts in the circuit’s learned response while remaining grounded in physically measurable observables. Performance. With a 50%/50% train/test split over 25 independent trials, the test accuracy ranges from 96.0% to 100.0%, with a mean of 97.7% ± 0.9% (20 epochs). See Fig. 5. Notably, this performance is achieved using only Gaussian operations and x-quadrature measurements, without invoking intrinsic photonic nonlinearities. The classifier’s expressivity arises from the geometry of the circuit’s input–output map and the conditioning of its Jacobian under parameter perturbations. The pseudoinverse update succeeds because its physically realizable response manifold supports well-conditioned Gauss–Newton–like corrections. This demonstrates that PCPL learning performance depends primarily on how a physical system structures parameter sensitivities, rather than on the presence of specific computational primitives. More generally,

12 PCPL is compatible with nonlinear circuitry, including squeezing, and, in principle, other quantum measurement schemes, including photon-number detection, so long as measurable contrasts between nearby responses can be used to generate update signals. Changing circuit architecture. To understand what aspects of the circuit enable effective Mode B learning, we systematically altered the circuit architecture and examined how these changes affected learning performance. First, removing the final mode-concentration stage did not significantly affect accuracy. This indicates that compressing information into a single mode is not essential for expressivity in this task; the measured observable remains a linear function of the encoded inputs regardless of where mode mixing occurs. In contrast, relocating the rotation gates to the end of the circuit reduced accuracy to approximately 92.0% (base accuracy 97.33%). This modification alters how parameter perturbations influence the measured quadrature, effectively distorting the sensitivity directions encoded in the Jacobian. Introducing squeezing operations at the beginning of the circuit further degraded performance to approximately 81.3%. Squeezing introduces nonlinearities and so reshapes the circuit’s response manifold and can lead to a poorly conditioned Jacobian, making the local inverse problem harder to solve using linearized pseudoinverse updates. Similarly, fixing the beam splitter parameters reduced accuracy to approximately 78.7%, demonstrating that modemixing flexibility is crucial for spanning independent sensitivity directions in parameter space. Reintroducing trainable displacement scaling partially restored performance to approximately 92.0%, showing that increasing the space of modifiable parameters can recover controllability of the output manifold. In sum, these architectural perturbations demonstrate that learning success depends primarily on how circuit structure shapes the conditioning and expressivity of the Jacobian. Operations that preserve a well-conditioned mapping between parameter perturbations and measurable outputs enable stable Gauss–Newton–like corrections, whereas operations that distort or restrict this mapping degrade learning. This underscores that PCPL performance is governed by response geometry rather than circuit complexity alone.

2.

Tunable linear optical multiplier: Towards an Autonomous Physical Learner

Mode B PCPL requires a local update signal proportional to the product of a sensitivity and an error.

(a) Best Classification (100.0%)

(b) Batch size N = 5, η = 0.1 FIG. 5. Classification accuracy via mini-batch pseudoinverse nudging. (a) distribution of photonic measurement outcomes by class showing clear separation between Class 0 (blue, centered at -1), Class 1 (orange, centered at 0), and Class 2 (green), achieving 100.0% classification accuracy. (b) Training convergence using mini-batch pseudoinverse nudging with batch size N = 5 and learning rate η = 0.1, reaching nearly perfect accuracy within 10 epochs.

In classical gradient descent learning, this appears as  
∂f ztarget − z , (46) ∆θ ∝ ∂θ so any fully physical implementation of PCPL must be able to generate such multiplicative signals from measurable contrasts. Here we address the complementary requirement: whether the multiplication of that sensitivity by an error signal can itself be realized within a physical optical system. The mechanism we exploit—interference between displacement-encoded optical modes followed by quadrature measurement—relies only on linear optical elements and coherent states. As a result, an equivalent multiplier could in principle be implemented using purely classical electromagnetic fields. Interference of classical optical amplitudes and square-law detection have long been used in analog optical processors to perform multiplicative and convolution-like operations. The present construction does not claim novelty at the level of op-

13 tical physics; rather, its significance lies in embedding such interference-based multiplication directly into the structure of a PCPL learning rule. In this sense, the circuit functions as a physical gradient unit, translating contrast measurements into update signals within the same substrate that performs inference. This represents a step toward partially autonomous physical learners, in which not only inference but also elements of the update rule are implemented directly in the physical substrate. We define a physical learning system as autonomous to the extent that components of the parameter-update rule are implemented within the same physical substrate that performs inference, rather than being computed by an external digital processor. Goal. Our goal is to construct a linear optical circuit whose measured output approximates the analog product δf · e, thereby embedding a fundamental part of the learning operation into the measurement process. We consider two approaches within linear optics. In the first approach, the signals δf and e are encoded into quadrature amplitudes via displacement gates and the encoded modes are coupled through a beam splitter. A linear fit is then applied to x-quadrature mean measurements outcomes to find parameters that minimize the approximation error. Note this gives you δf ± e, ie, δf and e separately, not their product. In the second approach, we exploit the intensitydifference property of a beam splitter. Encoding δf and e into the two input modes, the beam √ splitter produces output amplitudes α = (δf − e)/ 2 and 0 √ α1 = (δf + e)/ 2. The difference of squared direct detection then yields α2 − α02 δf ·e= 1 ≈ gk , ε 2·ε

(47)

where gk is the gradient estimate weighted by the error, obtained from intensity measurements or classical post-processing of x-quadrature outputs, scaled by the fixed perturbation size ε, so the product (δf /ε) · e can be recovered directly from the intensity measurement without requiring any nonlinear optical elements. Both approaches yield comparable final classification accuracy. Implementation. The parameter update is implemented by a second CV circuit in which each parameter mode is initialized via a displacement D(θk ) and the scaled gradient is added by a second displacement D(ηgk ), so that the updated value is read from the quadrature mean ⟨x̂k ⟩/2 = θk + ηgk . This minimal circuit thus represents a CV self-learning system in which gradient computation and parameter update all occur within integrated Gaussian circuits, with the exception of classical finite-difference perturbations used to evaluate δf and parameter clip-

Classical finite difference δfk

= e

f (x; θk + ε) − f (x; θ)

=

ytrue − f (x; θ) δfk , e

BS multiplier circuit (2 modes; fixed 50/50 BS; runs 12× per sample) (k)

D(δfk , 0)

|0⟩0

BS(π/4, 0)

(k)

⟨x̂1 ⟩(k)

D(e, 0)

|0⟩1

⟨x̂0 ⟩(k)

⟨x̂0 ⟩(k) , ⟨x̂1 ⟩(k)

Classical post-processing gk

=

α12 − α02 1 · 2 ε

where αi

=

⟨x̂i ⟩ 2

g1 , . . . , g12 D gate parameter update (12 modes; representative mode k) |0⟩k

D(θk , 0)

D(ηgk , 0)

⟨x̂k ⟩

θk ← ⟨x̂k ⟩/2 = θk + ηgk

FIG. 6. Simplest linear optical learning circuit for Iris dataset classification. The circuit consists of a classical finite-difference block, a two-mode CV BS multiplier circuit for gradient computation (run once per parameter per sample), a classical post-processing block for intensity-difference readout, and a 12-mode D gate parameter update circuit with x-quadrature mean readout.

ping applied to the updates. See Fig. 6. We now incorporate the linear optical multiplier into the larger linear optical learning circuit. To classify the Iris dataset, the system uses a 4-mode photonic classifier with 12 trainable parameters (4 rotation angles and 8 beam splitter parameters). The classifier encodes iris features through displacement gates D̂(xi ), applies learnable rotations R̂(θi ), and processes the state through two layers of parameterized beam splitters followed by an additional adaptive rotation layer and two final beam splitter operations. The forward pass produces a continuous x-quadrature mean measurement output. Performance. The minimal two-mode linear optical multiplier performs parameter updates through a dedicated two-mode displacement circuit: each parameter mode is initialized by D̂(θk ) and the scaled gradient is added by D̂(ηgk ), so the updated value is recovered as ⟨x̂k ⟩/2 = θk + ηgk . No additional processing stages are involved. The system achieves 96.2% ± 2.3% (Configuration A, 20 epochs)/training and testing accuracy for the Iris dataset (25 independent realizations), demonstrating that photonic

14 physical self-learning with photonic-computed gradients is feasible. See Appendix C for an example with squeezing as well as the effect of circuit architecture on the performance.

IV.

DISCUSSION

This work introduces Perturbative Contrastive Physical Learning (PCPL) as a unifying principle for learning in physical systems. Rather than defining learning through explicit gradients of a loss function, PCPL frames learning as the measurable contrast between nearby physical states produced by controlled perturbations. In this view, the fundamental primitive of learning is not symbolic differentiation but experimentally accessible response. By elevating contrast to the central role, PCPL provides a substrate-independent description of learning that applies equally to mechanical, optical, and photonic systems. This reframing shifts the question from “How does one compute gradients?” to “What contrasts can a physical system sense and translate into parameter change?” A key insight of the PCPL framework is the distinction between two learning geometries. Mode A learning probes local sensitivity: the system compares its response before and after a small internal perturbation, effectively measuring a Jacobian action. Mode B learning performs local inversion: contrasts between nudged responses and targets induce Gauss–Newton–like parameter updates governed by the natural metric of the system’s input–output map. Importantly, this inverse geometry need not be computed explicitly. In the physical implementations shown here, the effective metric emerges implicitly from the system’s response, demonstrating that metric-aware learning can be realized through local perturbation and measurement rather than global matrix inversion. In the two platforms, a spring network and a CV photonic gates based circuit, learning is successful, suggesting that learning performance in physical substrates depends less on whether a system is classical or photonic, and more on how its architecture shapes the conditioning and expressivity of its response manifold. The tunable continuous-variable photonic multiplier extends this idea by embedding part of the learning rule itself as part of the physical substrate. By implementing an analog multiplication between a sensitivity signal and an error signal through optical interference, we show that a core computational primitive of Gauss–Newton learning can migrate from classical post-processing into the physical substrate. This represents an initial step toward autonomous physical learners in which inference and

elements of parameter update emerge from the same substrate. In such systems, learning is no longer an algorithm imposed from outside, but an algorithmic process shaped by the system’s own response geometry. Previous work has demonstrated that physical systems can be trained through response-based adaptation, including the programming of stress and strain patterns in spring networks and cell packings in disordered materials [30, 31]. These approaches emphasize experimentally realizable update rules, often with a strong focus on strictly local adaptation. PCPL does not adopt locality as a defining principle. Instead, it is organized around the geometry of physically measurable contrasts and the contrast measurement functional M need not act on strictly local observables. In fact, the measurable contrasts may involve collective, nonlocal, interferometric, or globally encoded response variables, provided that they contain sufficient information about the underlying response geometry of the system. In Mode A, these contrasts probe local sensitivities, while in Mode B they can encode information about a global inverse problem through the system’s collective response. Importantly, this does not imply centralized learning: no external processor computes gradients or inverts system-wide Jacobians. Instead, the learning signals arise directly from the system’s own physical response to controlled perturbations. While the resulting updates can reflect global properties of the system, such as collective modes, longrange couplings, or inverse-problem geometry, these effects are mediated through the substrate’s intrinsic response rather than through an explicit global optimization algorithm. In other words, the system does not require a separate computational layer with access to a full model of its state; the effective learning geometry emerges implicitly from measurable contrasts between nearby physical configurations. In PCPL, by contrast, the “computation” of update directions is distributed across the physical degrees of freedom themselves, and is realized through experimentally accessible perturbation–response measurements rather than algorithmic inversion. PCPL therefore occupies an intermediate regime between strictly local adaptation and centrally computed gradient descent, demonstrating that nonlocal learning behavior can arise without centralized algorithmic control. As for experimental implications, the optical character of CV circuit brings about the distinct possibility of physical learning with low energy consumption compared to Joule-heating electronic circuits, even if one must take into account the amount of laser power necessary to drive the circuit and, specifically, to pump the optical parametric oscilla-

15 tors that produce squeezed light. Nevertheless, perspectives are enticing, especially in light of the constantly evolving state of the art of integrated photonics [32], as integration begets scalability. Learning with continuous quantum variables using Gaussian and non-Gaussian gates [33] has been studied extensively, including Continuous-Variable Quantum Neural Networks [34] and relevant applications to state classification (with classical feedforward postprocessing) [35] and time-series forecasting (with classical optimization) [36]. In addition, recent work has demonstrated that learning can be implemented in linear optical interferometers with coherent states and homodyne measurements [37]. These studies demonstrate that CV quantum circuits can be used for learning tasks; however, classical computational processing blocks are still required. In other words, to realize a fully autonomous quantum learning system, aka a quantum creature, perturbations, contrasts using single-shot measurements, multiplication/error couplings, and parameter updates are all implemented within the quantum substrate itself.

We are currently working towards this goal. More broadly, PCPL expands the notion of what constitutes a learnable parameter. In physical systems, modifiable degrees of freedom include not only coupling strengths but also geometry, boundary conditions, and even interpretive mappings between physical states and functional outputs. By defining learnability operationally—through what can be perturbed and sensed—PCPL shifts the focus from abstract optimization to experimentally accessible contrasts between physical states. This perspective opens pathways for designing adaptive materials, programmable photonic circuits, hybrid classical–quantum, and fully quantum, platforms, where learning is embedded directly in physical laws. Acknowledgements. We thank Benjamin Scellier for useful discussions. JMS acknowledges financial support from the National Science Foundation via DMR-2204312. OP was supported by National Science Foundation grants OSI-2531569 [NQVL:QSTD:Design: Quantum Computing Applications of Photonics (QCAP)], PHY-2514971, and ECCS-2530171.

[1] Donald O. Hebb. The Organization of Behavior: A Neuropsychological Theory. Wiley, New York, 1949. [2] Peter Dayan and Larry F. Abbott. Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems. MIT Press, Cambridge, MA, 2001. [3] Chris R. Reid, Tanya Latty, Audrey Dussutour, and Madeleine Beekman. Slime mold uses an externalized spatial memory to navigate in complex environments. Proceedings of the National Academy of Sciences, 113(47):E7051–E7058, 2016. [4] Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, 2016. [5] Benjamin Scellier and Yoshua Bengio. Equilibrium propagation: Bridging the gap between energybased models and backpropagation. Frontiers in Computational Neuroscience, Volume 11 - 2017, 2017. [6] Logan G. Wright, Tatsuhiro Onodera, Martin M. Stein, Tianyu Wang, Darren T. Schachter, Zoey Hu, and Peter L. McMahon. Deep physical neural networks trained with backpropagation. Nature, 601(7894):549–555, 2022. [7] Menachem Stern and Arvind Murugan. Learning without neurons in physical systems. Annual Review of Condensed Matter Physics, 14(Volume 14, 2023):417–441, 2023. [8] Ali Momeni, Babak Rahmani, Benjamin Scellier, Logan G. Wright, Peter L. McMahon, Clara C. Wanjura, Yuhang Li, Anas Skalli, Natalia G. Berloff, Tatsuhiro Onodera, et al. Training of physical neural networks. Nature, 645(8079):53–61, 2025.

[9] Jérémie Laydevant, Logan G Wright, Tianyu Wang, and Peter L McMahon. The hardware is the software. Neuron, 112(2):180–183, 2024. [10] Benjamin Scellier. A deep learning theory for neural networks grounded in physics. arXiv preprint arXiv:2103.09985, 2021. [11] Menachem Stern, Daniel Hexner, Jason W. Rocks, and Andrea J. Liu. Supervised learning in physical networks: From machine learning to learning machines. Physical Review X, 11(2):021045, 2021. [12] Menachem Stern, Marcelo Guzman, Felipe Martins, Andrea J. Liu, and Vijay Balasubramanian. Physical networks become what they learn. Physical Review Letters, 134(14):147402, 2025. [13] Sam Dillavou, Menachem Stern, Andrea J. Liu, and Douglas J. Durian. Demonstration of decentralized physics-driven learning. Physical Review Applied, 18(1):014040, 2022. [14] Lauren E. Altman, Menachem Stern, Andrea J. Liu, and Douglas J. Durian. Experimental demonstration of coupled learning in elastic networks. Phys. Rev. Appl., 22:024053, Aug 2024. [15] Sam Dillavou, Benjamin D. Beyer, Menachem Stern, Andrea J. Liu, Marc Z. Miskin, and Douglas J. Durian. Emergent learning in a nonlinear analog network. Proceedings of the National Academy of Sciences, 2024. [16] Vidyesh Rao Anisetti, Benjamin Scellier, and J. M. Schwarz. Learning by non-interfering feedback chemical signaling in physical networks. Physical Review Research, 5:023024, 2023. [17] Vidyesh Rao Anisetti, Ananth Kandala, Benjamin

16 Scellier, and J. M. Schwarz. Frequency propagation: Multimechanism learning in nonlinear physical networks. Neural Computation, 36(4):596–620, 2024. [18] Martin J. Falk, Adam T. Strupp, Benjamin Scellier, and Arvind Murugan. Temporal contrastive learning through implicit non-equilibrium memory. Nature Communications, 16(1):2163, 2025. [19] Benjamin Scellier. Quantum equilibrium propagation: Gradient-descent training of quantum systems. arXiv preprint arXiv:2406.00879, 2024. [20] Clara C. Wanjura and Florian Marquardt. Quantum equilibrium propagation for efficient training of quantum systems based on onsager reciprocity. Nature Communications, 16(1):6595, 2025. [21] Mile Gu, Christian Weedbrook, Nicolas C. Menicucci, Timothy C. Ralph, and Peter van Loock. Quantum computing with continuous-variable clusters. Physical Review A, 79:062318, 2009. [22] Nicolas C. Menicucci, Steven T. Flammia, and Olivier Pfister. One-way quantum computing in the optical frequency comb. Physical Review Letters, 101:130501, 2008. [23] Olivier Pfister. Continuous-variable quantum computing in the quantum optical frequency comb. Journal of Physics B: Atomic, Molecular and Optical Physics, 53(1):012001, 2020. [24] Nathan Killoran, Thomas R. Bromley, Juan Miguel Arrazola, Maria Schuld, Nicolás Quesada, and Seth Lloyd. Continuous-variable quantum neural networks. Physical Review Research, 1:033063, 2019. [25] Maria Schuld, Ville Bergholm, Christian Gogolin, Josh Izaac, and Nathan Killoran. Evaluating analytic gradients on quantum hardware. Physical Review A, 99:032331, 2019. [26] Dong-Hyun Lee, Saizheng Zhang, Asja Fischer, and Yoshua Bengio. Difference target propagation. In Joint european conference on machine learning and knowledge discovery in databases, pages 498–515. Springer, 2015. [27] Shuailong Li and Xinqiao Mao. Training allmechanical neural networks for task learning through in situ backpropagation. Nature Communications, 15(1):10528, 2024. [28] Heru Nugroho, Nugraha Priya Utama, and Kridanto Surendro. Performance evaluation for class centerbased missing data imputation algorithm. In Proceedings of the 2020 9th International Conference on Software and Computer Applications, ICSCA ’20, page 36–40, New York, NY, USA, 2020. Association for Computing Machinery. [29] Benjamin Scellier. Quantum equilibrium propagation: Gradient-descent training of quantum systems. arXiv:2406.00879 [quant-ph], 2024. [30] Daniel Hexner. Training precise stress patterns. Soft Matter, 19:2120–2126, 2023. [31] Shabeeb Ameen, Tao Zhang, and J. M. Schwarz. Training cell stress patterns in 3d cellular packings. arXiv preprint arXiv:2604.25439, 2026. [32] Wim Bogaerts, Daniel Pérez, José Capmany, David A. B. Miller, Joyce Poon, Dirk Englund, Francesco Morichetti, and Andrea Melloni. Programmable

photonic circuits. Nature, 586(7828):207–216, 2020. [33] S. L. Braunstein and P. van Loock. Quantum information with continuous variables. Rev. Mod. Phys., 77(2):513, 2005. [34] Nathan Killoran, Thomas R. Bromley, Juan Miguel Arrazola, Maria Schuld, Nicolás Quesada, and Seth Lloyd. Continuous-variable quantum neural networks. Phys. Rev. Res., 1:033063, Oct 2019. [35] Mohammad Rezaei Shokou and Hossein Davoodi Yeganeh. Hybrid quantum-classical learning of nonlinear entanglement witnesses via continuousvariable quantum neural networks, 2025. [36] Prabhat Anand, M Girish Chandra, and Ankit Khandelwal. Time-series forecasting using continuous variables-based quantum neural networks. In 2024 16th International Conference on COMmunication Systems & NETworkS (COMSNETS), pages 994–999, 2024. [37] Amanuel Anteneh, Kyungeun Kim, J. M. Schwarz, Israel Klich, and Olivier Pfister. Laser interferometry as a robust neuromorphic platform for machine learning. arXiv preprint arXiv:2601.18047, 2026.

APPENDIX A: MODE B GRADIENT DESCENT VS. PSEUDOINVERSE PERFORMANCE COMPARISON

We compare gradient descent (GD) and the pseudoinverse update rule on a linear classifier trained on the Iris dataset (50%/50% train–test split, 25 independent trials). Two experiments are conducted to isolate the effect of the update rule (Fig. 7). First, we sweep the learning rate η over a wide range to test sensitivity: GD requires careful tuning of η, whereas the pseudoinverse update is self-normalizing via (J ⊤ J)−1 and therefore robust to the choice of η. Second, we apply uneven feature scaling to introduce ill-conditioning: one feature is amplified by a factor of 102 (or 103 ) while another is suppressed by 10−2 (or 10−3 ), creating a highly anisotropic loss landscape in which GD oscillates or diverges while the pseudoinverse update converges stably. In both experiments, adaptive thresholds and parameter clipping are deliberately excluded to isolate the effect of the update rule alone; these techniques can improve accuracy by up to 4% in practice. Both methods use a 1-output linear model with fixed thresholds. Table I demonstrates this on a linear classifier trained on the Iris dataset: while GD and the pseudoinverse update achieve comparable accuracy at a well-tuned learning rate (η = 0.05–0.1), GD degrades rapidly for larger η or under feature ill-conditioning, whereas the pseudoinverse update maintains consistent performance throughout. The numerical results in Table I confirm that the pseudoinverse update of Mode B is robust to hyperparameter choice and input ill-conditioning, a conse-

17 Method η = 0.05 η = 1.0 scaling (102 , 10−2 ) scaling (103 , 10−3 ) GD 93.7% ± 1.4% 65.9% ± 2.7% 60.6% ± 1.1% 45.3% ± 13.5% Pseudoinverse 94.1% ± 2.9% 94.2% ± 2.2% 94.0% ± 1.6% 93.6% ± 1.8% TABLE I. Test accuracy for Iris classification(50%/50% split, 25 trials) comparing GD and pseudoinverse update under varying learning rates and feature ill-conditioning. Note that additional techniques such as adaptive thresholds and parameter clipping were deliberately excluded to isolate the effect of the update rule; these can improve accuracy by up to 4% in practice.

quence of its metric-aware correction in the natural geometry of the system’s response manifold. APPENDIX B: MORE DETAILS ON PCPL Least Squares Problem with Tikhonov regularization

(48)

where J ∈ RN ×p is the Jacobian of the output observable with respect to the system parameters. We seek the parameter update ∆θ that minimizes the regularized objective min J ∆θ − e ∆θ

2

+ λ ∥∆θ∥2 ,

G(y, θ ′ ) − G(y, θ) ≈ Jθ G(y, θ) · (θ ′ − θ),

(49)

where the second term λ ∥∆θ∥2 penalizes large parameter updates and λ > 0 controls the strength of regularization. Taking the derivative with respect to ∆θ and setting it to zero yields i ∂ h 2 J ∆θ − e + λ ∥∆θ∥2 ∂ ∆θ = 2J ⊤ (J ∆θ − e) + 2λ ∆θ = 0, (50) which gives (J ⊤ J + λI) ∆θ = J ⊤ e. The closed-form solution is thus the regularized pseudoinverse step
−1 ⊤ ∆θ = J ⊤ J + λI J e, λ>0 (51) to arrive at Eq. 10. Linear Response and the Contrast Measurement Functional M

To obtain a direct mapping from parameter perturbations and M , we start with the approximation X ∂G G(y, θ ′ ) − G(y, θ) ≈ (y, θ) · (θi′ − θi ). (52) ∂θ i i

(53)

where Jθ G(y, θ) is the Jacobian matrix with elements: [Jθ G(y, θ)]j,i =

For N input samples, we collect the measured quadrature outputs into a vector z ∈ RN and define z∗ as the desired class-dependent targets. Assuming linear response of the system with respect to parameters θ ∈ Rp , the residual e = z∗ − z satisfies e ≈ J ∆θ,

The RHS can be written more compactly using the Jacobian matrix Jθ G of partial derivatives:

∂Gj (y, θ), ∂θi

(54)

and Gj is the j-th component of the output vector G. Similarly, the contrast measurement functional M applied to the output contrast can be approximated using its own Jacobian: M (G(y, θ ′ ), G(y, θ)) ≈ JG M (G, G′ ) · (G(y, θ ′ ) − G(y, θ)),

(55)

where G′ = G(y, θ ′ ) and JG M (G, G′ ) is the Jacobian of the measurement operator with respect to changes in G. Combining these two approximations, we can derive the relationship between parameter perturbations and measurement outputs: M (G(y, θ ′ ), G(y, θ)) ≈ JG M (G, G′ ) · (G(y, θ ′ ) − G(y, θ)), ≈ JG M (G, G′ ) · Jθ G(y, θ) · (θ ′ − θ).

(56)

This gives us a direct mapping from parameter perturbations to measurement outputs through the composition of two Jacobian transformations. Gradient Descent Approximation and Beyond

We start with contrast measurement functionals that satisfy: JG M (G, G′ ) · Jθ G(y, θ) ≈ γ · ∇θ L(θ)T ,

(57)

for some positive constant γ and loss function L(θ). One way to achieve this is to design the measurement operator such that: JG M (G, G′ ) ≈ λ · (Jθ G(y, θ))+ ,

(58)

where (Jθ G)+ is the pseudoinverse of the Jacobian and λ is a scaling factor. This design choice enables

18

FIG. 7. Robustness comparison between gradient descent (GD) and pseudoinverse parameter updates for Iris classification. (a) Test accuracy as a function of learning rate η. GD maintains competitive performance for small η but degrades sharply above η = 0.5, whereas the pseudoinverse update remains stable across the full range tested. (b) Test accuracy under increasingly ill-conditioned feature scaling (s0 , s3 ), where s0 and s3 denote the scaling factors applied to the first and last features, respectively. GD accuracy drops substantially as the condition number increases, while the pseudoinverse update retains above 93% accuracy throughout. Error bars indicate one standard deviation over 25 independent trials.

us to construct a parameter update rule that approximates gradient descent. With this design, the parameter update becomes: θnew = θ + η · M (G(y, θ ′ ), G(y, θ))

(59)

(i)

≈ θ + η · JG M (G, G′ ) · Jθ G(y, θ) · (θ ′ − θ) (60)

(ii)

≈ θ + ηλ · (Jθ G)+ · Jθ G · (θ ′ − θ)

(61)

= θ + ηλ · (θ ′ − θ),

(62)

(iii)

where (i) applies the first-order linear approximation, (ii) substitutes the design from equation (2), and (iii) uses the pseudoinverse property (Jθ G)+ · Jθ G ≈ I. For nudging rules where θ ′ − θ ∝ ∇θ L(θ), this results in effective gradient descent toward minimizing the loss function. However, more generally one can treat G(y, θ) as a black box and measure linear response to update θ given the infinitesimal linear approximation. This approach works well for single-output problems and small multi-output cases (n ≤ 3). For higher-dimensional outputs (n > 3), on the other hand, the inverse mapping δθ ∝ G−1 (ztarget − G(y, θ)) relies heavily on infinitesimal linear approximations through Jacobianbased updates. In such cases, establishing explicit non-linear models for the intermediate y → z mappings will presumably improve training efficiency compared to relying solely on local linear approximations. The accuracy of such approximations depends on the size of the perturbation ∥θ ′ − θ∥ and the nonlin-

earity of G. The error term is: Error = G(y, θ ′ ) − G(y, θ) − Jθ G(y, θ) · (θ ′ − θ). (63) For sufficiently small perturbations, this error becomes negligible compared to the first-order term, which is the key assumption that enables effective parameter learning through PCPL. APPENDIX C: MORE ON TUNABLE MULTIPLIERS Tunable multiplier with squeezing

We consider a tunable multiplier using squeezing gates. In a single-mode Gaussian circuit, applying a displacement D̂(d · e) followed by squeezing Ŝ(s · δf ) to vacuum gives an x-quadrature mean that depends on both signals. Here d and s are trainable parameters, so d · e and s · δf are simply classical signals scaled by fixed gains before encoding into the circuit. However, we will show how to find suitable values based on the target system. The interaction between the two signals arises from the nonlinear dependence e−s·δf introduced by the squeezing gate acting on the already-displaced state. √ ⟨x̂⟩ = 2 d · e · e−s δf (64) since the squeezing scales the displacement mean by e−s δf , where rotation is not yet applied. In the regime s · δf ≪ 1 this approximates: √ ⟨x̂⟩ ≈ 2 d · e · (1 − s δf ) (65)

19 which contains the product term s · d · δf · e and the first term is fixed values for fixed e. We therefore define a single-mode circuit layer as: R̂(ϕ) Ŝ(s · δf ) D̂(d · e) |0⟩

(66)

where s, d, and ϕ are trainable parameters and the rotation R̂(ϕ) aligns the output quadrature with the √ measurement axis, giving ⟨x̂⟩ = 2 d · e · e−s δf cos ϕ. This circuit is trained to minimize 1 − R2 between its output and the classical gradient (δf /ε) · e over a dataset of realistic (δf, e) pairs generated from the CV linear classifier with random parameters. We refer to this as the DSR circuit. The full circuit with L layers has trainable parameters {sℓ , dℓ , ϕℓ }L ℓ=1 , totalling 3L parameters, with output ⟨x̂0 ⟩. Figure 8 shows the circuit diagram and the training performance. The inputs δf = f (x; θi + ε) − f (x; θ) and e = ytrue − f (x; θ) are computed classically before being encoded into the circuit. The circuit parameters {sℓ , dℓ , ϕℓ }L ℓ=1 are pre-trained offline and fixed during Iris classifier training. Specifically, the DSR multiplier is trained offline using a dataset of (δf, e, gclassical ) tuples generated by running the CV linear classifier with randomly sampled parameter vectors θ ∼ Uniform(−0.3, 0.3)12 on Iris training samples, with nθ = 50 random parameter sweeps per sample. For each sweep, a random parameter index i is selected, the finite difference δf = f (x; θi + ε) − f (x; θ) is computed with ε = 0.015, and the error e = ytrue −f (x; θ) is recorded. The circuit parameters are optimized by minimizing 1 − R2 between the circuit x-quadrature mean output and gclassical = (δf /ε) · e using a coordinate-wise random search optimizer over 400 iterations with initial learning rate 0.1 and decay factor 0.997. Inputs δf and e are normalized by their training-set standard deviations before encoding. Now we test for classification by combining the DSR multiplier (with pre-trained data) and the D gate parameter update circuit (same as the main text). Figure 9 shows the classification performance with this combined circuit. Learning is inconsistent across trials due to approximation error introduced by the nonlinear squeezing-based multiplication, which only approximates δf · e in the smallsignal regime s · δf ≪ 1. When this condition is violated, the gradient signal deviates from the true error-weighted finite difference, causing parameter updates to point in incorrect directions. Within individual trials, fluctuations in training accuracy are also observed; to mitigate this, we record the bestperforming parameters during training and restore them at the end (15 epochs, learning rate η = 1.0). Nevertheless, in successful trials the best accuracy reaches 96.0%, demonstrating that the DSR multi-

plier can provide sufficient gradient information for effective learning when the operating conditions are favorable. The inconsistency arises from two main sources: (1) parameter range clipping in the D gate update circuit, which truncates large parameter updates and prevents the classifier from reaching certain regions of parameter space; and (2) the nonlinearity of the squeezing-based multiplication, which introduces a systematic bias in the gradient estimate that grows with the magnitude of δf . Since both sources of error are data- and initializationdependent, some random seeds lead to favorable operating regimes where the approximation holds well throughout training, while others do not. This sensitivity to initialization is an inherent limitation of the pure CV approach without planned circuit design. A more carefully designed CV multiplier circuit, with explicit consideration of the operating range of δf and e and a circuit architecture that reduces the dependence on the small-signal approximation, would be expected to improve both consistency and accuracy across trials.

Tunable linear optical multiplier

: Configuration A (from main text) and classical update comparison To isolate the contribution of each CV component, we compare four configurations: • Configuration A pairs the linear optical multiplier gradient with the CV displacement update (the full CV system) • Configuration B replaces the CV gradient with a classical finite-difference gradient while retaining the CV displacement update • Configuration C retains the CV gradient but replaces the displacement update with classical gradient descent • Configuration D uses classical gradients and classical gradient descent throughout, serving as the full classical baseline. Table II summarizes the results across 25 independent trials. All four configurations achieve identical performance: 98.7% best test accuracy, 96.2% final test accuracy, and 2.3% standard deviation. The absence of any measurable difference across configurations indicates that the linear optical multiplier gradient is functionally equivalent to the classical finitedifference gradient, and that the CV displacement update is functionally equivalent to classical gradient descent, at least for this task. This is consistent

20 |0⟩

D(dℓ · e, 0)

S(sℓ · δf )

R(ϕℓ )

⟨x̂0 ⟩

(a) Trainable circuit diagram

(b) Training performance FIG. 8. DSR multiplier.

FIG. 9. Performance of the DSR multiplier across trials.

with the theoretical expectation that the linear optical multiplier circuit computes δf · e/ε to machine precision via the beam splitter intensity-difference identity, and that the displacement update circuit implements arithmetic addition exactly within the Gaussian backend. The CV components therefore

do not introduce approximation errors relative to their classical counterparts, and the full CV system of Configuration A achieves classical-equivalent accuracy while realizing both gradient computation and parameter update as physical analog operations. See Fig. 10.

21 TABLE II. Ablation results across 25 independent trials (20 epochs). Config Train (best) Train (final) Test (best) Test (final) Std (test) A 98.7% 96.3% 98.7% 96.2% 2.3% B 98.7% 96.3% 98.7% 96.2% 2.3% C 98.7% 96.3% 98.7% 96.2% 2.3% D 98.7% 96.3% 98.7% 96.2% 2.3%

FIG. 10. Performance comparison among different combinations of classical and CV learning components using the same seeds.