# Deep physical neural networks enabled by a backpropagation algorithm for
arbitrary physical systems: supplementary material

###### Contents

  1. 1 Physics-aware training
    1. 1.1 Intuition for why physics-aware training works
    2. 1.2 General formulation of physics-aware training
    3. 1.3 Motivating physics-aware training with multilayer feedforward architecture
    4. 1.4 Numerical example
  2. 2 Supplementary methods
    1. 2.1 Femtosecond second-harmonic generation PNN
      1. 2.1.1 Experimental setup
      2. 2.1.2 Input-output transformation characterization
    2. 2.2 Analog transistor PNN
      1. 2.2.1 Experimental setup
      2. 2.2.2 Input-output transformation characterization
    3. 2.3 Oscillating plate PNN
      1. 2.3.1 Experimental setup
      2. 2.3.2 Input and output encoding
      3. 2.3.3 Characterization of input-output transformation
    4. 2.4 Data-driven differentiable models
      1. 2.4.1 Digital model for the mean
      2. 2.4.2 Noise digital model
    5. 2.5 Descriptions of PNN architectures
      1. 2.5.1 Figure 2: Femtosecond SHG-vowels PNN
      2. 2.5.2 Figure 4: Analog transistor-MNIST PNN
      3. 2.5.3 Figure 4: Oscillating plate MNIST PNN
      4. 2.5.4 Figure 4: Femtosecond SHG-MNIST PNN
  3. 3 Simulation analysis of candidate PNN systems
    1. 3.1 General case
    2. 3.2 Single-mode nonlinear optical pulse propagation: scaling
      1. 3.2.1 Quantitative example
    3. 3.3 Self-simulation analysis of multimode wave propagation
      1. 3.3.1 Linear propagation
      2. 3.3.2 Nonlinear propagation
    4. 3.4 Multimode nonlinear oscillators/oscillator networks
    5. 3.5 Nonlinear electrical circuit dynamics/analog transistor networks
  4. 4 Supplementary discussion: when can PNNs provide technological benefits?

##  1 Physics-aware training

In this section, we describe physics-aware training (PAT) and explain its
effectiveness at training physical neural networks (PNNs). In section 1.2, we
develop the algorithm in full mathematical detail, comparing it to the
conventional backpropagation algorithm. In Sec. 1.3, this general case is
specified to a generic multi-layer physical neural network. In that section,
we introduce and compare physics-aware training to the most relevant
alternative for training sequences of physical systems with backpropagation,
_in silico_ training, in which the backpropagation algorithm is applied to a
simulation of the physical system(s). Finally, in Sec. 1.4, we present a
numerical example to demonstrate how PAT efficiently trains PNNs, by training
a simulated feedforward PNN to perform the vowel classification task with the
different training algorithms introduced in this section.

Before delving into more detailed subsections, however, the next subsection
provides a high-level intuition for why physics-aware training works, and why
alternatives fail. This intuition is backed up rigorously in subsequent
subsections.

###  1.1 Intuition for why physics-aware training works

The backpropagation algorithm solves the problem of efficiently training the
functionality of sequences of mathematical operations with many trainable
parameters to perform desired mathematical functions. Its key ingredient is
the analytic differentiation of the mathematical operations to efficiently
compute the gradient, and hence to determine the optimal parameter updates to
improve performance.

Physics-aware training (PAT) solves the problem of applying the
backpropagation algorithms to train sequences of real physical operations to
perform desired physical functions – which includes performing physically-
implemented computations, i.e. physical neural networks.

There are two natural alternatives to PAT. The backpropagation algorithm
cannot be applied directly – real physical systems cannot be analytically
differentiated. Instead, one could estimate derivatives by using a finite-
difference estimator. But this requires $n$ uses of the physical system, where
$n$ is the number of trainable parameters. Since $n$ is often $10^{6}$ or even
$10^{9}$, and keeps increasing as deep neural network models are scaled up,
finite-difference just isn’t a practical solution, even for very fast physical
systems. The other alternative is to train the physical system using numerical
simulations, which are ultimately made up of mathematical operations that one
can apply the backpropagation algorithm to. The problem with this approach,
however, is that a simulation is not reality – even the most finely-tuned
simulation will only describe reality to within some finite precision. The
existence of this inevitable simulation-reality gap leads to two problems.

The first problem is the build-up of the simulation-reality gap through
layers. Even if the simulation-reality gap is very small (i.e., the model is
very accurate), it can matter when we have a process that is composed of a
sequence of physical layers, where each layer after the first one takes in the
output of the previous layer as an input. Basic error propagation shows us
that a small simulation-reality gap (i.e., a small difference between reality
and the simulation of reality) is going to build up, usually very rapidly,
through layers.

This is a very simple problem, so simple that it can be easily misunderstood.
Consider a toy example. Suppose we have a function $f(x)=2x^{1.1}$ (the
‘reality’) and a very good ‘model’ for that function $f_{m}(x)=2.01x^{1.1}$.
If we take a composition of $f(x)$, $f(f(f(f(f(…f(x))))))=f^{n}(x)$, as occurs
in a deep network or general hierarchical process, the small inaccuracies of
the model blow-up. For $n=1$, the simulation-reality gap (the relative error
in the model’s output compared to the reality) is about 0.5%. For $n=5$, it
has grown to 3%. By $n=20$, the gap is 30%. Of course, this is only a toy
example: the blow-up depends, among other things, on the function, where it is
evaluated, and how and where the model is inaccurate. The key point is that
the model is increasingly wrong as the depth increases.

The second problem is even bigger – the build-up of the simulation-reality gap
through training. Training is a sequential process that involves $10^{5}$, or
even many more steps, each of which depends on the previous one. Each step of
training results in new parameters, $\theta_{i}=T(\theta_{i-1},f)$, where $T$
is the function(al) that performs one training step and $\theta_{i}$ is the
parameter vector on the $i$th training step. The final parameters at training
step $n$ are thus a composition of many applications of $T$,
$\theta_{n}=T^{n}(\theta_{0},f)$. Just as the simulation-reality gap grows
with the number of layers, so too does it grow with training steps. Since
there are usually many more training steps than layers, and because each
training step is a function of the full-depth network (so it is fully affected
by the first problem), this second problem is the most significant, and means
that even extremely accurate models that can predict deep networks well, are
still not good enough to facilitate training of real physical systems.

Physics-aware training mitigates these problems. It solves the first problem,
the depth problem, by simply using the true $f(x)$. The second problem, the
training-step problem, is mitigated by PAT both because the simulation-reality
gap for $f(x)$ is zero, and because the parameter updates are estimated more
accurately than if backpropagation were just applied directly to the model.
This second point is more subtle and will be elaborated on in rigorous detail
in subsequent sections. An approximate intuition, however, is that PAT assures
that when we adjust parameters to change the output of a sequence of physical
operations, we are accurately estimating which direction the output needs to
change to improve its true performance. In contrast, when just performing
backpropagation on a model, we may inadvertently think that the network is
overshooting the desired output, when in fact, if the simulation’s parameters
were applied to the true physical system, the output would be undershooting. A
second aspect of this improved parameter update estimation is that, even
though PAT uses a model to estimate gradients, it uses that model at the
correct locations for each layer. For example, in the toy example above, if
one is performing backpropagation on the model, the location at which
derivatives of the model are evaluated is off by only 0.5% for the second
layer (which receives the first layer’s output as an input), but is off by 30%
for the 21st layer. With PAT, the gradients are always estimated assuming the
correct inputs. The net effect is that, with PAT, the simulation-reality gap
grows much more slowly than if training were performed solely with
simulations. Since PAT is always grounded in the true loss, training can never
become misled about how well the physical computation is performing. As a
result, even if a small simulation-reality gap may still develop with PAT, it
still generally finds parameters that lead to good performance, even if those
parameters are sometimes different from those that would have been obtained if
ideal backpropagation could be applied to the physical system.

In sum, PAT corrects backpropagation of physical systems, compared to
performing backpropagation exclusively with a model, by preventing the build-
up of a simulation-reality gap in several different ways. First, it exactly
calculates the output of a sequence of physical operations, as well as each
intermediate output in the sequence, because it uses the actual physical
system to do this. Second, it improves the estimation of parameter updates
based on gradient descent, because the deviation of the true physical output
from the intended output is exactly known, and because gradients are always
estimated with respect to the correct inputs.

PAT seems to do many other things that are helpful, that are mostly beyond the
scope of this work, but that are worth mentioning briefly here. These other
helpful things include: making sure the trained model works despite the
existence of noise in the forward pass, allowing generic models to be used
even when there are device-device variations, and allowing training to proceed
accurately even when the model fails to account for parts of the true physics.
The first occurs because PAT includes a perfect model of the physical system’s
noise, inherently (the physical system itself). During training, if the
physical system’s output randomly fluctuates, estimated parameter updates will
also fluctuate. Training’s progress along ‘noisy’ dimensions (i.e., dimensions
for which the loss strongly fluctuates due to noise) in parameter space will
therefore be slowed depending on how noisy those dimensions are. This is
essentially a generalization of how quantization-aware training [1]
facilitates the training of deep neural networks that can tolerate low-
precision weights and low-precision operations. PAT allows generic models to
be used on heterogeneous devices because it does not require a perfect model
for the physical system (though it does require access to the true physical
system). This may be helpful when training mass-manufactured devices whose
individual parameters may vary. Lastly, when the model is not accurate in
particular regimes, PAT seems to see this as noise, since the estimated
parameter updates will average to approximately zero along ‘mis-modelled’
directions, due to the fluctuating disagreement between the physics and model.
This is probably a generalization of the phenomenon of feedback alignment
encountered when training deep neural networks with random synaptic feedback
[2].

###  1.2 General formulation of physics-aware training

Physics-aware training (PAT) is a gradient-based learning algorithm. The
algorithm computes the gradients of the loss w.r.t. the parameters of the
network. Since the loss indicates how well the network is performing at its
machine learning task, the gradients of the loss are subsequently used to
update the parameters of the network. The gradients are computed efficiently
via the backpropagation algorithm, which we briefly review below.

The backpropagation algorithm is commonly applied to neural networks composed
of differentiable functions. It involves two key steps: a forward pass to
compute the loss and a backward pass to compute the gradients with respect to
the loss. The mathematical technique underpinning this algorithm is reverse-
mode automatic differentiation (autodiff). Here, each differentiable function
in the network is an autodiff function which specifies how signals propagate
forward through the network and how error signals propagate backward. Given
the constituent autodiff functions and a prescription for how these different
functions are connected to each other, i.e., the network architecture,
reverse-mode autodiff is able to compute the desired gradients iteratively
from the outputs towards the inputs and parameters (heuristically termed
“backward”) in an efficient manner.

For example, the output of a conventional deep neural network is given by
$f(f(\ldots
f(f(\bm{x},\bm{\theta}^{[1]}),\bm{\theta}^{[2]})\ldots,\bm{\theta}^{[N-1]}),\bm{\theta}^{[N]})$.
$f$ denotes the constituent autodiff function and is given by
$f(\bm{x},\bm{\theta})=\text{Relu}(W\bm{x}+\bm{b})$, where the weight matrix
$\bm{W}$ and the bias $\bm{b}$ are the parameters of a given layer, and Relu
is the rectifying linear unit activation function (other activation functions
may be used). Given a prescription for how the forward and backward pass is
performed for $f$, the autograd algorithm is able to compute the overall loss
of the network. In Sec. 1.3, the full training procedure for simple
feedforward NN of this kind is described in more detail.

The conventional backpropagation algorithm uses the same function for the
forward and backward passes (see Fig. S1). Mathematically, the forward pass,
which maps the input and parameters into the output, is given by

$\displaystyle\bm{y}=f(\bm{x},\bm{\theta}),$ (1)

where $\bm{x}\in\mathbb{R}^{n}$ is the input, $\bm{\theta}\in\mathbb{R}^{p}$
are the parameters, $\bm{y}\in\mathbb{R}^{m}$ is the output of the map, and
$f:\mathbb{R}^{n}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{m}$ represents
some general function that is an constituent operation which is applied in the
overall neural network.

The backward pass, which maps the gradients w.r.t. the output into gradients
w.r.t. the input and parameters, is given by the following Jacobian-vector
product:

$\displaystyle\frac{\partial L}{\partial\bm{x}}$
$\displaystyle=\left(\frac{\partial
f}{\partial\bm{x}}(\bm{x},\bm{\theta})\right)^{\mathrm{T}}\frac{\partial
L}{\partial\bm{y}},$ (2) $\displaystyle\frac{\partial L}{\partial\bm{\theta}}$
$\displaystyle=\left(\frac{\partial
f}{\partial\bm{\theta}}(\bm{x},\bm{\theta})\right)^{\mathrm{T}}\frac{\partial
L}{\partial\bm{y}},$ (3)

where $\frac{\partial L}{\partial\bm{y}}\in\mathbb{R}^{m},\frac{\partial
L}{\partial\bm{x}}\in\mathbb{R}^{n},\frac{\partial
L}{\partial\bm{\theta}}\in\mathbb{R}^{p}$ are the gradients of the loss w.r.t.
the output, input and parameters respectively. $\frac{\partial
f}{\partial\bm{x}}(\bm{x},\bm{\theta})\in\mathbb{R}^{m\times n}$ denotes the
Jacobian matrix of the function $f$ w.r.t. $\bm{x}$ evaluated at
$(\bm{x},\bm{\theta})$, i.e. $\left(\frac{\partial
f}{\partial\bm{x}}(\bm{x},\bm{\theta})\right)_{ij}=\frac{\partial
f_{i}}{\partial x_{j}}(\bm{x},\bm{\theta})$. Analogously,
$\left(\frac{\partial
f}{\partial\bm{\theta}}(\bm{x},\bm{\theta})\right)_{ij}=\frac{\partial
f_{i}}{\partial\theta_{j}}(\bm{x},\bm{\theta})$.

Though the backpropagation algorithm given by (1-3) works well for deep
learning with digital computers - it cannot be naively applied to physical
neural networks. This is because the analytic gradients of input-output
transformations performed by physical systems (such as $\frac{\partial
f}{\partial\bm{x}}(\bm{x},\bm{\theta})$) cannot be computed. The derivatives
could be estimated via a finite-difference approach - requiring $n+p$ number
of calls to the physical system per backward pass. The resulting algorithm
would be $n+p$ times slower - this is significant for reasonably sized
networks with thousands or millions of trainable parameters.

In physics-aware training, we propose an alternative implementation of the
backpropagation algorithm that is suitable for arbitrary real physical input-
output transformations. This variant employs autodiff functions that utilizes
different functions for the forward and backward pass. As shown schematically
in Fig. S1B, we use the non-differentiable transformation of the physical
system $(f_{\text{p}})$ in the forward pass and use the differentiable digital
model of the physical transformation $(f_{\text{m}})$ in the backward pass.
Since physics-aware training is a backpropagation algorithm, which efficiently
computes gradients in an iterative manner, it is able to train the physical
system efficiently. Moreover, as it is formulated in the same paradigm of
reverse-mode automatic differentiation, a user can define these custom
autodiff functions in any conventional machine learning library (such as
PyTorch) to design and train physical neural networks with the same workflow
adopted for conventional neural networks.

Formally, physics-aware training is described as follows. For a constituent
physical transformation in the overall physical neural network, the forward
pass operation of this constituent is given by

$\displaystyle\bm{y}=f_{\text{p}}(\bm{x},\bm{\theta}).$ (4)

As a different function is used in the backward pass than in the forward pass,
the autodiff function is no longer able to backpropagate the gradients at the
output layer to the exact gradients at the input layer. Instead it strives to
approximate the backpropagation of the exact gradients. Thus, the backward
pass is given by

$\displaystyle g_{\bm{x}}$ $\displaystyle=\left(\frac{\partial
f_{\text{m}}}{\partial\bm{x}}(\bm{x},\bm{\theta})\right)^{\mathrm{T}}g_{\bm{y}},$
(5) $\displaystyle g_{\bm{\theta}}$ $\displaystyle=\left(\frac{\partial
f_{\text{m}}}{\partial\bm{\theta}}(\bm{x},\bm{\theta})\right)^{\mathrm{T}}g_{\bm{y}},$
(6)

where $g_{\bm{y}}$, $g_{\bm{x}}$, and $g_{\bm{\theta}}$ are estimators of the
gradients $\frac{\partial L}{\partial\bm{y}}$, $\frac{\partial
L}{\partial\bm{x}}$, and $\frac{\partial L}{\partial\bm{\theta}}$
respectively.

In listing 1, we show how to construct these custom autograd function for
physics-aware training in PyTorch, a popular library for machine learning.
With these functions at hand, the user can proceed to build physical neural
network with complex architectures in PyTorch. (An implementation of physics-
aware training in PyTorch is available at: https://github.com/mcmahon-
lab/Physics-Aware-Training.)

Figure S1: Schematic for the autodiff function in conventional backpropagation
(A) and in physics-aware training (B). In the forward pass, input and
parameters are propagated to the output. In the backward pass, gradients
w.r.t. the output are backpropagated to gradients w.r.t. the input and
parameters. In physics-aware training, the backward pass is estimated using a
differentiable digital model, $f_{\text{m}}(\bm{x},\bm{\theta})$, while the
forward pass is implemented by the physical system,
$f_{\text{p}}(\bm{x},\bm{\theta}$).

⬇

1 import torch

2

3 def make_pat_func(forward_f, backward_f):

4 ”””

5 A function that constructs and returns the custom autograd function for
physics-aware training.

6

7 Parameters

8 ———-

9 f_forward: The function that is applied in the forward pass.

10 Typically, the computation is performed on a physical system that is
connected and controlled by

11 the computer that performs the training loop.

12 It takes in an input PyTorch tensor x, a parameter PyTorch tensor theta and
returns the output PyTorch tensor y.

13 f_backward: The function that is employed in the backward pass to propagate
estimators of gradients.

14 It takes in an input PyTorch tensor x, an input PyTorch theta and returns
the output PyTorch tensor y.

15

16 Returns

17 ——-

18 f_pat: The custom autograd function for physics-aware training.

19 It takes in an input PyTorch tensor x, an input PyTorch theta and returns
the output PyTorch tensor y.

20

21 ”””

22 class func(torch.autograd.Function):

23 @staticmethod

24 def forward(ctx, x, theta): #here ctx is an object that stores data for the
backward computation.

25 ctx.save_for_backward(x, theta) #ctx is used to save the input and
parameter of this function.

26 return f_forward(x, theta)

27

28 @staticmethod

29 def backward(ctx, grad_output):

30 x, theta = ctx.saved_tensors #load the input and parameters that are stored
in the forward pass

31 torch.set_grad_enabled(True) #autograd has to be enabled to perform
jacobian computation.

32 # Perform vector jacobian product of the backward function with PyTorch.

33 y = torch.autograd.functional.vjp(f_backward, (x, theta), v=grad_output)

34 torch.set_grad_enabled(False) #autograd should be restored to off state
after jacobian computation.

35 return y[1]

36 f_pat = func.apply

37 return f_pat

Listing 1: A helper function that generates the custom autograd function for
physics-aware training. This code was designed to run on Python 3.7 and
PyTorch 1.6 .

Figure S2: Schematic of feedforward physical neural network (PNN). In this
architecture, the output of a physical system is passed as an input to the
subsequent physical system, i.e. $\bm{x}^{[l+1]}=\bm{y}^{[l]}$. A parameter
$\bm{\theta}^{[l]}$ is fed into each layer of the PNN to train the
transformation that the PNN applies to the input $\bm{x}^{[l]}$. Figure S3:
Schematic of the full training loop for the different training algorithms
applied to feedward PNNs. All three algorithms are backpropagation algorithms
and involves four key steps, namely: the forward pass, computing the error
vector, the backward pass and the update of the parameters. The main
difference between the algorithms is whether the physical transformation
$f_{\text{p}}$ or $f_{\text{m}}$ is used for the forward and backward passes.
For instance for _in silico_ training (panel B), $f_{\text{m}}$ is used in the
forward and backward passes. Thus, the intermediate quantities in this
algorithm $\tilde{\bm{x}}^{[l]}$ are different from that of the ideal
backpropagation algorithm (panel A) and PAT (panel C). The blue dashed arrows
indicate that physical system outputs are saved during the forward pass and
the Jacobians will be evaluated at the saved locations.

###  1.3 Motivating physics-aware training with multilayer feedforward
architecture

In the previous subsection, we were able to describe physics-aware training in
full generality by casting it in the language of automatic differentiation.
Thus, we assumed familiarity with automatic differentiation and did not
provide a full prescription of the training loop that is performed during
physics-aware training. In this subsection, we pedagogically walk-through all
the steps involved in physics-aware training for a specific PNN architecture.
Here, we specialize to the feedforward PNN architecture as it is the PNN
counterpart of the canonical feedforward neural network (multilayer
perceptrons) in deep learning. Similar to Sec. 1.2, we begin by reviewing the
full training loop for the conventional backpropagation algorithm. We show why
it is ineffective for PNNs, before introducing the alternative training
algorithms: _in silico_ training and physics-aware training.

From this section onwards, we call the conventional backpropagation algorithm
applied to the physical transformations in PNNs the “ideal backpropagation
algorithm”. This renaming is motivated by two reasons. First, as was pointed
out in the last section, physical systems cannot be autodifferentiated, so it
is impractical to apply this algorithm to real PNNs. Second, to avoid
nomenclature confusion with in silico training, which applies the conventional
backpropagation algorithm to a model of the system.

In supervised machine learning, a training dataset defined by the set of pairs
of examples and targets
$\set{(\bm{x}^{(k)},\bm{y}^{(k)})}_{k=1}^{N_{\text{data}}}$ is provided.
During training, the goal is to obtain parameters of the network that maps any
given example $\bm{x}^{(k)}$ to an output (prediction) $\hat{\bm{y}}^{(k)}$
that is as close as possible to the desired target $\bm{y}^{(k)}$. In
gradient-based learning algorithms, this training is performed by minimizing
the following overall loss function $L_{\text{total}}$ via gradient descent:

$\displaystyle
L_{\text{overall}}=\frac{1}{N_{\text{data}}}\sum_{k}L^{(k)}=\frac{1}{N_{\text{data}}}\sum_{k}\ell(\hat{\bm{y}}^{(k)},\bm{y}^{(k)}),$
(7)

where $\ell$ is the loss function of choice. For instance, in the case of
mean-square error loss, $\ell(\hat{\bm{y}},\bm{y})=|\hat{\bm{y}}-\bm{y}|^{2}$.
As the overall loss function is an average over the different individual
examples, the gradient of the loss with respect to parameters of the networks
are also averaged over individual examples. Thus, following the approach taken
in Ref. [3], the following text will focus on a single example for notational
simplicity; we will on occasions present expressions encompassing all examples
for completeness.

As shown in Fig. S2, the feedforward PNN is composed of $N$ physical systems.
A given example $\bm{x}$ is fed into the input of the physical system at the
first layer $\bm{x}^{[1]}=\bm{x}$. Subsequent physical systems’ input are
taken from the output of the previous layer. Thus, the forward pass is
mathematically given by the following iterative expression:

$\displaystyle\text{Perform Forward
Pass}:\qquad\bm{x}^{[l+1]}=\bm{y}^{[l]}={f}_{\text{p}}(\bm{x}^{[l]},\bm{\theta}^{[l]}).$
(8)

In this architecture, the output of the physical system at the final layer is
the prediction produced by the PNN $\hat{\bm{y}}=\bm{y}^{[N]}$.

To train the parameters of the PNN via gradient descent, we require the
gradient of the loss $L=\ell(\hat{\bm{y}},\bm{y})$ with respect to parameters
of each physical system $\frac{\partial L}{\partial\bm{\theta}^{[l]}}$. These
gradients are computed via the backpropagation algorithm, which iteratively
computes the gradients from the “back”, i.e. from the last layer towards the
first layer. This iterative procedure is more efficient than if
$\frac{\partial L}{\partial\bm{\theta}^{[l]}}$ is computed from scratch for
each layer, which wastes computation on redundant calculations of many
intermediate terms that result from the chain rule. Thus, to begin this
process, we first require the error vector $(\frac{\partial
L}{\partial\hat{\bm{y}}}=\frac{\partial L}{\partial\bm{y}^{[N]}})$, which is
computed via

$\displaystyle\text{Compute Error Vector}:\qquad\frac{\partial
L}{\partial\bm{y}^{[N]}}=\frac{\partial\ell}{\partial{{\bm{y}^{[N]}}}}({{\bm{y}^{[N]}}},\bm{y}).$
(9)

The error vector points out the direction the output of the PNN should change
in to minimize the loss. (Since the goal is to minimize the loss, in practice
one updates parameters along the opposite direction of the gradient.) In the
case of the MSE loss, the error vector takes a particularly intuitive form
given by the difference between the output and the target ($\frac{\partial
L}{\partial\bm{y}^{[N]}}=2(\hat{\bm{y}}-\bm{y})$). This error vector serves as
the “initial” condition for the following backward pass computation:

$\displaystyle[left={\text{Perform Backward Pass}:\qquad}]\frac{\partial
L}{\partial\bm{y}^{[l-1]}}$ $\displaystyle=\left(\frac{\partial
f_{\text{p}}}{\partial\bm{x}}(\bm{x}^{[l]},\bm{\theta}^{[l]})\right)^{\mathrm{T}}\frac{\partial
L}{\partial\bm{y}^{[l]}},$ (10a) $\displaystyle\frac{\partial
L}{\partial\bm{\theta}^{[l]}}$ $\displaystyle=\left(\frac{\partial
f_{\text{p}}}{\partial\bm{\theta}}(\bm{x}^{[l]},\bm{\theta}^{[l]})\right)^{\mathrm{T}}\frac{\partial
L}{\partial\bm{y}^{[l]}}.$ (10b)

To gain an intuitive understanding of Eq. 10, it is important to internalize
that just like the error vector, the gradients of the loss w.r.t. to the
outputs $\frac{\partial L}{\partial\bm{y}^{[l]}}$ (parameters $\frac{\partial
L}{\partial\bm{\theta}^{[l]}}$) points out how the output (parameters) of each
layer should be changed to decrease the loss. Thus, Eq. 10a iteratively
prescribes how an output in the previous layer $l-1$ should change (to reduce
the loss) given information about how an output in the current layer $l$
should change (to reduce the loss). Eq. 10b then takes this computed
information about how outputs should be changed, to compute how the parameters
in each layer should be changed to induce that desired change in the output of
that layer.

With the gradient of the loss w.r.t. parameters, the parameters can be updated
with any gradient-descent optimizer. Though advanced optimizers such as Adam
[4] and AdaDelta [5] typically results in faster training (and are used in our
work), this subsection will focus on the canonical gradient descent update
rule (without momentum) for simplicity. In this framework, the parameters are
updated as follows:

$\displaystyle\text{Update
parameters}:\qquad\bm{\theta}^{[l]}\rightarrow\bm{\theta}^{[l]}-\eta\frac{1}{N_{\text{data}}}\sum_{k}\frac{\partial
L^{(k)}}{\partial\bm{\theta}^{[l]}},$ (11)

where $\eta$ is the learning rate and $\frac{\partial
L^{(k)}}{\partial\bm{\theta}^{[l]}}$ denotes the gradient vector computed for
the $k$th example $\bm{x}^{(k)}$. Consistent with (7), the parameter update
involves an average over the different examples.

Although the ideal backpropagation algorithm (8 \- 11) would in principle
arrive at good parameters of the physical systems that correspond to low loss,
it is not practical in practice as it cannot be performed efficiently. As
mentioned in Sec. 1.2, this is because the algorithm requires analytic
gradients of the physical transformation $f_{\text{p}}$, which can only be
approximated via a finite-difference approach. More specifically, each
backward call requires $n+p$ number of repeated calls to the physical system,
making the training prohibitively slow for PNNs with large number of
parameters.

_In silico_ training circumvents this issue by training the parameters of a
network with a differentiable digital model $f_{\text{m}}$ of the physical
transformation. As the digital model is differentiable, the analytic Jacobian
required for backpropagation can be computed in roughly the same time required
for the forward pass (without repetition). Hence, _in silico_ training could
also be referred to as training performed in numerical simulations, or
training performed exclusively with a numerical model. In _in silico_
training, the training happens solely on a digital computer, to arrive at a
set of parameters for the PNNs. After training, these parameters are loaded on
to physical system for evaluation. Thus, the algorithm is equivalent to
performing conventional backpropagation with $f_{\text{m}}$ as the autodiff
function. Thus, the full training loop is given by

Perform Forward Pass
$\displaystyle:\qquad\tilde{\bm{x}}^{[l+1]}=\tilde{\bm{y}}^{[l]}=f_{\text{m}}(\tilde{\bm{x}}^{[l]},{\bm{\theta}}^{[l]}),$
(12) Compute Error Vector
$\displaystyle:\qquad\frac{\partial\tilde{L}}{\partial\tilde{\bm{y}}^{[N]}}=\frac{\partial\ell}{\partial{\tilde{\bm{y}}}^{[N]}}({\tilde{\bm{y}}}^{[N]},\bm{y}),$
(13) Perform Backward Pass $\displaystyle:\qquad\begin{aligned}
\frac{\partial\tilde{L}}{\partial\tilde{\bm{y}}^{[l-1]}}&=\left(\frac{\partial
f_{\text{m}}}{\partial\bm{x}}(\tilde{\bm{x}}^{[l]},{\bm{\theta}}^{[l]})\right)^{\mathrm{T}}\frac{\partial\tilde{L}}{\partial\tilde{\bm{y}}^{[l]}},\\\
\frac{\partial\tilde{L}}{\partial{\bm{\theta}}^{[l]}}&=\left(\frac{\partial
f_{\text{m}}}{\partial\bm{\theta}}(\tilde{\bm{x}}^{[l]},{\bm{\theta}}^{[l]})\right)^{\mathrm{T}}\frac{\partial\tilde{L}}{\partial\tilde{\bm{y}}^{[l]}},\end{aligned}$
(14) Update parameters
$\displaystyle:\qquad{\bm{\theta}}^{[l]}\rightarrow{\bm{\theta}}^{[l]}-\eta\frac{1}{N_{\text{data}}}\sum_{k}\frac{\partial\tilde{L}^{(k)}}{\partial{\bm{\theta}}^{[l]}},$
(15)

where $\tilde{\bm{x}}^{[l]},\tilde{\bm{y}}^{[l]}$ denotes the input and output
of each layer that is predicted by the differentiable digital model. With the
predicted output of the PNN ($\tilde{\bm{y}}^{[N]}$), the predicted loss
$\tilde{L}$ is computed and minimized in _in silico_ training via parameter
updates governed by gradients of predicted loss $\tilde{L}$
($\frac{\partial\tilde{L}}{\partial{\bm{\theta}}^{[l]}}$) instead of the true
loss $L$ ($\frac{\partial L}{\partial\bm{{\theta}}^{[l]}}$). Therefore, the
angle between the gradients vectors $\measuredangle(\frac{\partial
L}{\partial\bm{{\theta}}^{[l]}},\frac{\partial\tilde{L}}{\partial{\bm{\theta}}^{[l]}})$
conceptually characterizes the performance of the _in silico_ training
algorithm. When this angle is above $90^{\circ}$, training the PNN (via a
parameter update) would decrease the predicted loss, but actually increase the
true error, and _degrade_ the performance of the PNN.

We use the angle between gradients and the gradient obtained by ideal
backpropagation as a conceptual tool here, rather than a readily computable
metric of performance since the ideal backpropagation algorithm is impractical
to implement in experiment. Ultimately, the best metric of training algorithms
is their achieved performance. In addition, since there are usually many
equivalently-performing local optima in training large neural networks, it is
not necessarily an indication that a backpropagation algorithm performs poorly
if it does not follow the exact same training trajectory as ideal
backpropagation. This angle is only evaluated in Sec. 1.4 as we treat a
hypothetical numerical example, which allows for the “true” gradient to be
computed via autograd. In addition, the gradient vectors of the two differing
algorithm must be evaluated for the same parameters and the same PNN input
$\bm{x}$ (example of the ML dataset). Finally, in addition to the angle, the
difference between the magnitude of the gradient vectors $(|\frac{\partial
L}{\partial\bm{{\theta}}^{[l]}}|-|\frac{\partial\tilde{L}}{\partial{\bm{\theta}}^{[l]}}|)$
is also an important metric as the magnitude of the gradient vector determines
how big a step is taken during training. Thus, if the magnitude of the
predicted gradient is vastly off, it can lead to training instabilities.

With the angle metric in mind, it becomes apparent that _in silico_ training
has a decisive flaw; the error vector
$\frac{\partial\tilde{L}}{\partial\tilde{\bm{y}}^{[N]}}$, which is the initial
condition of the backward pass, is computed via the predicted output
$\tilde{\bm{y}}^{[N]}$ instead of the true output $\bm{y}^{[N]}$. Thus,
without an exceptionally precise digital model, this predicted error vector
would point in a significantly different direction from the true error vector.
This inaccurate error vector would then be backpropagated via (14), to result
in bad parameter updates. This problem is exacerbated for deep PNNs, as the
difference between the predicted output of the digital model and the physical
system scales exponentially w.r.t. to depth. In addition to the inaccurate
error vector, the Jacobian matrices that are used to backpropagate the
gradients in the intermediate layers are also inaccurate as they evaluated at
the inaccurate points, at $\tilde{\bm{x}}^{[l]}$ instead of $\bm{x}^{[l]}$.

Physics-aware training is specifically designed to avoid the pitfalls of both
the ideal backpropagation algorithm and _in silico_ training. It is a hybrid
algorithm involving computation in both the physical and digital domain. More
specifically, the physical system is used to perform the forward pass, which
alleviates the burden of having the differential digital models be
exceptionally accurate (as in _in silico_ training). The differentiable
digital model is only utilized in the backward pass to complement parts of the
training loop that the physical system cannot perform.

As described in Section. 1.2, physics-aware training can be formalized by the
use of custom constituent autodiff functions in an overall network
architecture. In the case of the feedforward PNN, the autodiff algorithm with
these custom functions results in and simplifies to the following training
loop:

Perform Forward Pass
$\displaystyle:\qquad\bm{x}^{[l+1]}=\bm{y}^{[l]}={f}_{\text{p}}(\bm{x}^{[l]},\bm{\theta}^{[l]}),$
(16) Compute Error Vector $\displaystyle:\qquad
g_{\bm{y}^{[N]}}=\frac{\partial
L}{\partial\bm{y}^{[N]}}=\frac{\partial\ell}{\partial{{\bm{y}^{[N]}}}}({{\bm{y}^{[N]}}},\bm{y}),$
(17) Perform Backward Pass $\displaystyle:\qquad\begin{aligned}
g_{\bm{y}^{[l-1]}}&=\left(\frac{\partial
f_{\text{m}}}{\partial\bm{x}}(\bm{x}^{[l]},\bm{\theta}^{[l]})\right)^{\mathrm{T}}g_{\bm{y}^{[l]}},\\\
g_{\bm{\theta}^{[l]}}&=\left(\frac{\partial
f_{\text{m}}}{\partial\bm{\theta}}(\bm{x}^{[l]},\bm{\theta}^{[l]})\right)^{\mathrm{T}}g_{\bm{y}^{[l]}},\end{aligned}$
(18) Update parameters
$\displaystyle:\qquad\bm{\theta}^{[l]}\rightarrow\bm{\theta}^{[l]}-\eta\frac{1}{N_{\text{data}}}\sum_{k}g_{\bm{\theta}^{[l]}}^{(k)},$
(19)

where $g_{\bm{\theta}^{[l]}},g_{\bm{y}^{[l]}}$ are estimators of the “exact”
gradients $\frac{\partial L}{\partial\bm{\theta}^{[l]}},\frac{\partial
L}{\partial\bm{y}^{[l]}}$ respectively. Thus, analogous to the metric defined
for _in silico_ training, the angle between the estimated gradient vector and
the exact gradient vector $\measuredangle(\frac{\partial
L}{\partial\bm{{\theta}}^{[l]}},g_{\bm{\theta}^{[l]}})$ conceptually
characterizes the performance of physics-aware training. In physics-aware
training, the error vector is exact $(g_{\bm{y}^{[N]}}=\frac{\partial
L}{\partial\bm{y}^{[N]}})$ as the forward pass is performed by the physical
system. This error vector is then backpropagated via (18), which involves
Jacobian matrices of the differential digital model evaluated at the “correct”
inputs ($\bm{x}^{[l]}$ instead of $\tilde{\bm{x}}^{[l]}$) at each layer. Thus,
in addition to utilizing the output of the PNN ($\bm{y}^{[N]}$) via physical
computations in the forward pass, intermediate outputs ($\bm{y}^{[l]}$) are
also utilized to facilitate the computation of accurate gradients in physics-
aware training.

###  1.4 Numerical example

In this subsection, we present a numerical example to demonstrate how PAT
efficiently trains PNNs. The numerical example presented here is a simplified,
simulated emulation of the experimental results presented in Figure 2 of the
main manuscript. Thus, we train a PNN to perform the vowel classification task
with a feedforward mult-layer architecture, where the constituent hypothetical
physical transformation is an autocorrelation of the input and parameters (see
(21)), which is a crude approximation for the physical transformation realized
with broadband optical second harmonic generation. We construct a
differentiable digital model for this hypothetical physical transformation,
which is then employed in PAT. The performance of PAT is compared against the
performance achieved by training with ideal backpropagation and _in silico_
training. We show the training curves and how the computation is performed in
the forward and backward pass, for the different algorithms.

In the vowel task, the input vowel formant vectors are 12-dimensional, while
the output required for predictions is 7-dimensional, since we consider 7
vowel classes and use the one-hot encoding. To accommodate this, we consider
the feedforward PNN introduced previously, adding some minor pre-processing
and post-processesing of the initial input and final output respectively.

More concretely, we consider a feedforward PNN with a constituent physical
function $f_{\text{p}}$ that has 24-dimensional inputs and outputs. Hence, the
input vowel vector is repeated twice to produce a 24-dimensional input for the
first layer of the PNN and selected portions of the final output are summed to
produce a 7-dimensional vector. Mathematically, the PNN is given by

$\displaystyle{\bm{x}}^{[1]}$
$\displaystyle=[x_{1},x_{1},x_{2},x_{2},\ldots,x_{12},x_{12}]^{\mathrm{T}},$
(20a) $\displaystyle\bm{x}^{[l+1]}$
$\displaystyle=\bm{y}^{[l]}={f}_{\text{p}}(\bm{x}^{[l]},\bm{\theta}^{[l]})\qquad\text{for}\
l=1\ldots N,$ (20b) $\displaystyle\bm{o}$
$\displaystyle=[y^{[N]}_{5}+y^{[N]}_{6},y^{[N]}_{7}+y^{[N]}_{8},\ldots,y^{[N]}_{17}+y^{[N]}_{18}]^{\mathrm{T}}.$
(20c)

where $\bm{o}$ denotes the final output of the PNN 111Here we denote the
output by $\bm{o}$ instead of $\hat{\bm{y}}$ as the output of the PNN has not
been made to approximate the target $\bm{y}$ since it has not been converted
to a probability distribution via the softmax function. and $y^{[l]}_{i}$
denotes the value of the neuron $i$ ($i$th index of vector $\bm{y}^{[l]}$) at
layer $l$. As is typical in conventional NNs, the predicted label for this PNN
is given by the index of the final output ${\bm{o}}$ that results in the
largest value. For the constituent physical function $f_{\text{p}}$, we choose
a mathematical form that roughly mirrors the transformation for the SHG
experiment conducted in Figure 2 of the main manuscript. Hence,
$f_{\text{p}}(\bm{x},\bm{\theta})=a\bm{s}/\max(\bm{s})+c$, where $a$ and $c$
are trainable factors and offsets. Here $\bm{s}$ is an autocorrelation of the
control signal,

$\displaystyle s_{j}=\sum_{i=0}^{\min(j,N-j-1)}q_{j-i}q_{j+i}$ (21)

where
$\bm{q}=\text{clamp}([x_{1},\ldots,x_{N},\theta_{1},\ldots,\theta_{N}]^{T})$
represents the overall control applied to a given layer, which is clamped
between 0 and 1 via the clamp operation ($\text{clamp}(v)=\max(0,\min(v,1))$).
With $f_{\text{p}}$ at hand, we train a differentiable model $f_{\text{m}}$ of
the physical transformation via a data driven approach. A deep neural network
(DNN) with 3 hidden layers that have 1000, 500, 300 hidden neurons
respectively was trained on 2000 instances of physical input vectors (i.e.,
the vector containing inputs and parameters) and output pairs generated by
repeated application of $f_{\text{p}}$.

The PNN is trained by minimizing the following loss function:

$\displaystyle L=\underbrace{H(\bm{y},\hat{\bm{y}})}_{L_{\text{cross-
entropy}}}+\underbrace{\lambda\sum_{l}\sum_{i}\max(0,q_{i}^{[l]}-v_{\text{max}})-\min(0,q_{i}^{[l]}-v_{\text{min}})}_{L_{\text{constraint}}}$
(22)

where $H$ is the cross-entropy function
($H(\bm{p},\bm{q})=-\sum_{i}p_{i}\log(q_{i})$) and
$\hat{\bm{y}}=\text{Softmax}(\bm{o})$. The first term ($L_{\text{cross-
entropy}}$) is the cross-entropy loss responsible for training the output of
the PNN to approach the desired target, while the second term
($L_{\text{constraint}}$) is a regularization term, with regularization
coefficient $\lambda$, to constrain the controls $\bm{q}^{[l]}$ of each layer
to be between $v_{\text{min}}$ and $v_{\text{max}}$. For this numerical
example $\lambda=0.02$, $v_{\text{min}}=0$, and $v_{\text{max}}=1$.

With this loss function, the PNN was trained with the ideal backpropagation
algorithm, _in silico_ training, and physics-aware training. In Fig. S4, we
show the training curves for the 3 training algorithms. Initially, all 3
algorithms are able to train the parameters to reduce the cross-entropy loss
(see inset of Fig. S4B). However, by about $\sim$350 epochs, _in silico_
training is no longer able to train the physical system effectively due to the
increasing differences between the outputs of $f_{\text{p}}$ and
$f_{\text{m}}$. PAT initially follows the training curve of the ideal backprop
algorithm quite closely. By about $\sim 500$ epochs, the effective rate of
learning for PAT slows down but nevertheless it is able to lower the loss
consistently until the test classification accuracy is $\sim$96%, which is
close to that of ideal backprop. After $\sim$2000 epochs, the cross-entropy
loss and classification accuracies exhibits fluctuations (this is also seen in
Figure 2 of the main manuscript). We speculate that minor differences in the
Jacobian of $f_{\text{p}}$ and $f_{\text{m}}$ causes the parameters to drift
around when they are located near an optimal parameter region of the overall
loss function.

We also characterize the inner workings of the PNNs for the different
algorithms. To do so, we consider two perspectives. The first looks at how the
algorithms differ over just one training step. This perspective, which is
presented in Fig. S5, assumes all algorithms initially have the exact same
parameters, and we then take 1 training step (i.e., 1 parameter update based
on 1 mini-batch of training examples). The second perspective analyzes
divergence of training trajectories over many training steps. This
perspective, which is shown in Figs. S6 and S7, shows how the single-step
inaccuracy of in silico training compounds through training leading to
increasingly worse training of the physical system, and shows how PAT avoids
this divergence.

In Fig. S5, we show how the forward and backward passes of the PNN differ
across the different algorithms when the parameters are set to initially be
the same. Here we observe that the prediction error between the outputs of
$f_{\text{m}}$ and $f_{\text{p}}$ build up progressively through layers during
the forward pass. This results in an incorrect error vector for _in silico_
training, that is backpropagated to result in incorrect gradients of the loss
w.r.t. parameters. As shown in the last row of Fig. S5, the gradients of _in
silico_ training do not resemble the “true” gradient of the ideal
backpropagation algorithm (which can only be obtained from a finite-difference
approach in a laboratory setting), in both direction and magnitude. In
contrast, PAT is able to approximate the gradients fairly accurately, as it
utilizes the experimental transformation $f_{\text{p}}$ in the forward pass.
This may be surprising because PAT uses the exact same model as in silico
training to estimate gradients. The key is that the gradients with PAT are
estimated at the correct locations in parameter space (i.e., the inputs to
each physical layer are known exactly), whereas in in silico training, the
gradients are calculated at locations that are increasingly mismatched with
reality as one goes deeper into the network.

In Fig. S6 and S7, we characterize the PNNs after training for 300 and 800
epochs respectively, for the different training algorithms. We show how the
vowel data propagates through the different layers of the PNNs to arrive at
the final prediction $\hat{\bm{y}}$ for 3 different vowel examples. In Fig.
S6, we see that the ideal backpropagation algorithm learns the fastest as it
is already able to correct classify 2 of out the 3 vowels. The parameters
learned by PAT are also reasonably similar to those obtained with ideal
backpropagation, and the trained model performs similarly well. In silico
training meanwhile trains slowly, but by 300 epochs has made some small
positive progress. When the PNN is executed with the parameters obtained by in
silico training after 300 epochs, the PNN’s output is shifted slightly towards
the target vowel. Thus, consistent with Fig. S4, learning has occured for all
3 algorithms after 300 epochs of training. As training proceeds, however, the
simulation-reality mismatch from each training step compounds, causing in
silico training to diverge rapidly away from ideal backpropagation and PAT.
This can be exacerbated because in silico training can often proceed into
regions of parameter space where its predictions are especially poor, leading
to an even faster divergence from ideal backpropagation and PAT. By epoch 800
(see Fig. S7), the parameters predicted by in silico training have drifted
into a region where the digital model’s predictions are especially inaccurate.
This is shown in the rightmost column of Fig. S7, where the output of the
experiment with parameter trained via _in silico_ training classifies the
vowel wrongly 2 out of 3 times, while the PAT and ideal backpropagation is
able to classify all 3 vowels accurately with a reasonable high confidence,
which is indicated by the value of $\hat{\bm{y}}$ being close to 1 for the
correct vowel.

Figure S4: Comparison of the test classification accuracy (A) and test cross-
entropy loss (B) versus training epoch for ideal backpropagation, _in silico_
training, and physics-aware training. Since this numerical example attempts to
mirror the experiment performed in Figure 2 of the main manuscript, it
considers the vowel classification task with a feedforward PNN, specified by
(20). The inset in B zooms in on the initial training phase, where the cross-
entropy loss for the _in silico_ training briefly decreases for about 300
epoch before increasing due to compounding mismatch between the $f_{\text{p}}$
and $f_{\text{m}}$.  Figure S5: Comparing 3 training algorithms in a numerical
example, for a single training step. More specifically, the inputs, outputs,
parameters and gradients of the feedforward PNN applied to the vowel
classification task are shown in each row of plots respectively, with each
physical layer in the first, second and third columns of plots respectively.
Here, we transfer the same parameter to all 3 algorithm, to evaluate how
similar the gradient vectors are for the different algorithms. Figure S6:
Comparing the training of 3 training algorithms over many training steps.
Here, the training has proceeded for 300 steps. More specifically, we show the
physical inputs at the different layers, the output and the final PNN
prediction of PNNs for three example vowels. In contrast to Fig. S5, each
algorithms here is associated with different parameters that result from each
of their own distinct training trajectories shown in Fig. S4. Plots associated
with “ _In silico training_ evaluated in experiments” refer to
characterizations of the PNN when the parameters trained via _in silico_
training is transferred over and evaluated with the experiment (with the
function $f_{\text{p}}$ instead of $f_{\text{m}}$). Figure S7: Comparing the
training of 3 training algorithms over many training steps. Here, the training
has proceeded for 800 steps. More specifically, we show the physical inputs at
the different layers, the output and the final PNN prediction of PNNs for
three example vowels. In contrast to Fig. S5, each algorithms here is
associated with different parameters that result from each of their own
distinct training trajectories shown in Fig. S4. Plots associated with “ _In
silico training_ evaluated in experiments” refer to characterizations of the
PNN when the parameters trained via _in silico_ training is transferred over
and evaluated with the experiment (with the function $f_{\text{p}}$ instead of
$f_{\text{m}}$).

##  2 Supplementary methods

In this section, we first describe the experimental setups for each of the
three physical systems used to implement PNNs. We then describe our data-
driven differentiable models, as well as the noise models used to improve the
performance of in silico training. Finally, the PNN architectures used in the
main text for classifying vowels and MNIST digits are described and typical
results are presented to illustrate how the PNNs perform the learned
classification tasks.

###  2.1 Femtosecond second-harmonic generation PNN

####  2.1.1 Experimental setup

Figure S8: Labelled photos of the pulse shaper and second harmonic generation
experimental setups.

Our femtosecond second-harmonic generation (SHG) PNN is based on a mode-locked
Titanium:sapphire oscillator (Spectra Physics Tsunami) which is operated on by
a digital micromirror device (DMD, Vialux V-650L)-based pulse shaper to
produce complex, amplitude-modulated pulses according to input data and
parameters. The spectrally-modulated pulses, centered around 780 nm, are then
focused into an 0.5 mm long BBO crystal where SHG occurs. After filtering the
residual 780 nm light, the resulting blue ($\sim$ 390 nm) light is diffracted
by a ruled reflection grating and focused onto a fast camera (Basler ace
acA800-510u) to measure the SHG spectrum. Annotated photos of the setup are
shown in Supplementary Figure S8.

Training PNNs requires executing the physical evolution millions of times, so
it was prudent to design the setup to allow for high-throughput operation.
Based on the design described above, we were able to operate the SHG-PNN at a
sustained rate of 200 Hz, reliably over months (including remotely during the
COVID-19 pandemic).

The DMD-based pulse shaper achieves continuous spectral amplitude modulation
in the following way. Light diffracted by a transmission grating (Ibsen
Photonics PCG-1765-808-981) is focused by a cylindrical lens (150 mm,
Thorlabs) onto the DMD, which is mounted at 45 degrees from the vertical so
the micromirrors’ rotation axis is vertical with respect to the table. The DMD
is then rotated slightly so that the light reflected by the mirrors in the ‘0’
state retroflects backward through the cylindrical lens and grating, following
an identical path except with a slight downward angle so as to be separated by
a pick-off mirror. By writing vertically-oriented (with respect to the table,
45 degrees to the DMD’s axes) gratings of varying duty cycle to the DMD, we
can achieve multi-level amplitude modulation at each spectral position (the
frequency components are dispersed across the DMD in the direction parallel to
the table by the grating). The design of this pulse shaper was inspired by
Ref. [6], and to control the DMD using Python we adapted code developed for
Ref. [7] available at Ref. [8].

While such a design allows continuous amplitude modulation at significantly
higher speeds than most spatial light modulator-based pulse shapers, it
results typically in pulses with complex spatiotemporal coupling which would
be challenging to account for with traditional simulations. In this regard,
our data-driven digital models were helpful in allowing the physical apparatus
to be designed to optimize the reliability and speed of the experiment, rather
than to ensure for the physics to easily be modelled by traditional
simulations.

Figure S9: Characterization of the femtosecond second harmonic generation
(SHG) experimental setup with simple inputs. A-C. Sharply peaked input signals
to the pulse shaper, which results in continous-wave like input waves for SHG.
D-F. The resulting output spectrum of femtosecond SHG for their respective
input signals.  Figure S10: Response of femtosecond SHG device to a complex-
waveform input whose spectrum has a small flat region with an amplitude that
is varied. A. Waveform for pulse shaper to apply to optical pulse, which are
the physical inputs of the SHG system. Different colors correspond to
different values $(A)$ that the small flat region takes on. B. Output spectrum
from SHG. Each curve of a particular color, is the output of the SHG system
subject to input signals with the same color in subpanel A. C. Slices through
output spectrum as a function of the input section amplitude $A$, which shows
the nonlinearity of the transformation implemented by SHG.

####  2.1.2 Input-output transformation characterization

Here, we characterize the input-output transformation of femtosecond SHG. In
Fig. S9, we characterize the SHG setup with simple quasi-cw inputs to
elucidate the basic physics of the device.

In Fig. S9 A and B, we consider continuous-wave-like inputs that are sharply
peaked in wavelength. The corresponding SHG outputs in Fig. S9 D and E are
also sharply peaked in wavelength, with half the wavelength. In panel C, we
consider an input that is the sum of input A and B; the output shown in panel
F, illustrates the nonlinear nature of SHG, as it is not the sum of output D
and E. There is an additional strong output at a center wavelength resulting
from combining both input beams at
$\mathrm{7}\mathrm{7}\mathrm{2}\mathrm{n}\mathrm{m}$ and
$\mathrm{7}\mathrm{7}\mathrm{7}\mathrm{n}\mathrm{m}$.

Having characterized its basic behavior, we illustrate the complex
controllable transformation of the SHG device in Fig. S10, by showing its
response to a complex input waveform. Figure S10A shows input waveforms in
which a small section has been held constant to a fixed amplitude $A$, for
varying values of $A$ between 0 and 1. Figure S10B shows the resulting output
spectrum from the SHG, which illustrates the complexity of the spectrum; The
output can be significantly different across a wide portion of the spectrum
despite the changed portion in the input being relatively small. Figure S10C
plots the output spectrum at a particular wavelength for varying values of the
amplitude $A$, which also shows that the nonlinearity is strongest for large
amplitude $(A\approx 1)$.

###  2.2 Analog transistor PNN

####  2.2.1 Experimental setup

The electronic PNN experiments were carried out with standard bulk electrical
components, a hobbyist circuit breadboard, and a high-speed DAQ device
(Measurement Computing USB-1208-HS-4AO), which allowed for a single 1 MS/s
analog input and output. A schematic of the setup is shown in Supplementary
Figure S11.

A variety of circuits were explored during the course of this work. The one
used ultimately for the MNIST digit classification task was designed to be as
simple as possible while still exhibiting very strongly nonlinear dynamics.
The inclusion of a RLC oscillator within the circuit was inspired by the
relative success of the oscillating plate (described in the next section),
which we found was easily able to achieve good performance on the mostly-
linear MNIST task.

Figure S11: Schematic of the analog electronic circuit used for the electronic
PNN.

####  2.2.2 Input-output transformation characterization

Here we briefly discuss characterization of the circuit’s response.
Supplementary Figure S13 shows the impulse response of the circuit to an
impulse with a varying amplitude. Supplementary Figure S12A shows a complex
input waveform in which a section has been held constant to a fixed amplitude
$A$, for varying values of $A$ from -10 to 10 V. Supplementary Figure S12B
shows the resulting output time series, illustrating both the nonlinearity of
the response as well as the complex coupling between outputs after the changed
value. Finally, Supplementary Figure S12C illustrates the nonlinearity of the
circuit’s response more clearly by plotting the output voltage of the
indicated output time bins as a function of the input amplitude voltage $A$.
The individual curves are highly nonlinear with respect to the input voltages.

At the highest input amplitudes in Supplementary Figure S12C, an additional
nonlinearity due to the limiting of the DAQ can be observed. While this
nonlinearity is not strictly a part of the electronic circuit, we opted to
allow the PNN to operate in this regime to illustrate how PNNs can inherently
utilize such “undesired” nonlinearites, and because the circuit’s response was
already highly nonlinear, so we did not expect this additional nonlinearity to
significantly change the circuit’s computational capabilities (in contrast, it
was important to avoid this DAQ-based nonlinearity in our oscillating plate
experiments described next, since the oscillating plate’s input-output
response was determined to be very linear).

Figure S12: Response of nonlinear analog transistor circuit to a complex
waveform with one varied constant section. A. Input waveforms (different
colors correspond to different input waveforms, where a section in the middle
has been held constant to different fixed voltages). B. Output waveforms
(different colors correspond to the different input traces in A). C. Slices
through output waveforms as a function of the input section amplitude show the
strong nonlinearity of the circuit’s transformation of the input. For the last
figure, the x-axis is plotted in the calibrated, normalized units used within
the PNN PyTorch code, which approximately correspond to -10 to 10 V range from
-1 to 1. Figure S13: Impulse response of the nonlinear analog transistor
circuit to impulses of varying amplitude. The circuit’s response is that of a
nonlinear oscillator. In B, the different curves are offset artificially to
improve visual clarity.

###  2.3 Oscillating plate PNN

####  2.3.1 Experimental setup

The setup of the oscillating plate PNN consists of an audio amplifier (Kinter
K2020A+), a commercially available full-range speaker, a microphone (Audio-
Technica ATR2100x-USB Cardioid Dynamic Microphone), and a computer controlling
the setup.

Figure S14: A schematic of the information flow in the oscillating plate
experimental setup. An input and control signal from the computer is amplified
and applied to the mechanical oscillator (realized by the voice coil of an
acoustic speaker). A microphone records the sound produced by the oscillating
plate and returns it to the computer.

We use the speaker to drive mechanical oscillations of a $3.2\text{ cm}\times
3.2\text{ cm}\times 1\text{ mm}$ titanium plate that is mounted on the
speaker’s voicecoil. The diaphragm has been removed from the speaker such that
most of the sound produced stems from the oscillations of the titanium plate.
In Supplementary Figure S15, the steps taken to construct the oscillator with
the mounted plate are shown.

Figure S15: Photographs of steps taken to construct an audio-frequency
mechanical oscillator from a commercially available speaker. A. First, we
remove the diaphragm and dust cap from a speaker with a precision knife to
expose the voice coil. B. Next, we glue a screw and washer to the voice coil
with commercially available two component glue. C. After letting the glue dry
for 24 hours, we attach the titanium plate on the screw. Finally, we securely
mount the speaker on a stable surface to suppress vibrations of the whole
device. Different speakers were used over the course of the experiment and not
all speakers shown above are the same model.

####  2.3.2 Input and output encoding

Similar to the electronical PNN, we encode the physical inputs in the time
domain. Supplementary Figure S16 shows the different steps of encoding and
decoding signals in the oscillating plate PNN. First, inputs are encoded in a
time-series of rectangular pulses that are transformed into an analog signal
at 192 kS/s by the laptop’s soundcard DAC. The signal is amplified by an audio
amplifier and drives the voicecoil of a speaker that in turn drives the
oscillating titanium plate. We ensured that both the soundcard DAC and the
pre-amplifier are operated in a regime where their input-output relation is
completely linear. The signal arriving at the oscillating plate was therefore
a linearly amplified and slightly low-pass filtered version of the input and
control signal created by the computer. The oscillations of the plate produce
soundwaves which are recorded by a microphone and converted into digital
signals at 192 kS/s. The signal is further downsampled by partitioning it into
a number of equally long subdivisions and averaging the signal’s amplitudes
over this window. The length of the window is determined by the desired output
dimension of the PNN layer.

Inputs and outputs are synchronized by a repeating trigger signal which
precedes every sample that is played on the speaker. By overlapping the
trigger signals we can synchronize samples to about
$5\text{\,}\mathrm{\SIUnitSymbolMicro s}$ (1/sampling frequency).

Figure S16: An example of a 48-dimensional digital input to the oscillating
plate setup (blue). After driving the oscillating plate, the signal is
recorded by the microphone (orange) and digitized in time to the desired
output dimension, here 24 (orange). Y-axis units are arbitrary normalized
amplitudes.

####  2.3.3 Characterization of input-output transformation

While we initially aimed to create high-amplitude, nonlinear oscillations on
the titanium plate, we realized that we could not reach such amplitudes with
our setup. In Supplementary Figure S17, C, we observe that in the time-domain
the output voltage’s response is overwhelmingly linear with respect to the
input voltage.

Neglecting noise and assuming the oscillations of the plate are completely
linear, each input to the plate can be regarded as exciting a band of modes
between the frequency of the input and the Nyquist frequency of the digital
sampler (192 kHz). Following their excitation, the modes decay in a transient
oscillation depending on the damping of the material, thereby affecting the
following outputs in a way that can be expressed as a convolution:
$y_{k}=\sum_{0<j<k}^{N}c_{k-j}x_{j+1}$. Here $y_{k}$ is the $k$th output,
$x_{j}$ is the $j$th input, and the $c_{j}$ coefficients are constants
determined by the medium of the oscillation and the excited modes.

As a matrix operation, the map from inputs to outputs can be approximated by:

$\begin{pmatrix}y_{1}\\\ y_{2}\\\ y_{3}\\\ \vdots\\\
y_{d}\end{pmatrix}=\begin{pmatrix}c_{1}&0&0&\dots&0\\\
c_{2}&c_{1}&0&\dots&0\\\ c_{3}&c_{2}&c_{1}&\dots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
0&0&0&\dots&c_{1}\end{pmatrix}\begin{pmatrix}x_{1}\\\ x_{2}\\\ x_{3}\\\
\vdots\\\ x_{d}\end{pmatrix}$ (23)

For a typical input signal to the oscillating plate PNN, the resulting output
is shown in Fig. S17, A and B. In Fig. S17A, eleven input time-series that are
identical except for 20 inputs at around $0.3$\mathrm{m}\mathrm{s}$$ are
plotted. Inputs change with a frequency of $192$
$\mathrm{k}\mathrm{H}\mathrm{z}$. The transient oscillation set in motion by
the inputs at $0.3$\mathrm{m}\mathrm{s}$$ persists for about
$2$\mathrm{m}\mathrm{s}$$ (not completely shown in the figure). Hence, at this
input frequency, an output of the oscillating-plate PNN is a convolution of
approximately $N=384$ previous inputs. If we encoded inputs at a slower
frequency, for example $5$ $\mathrm{k}\mathrm{H}\mathrm{z}$, we would only
convolve the previous $N=10$ inputs. Therefore, by controlling the frequency
of inputs, the oscillating plate emulates the operation of a convolution with
variable kernel size and fixed kernel coefficients. This feature was utilized
by designing our PNN for MNIST handwritten digit classification.

Figure S17: Response of oscillating plate to a complex waveform with one
varied constant section. A. Input waveforms. B. Output waveforms. Due to
causality, all outputs before $50\text{\,}\mathrm{m}\mathrm{s}$ are the same
(up to noise). Afterwards, different input voltages produce different
transient oscillations. C. Slices through output waveforms as a function of
the input voltage (in arbitrary units) show the overwhelmingly linear response
of the output to the input voltage.

###  2.4 Data-driven differentiable models

####  2.4.1 Digital model for the mean

As explained in Section 1, to approximate the gradients of physical systems,
one needs a differentiable digital model to perform autodifferentiation on in
the backward pass of PAT. In principle, any differentiable model could be
used, such as a traditional differential equation physics model, some other
analytic model, a data-driven model like a deep neural network, or a physics-
informed neural network. In this work, we found that using neural architecture
search [9] to find optimal deep neural network architectures allowed for the
most accurate digital models across the physical systems we considered, and
thus this was our default approach.

Of course, we expect that including physical insight into the differentiable
digital model will usually result in better performance for any specific
system, or at least allow for faster training of the digital model. Since any
system used as a PNN will be ideally capable of being executed many times very
quickly, we expect that it will usually be beneficial to include significant
data-driven components in the digital model, since obtaining a large volume of
training data will be straightforward. In some cases conceivable in the
future, however, it may be more challenging to rapidly change the physical
system’s parameters, and the input-output dimension of the physical system
used may be extremely high, such that even a large training data set
significantly undersamples the space. In these and other cases, it is obvious
that using physical insight to permit accurate differentiable digital models
with relatively little training data will become essential.

The accuracy of the DNN models we obtained for the SHG and nonlinear analog
circuit systems will be shown and commented on in later sections. Here, we
will give a summary overview of our methods and the models ultimately used.

We first collected training data by recording the outputs of the physical
system for a selection of $N$ inputs/parameter vectors, where $N$ was usually
$\sim 10^{4}$ and the input/parameter vectors were usually drawn from a
uniform distribution. We also investigated more sophisticated input
distributions, but found these were generally unnecessary to produce good
results (at least for the results considered in this work). This data set,
consisting of $N$ input/parameter vectors and their corresponding output
vectors, was split into training and validation sets and used to train neural
networks whose hyperparameters were found through a neural architecture
search.

To perform our neural architecture searches, we used the Optuna package [9]
with the following typical search space:

  1. 1.

Number of layers between 1 and 5

  2. 2.

Each layer width between 20 and 1500 units

  3. 3.

Initial learning rates $10^{-5}$ and $10^{-2}$

The swish activation function was found to perform best, and we used the Adam
optimizer.

For the SHG system, the optimal architecture was found to be a deep neural
network with the following layer dimensions: [100, 500, 1000, 300, 50].

For the electronic circuit, the optimal architecture was found to be a deep
neural network with the following layer dimensions: [70, 217, 427, 167, 197,
70].

For the oscillating plate, we initially tried to find a deep neural network
architecture. However, once we realized the system dynamics were linear, we
instead adopted a fully linear model, which was implemented in our deep neural
network training framework by creating a network with no hidden layers.

####  2.4.2 Noise digital model

The digital models we used also include a model of the physical system’s noise
(including its dependence on the particular input/parameter vector to the
system). We used these noise models for two key reasons. First, we were often
interested in simulating PNN architectures and wanted to have a reasonably
accurate prediction of their performance in the experiment. This was
especially useful for slow physical systems, like the oscillating plate, for
which MNIST experiments were much too slow to allow us to perform systematic
examination of different architectures. More importantly, we wanted to ensure
that we maximized the performance of in silico training, and providing an
accurate model of the system’s noise is a logical, standard technique to help
mitigate the simulation-reality gap. All the results for training with
simulation presented in this work (including in Figure 3) were obtained with a
simulation model of the physical system noise.

We find however that while including a sophisticated model of physical
system’s noise in a simulation does improve the ability of in silico training
to estimate the approximate performance level of results obtained in the lab,
including noise in the simulation does not allow in silico training to perform
as well as PAT for training the actual physical system.

A total digital model incorporating both mean-field and noise predictions
needs to learn the following stochastic map:
$y=f_{\text{mean}}(x)+n_{\text{noise}}(x)$, where
$x\in\mathrm{R}^{N_{\text{in}}}$ is the input to the PNN, and $y$ is a non-
deterministic output from the PNN. This output has a deterministic component
(given by the mean-field digital model $f_{\text{mean}}$), and a “noise”
component $n_{\text{noise}}(x)\in\mathrm{R}^{N_{\text{out}}}$ whose
distribution must agree with the true physical system’s input-dependent noise.

Note that the digital model needs to learn the noise distribution from _any_
given input/parameter vector $x$ (for brevity, we will refer to $x$ as just
the input vector from here onward, but note that it includes both the
parameters and input data to the physical system). For simplicity, we first
outline how to learn the noise distribution for a given fixed input vector,
denoted as $x^{\text{(fixed)}}$ before addressing the more general case of
predicting the noise for any input vector.

To learn the noise distribution for this fixed input, we first collect samples
of outputs for the same fixed input. We execute the physical system
$N_{\text{repeat}}$ times with the same input vector $x^{\text{(fixed)}}$. The
output vectors will be denoted $y^{(fixed,k)}$ where
$k=1,...,N_{\text{repeat}}$. From experimental measurements, we have
ascertained that the output distribution of the physical systems considered in
this work are reasonably well described by a multivariate normal distribution.
Thus, their entire noise distribution is well characterized by the mean vector
$\mu_{i}^{(fixed)}=\sum_{k}y^{(fixed,k)}_{i}$ predicted by the mean-field
digital model, and a covariance matrix
$\Sigma^{(fixed)}_{ij}=\frac{1}{N_{\text{repeat}}-1}\sum_{k}(y^{(fixed,k)}_{i}-\mu_{i}^{(fixed)})(y^{(fixed,k)}_{j}-\mu_{j}^{(fixed)})$.

How can we generate typical additive noise on top of the mean-field output
vector, given the covariance matrix? The answer is simple: we need to find a
matrix $A^{(fixed)}$ that satisfies
$A^{(fixed)}{A^{(fixed)}}^{\mathrm{T}}=\Sigma^{(fixed)}$. Once we have the $A$
matrix, a sample of the noise $n\in\mathrm{R}^{N_{\text{out}}}$, which has the
correct distribution, can be obtained from a matrix multiplication between
$A^{(fixed)}$ and a vector of i.i.d normal noise $z_{i}\sim
N(0,\sigma^{2}=1)$, i.e. $n=A^{(fixed)}z$.

So far, the dimensions of $A$ has been left unspecified. In principle, if the
noise is sufficiently simple (e.g., the covariance matrix is rank 1), the
covariance matrix can be decomposed as an outer product between a single
vector and itself, i.e. $A^{(fixed)}\in\mathrm{R}^{N_{\text{out}}\times 1}$.
We can generalize this concept of compressing the representation of the noise
via the following procedure. First perform the spectral decomposition on
$\Sigma^{(fixed)}$ to arrive at $UDU^{\mathrm{T}}$. Then, by looking at the
eigenvalue spectra, ascertain how many dominant eigenvalues are present in the
covariance matrix and denote that number as $N_{\lambda}$. Then let
$A^{(fixed)}=U\sqrt{D}[:,1:N_{\lambda}]$, to arrive at a compressed
$A^{(fixed)}$ matrix with dimension $N_{\text{out}}\times N_{\lambda}$. For
training our noise digital models, we find it is helpful to limit the number
of eigenvectors of the covariance matrix, so that the dimension of the input-
output map that must be learned is as small as possible (and because higher-
order eigenvectors are not usually important).

To summarize, for a given input vector $x^{(fixed)}$, we only need to learn
$A^{(fixed)}\in\mathrm{R}^{N_{\text{out}}\times N_{\lambda}}$ to output
samples of the noise. Thus, to learn a noise digital model for general inputs,
we simply train neural network to, given many different inputs $x^{(q)}$,
predict the corresponding $A^{(q)}$ for any given $x^{(q)}$. This corresponds
to a neural network that outputs a vector corresponding to a flattened $A$ (a
vector of dimension $N_{\text{out}}N_{\lambda}$). Training data for this
network is obtained by measuring the output of the physical system many times
(we found 40 repeats was usually sufficient) for a nominally-fixed $x^{(q)}$
in order to compute the corresponding $A^{(q)}$ (repeating the above procedure
for each training example). When used in the forward pass of the digital model
for training in simulation, the noise model NN’s predicted vector $A$ is
reshaped to the matrix form and then used to generate a noise sample that is
added to the mean-field digital model’s prediction.

To provide examples of typical noise digital models, Supplementary Figures S18
and S19 compare, for the validation set, the performance of the noise digital
model’s prediction of the covariance eigenvectors and the covariance matrix to
the experimental measurements of these quantities. The agreement is somewhat
more approximate than what we obtain for the mean-field digital models (shown
in the next subsections), but is overall excellent.

Figure S18: Characterization of the noise digital model. Here, the digital
model (in black) predicts the 2 dominant noise eigenvectors $\vec{A}_{1}$ and
$\vec{A}_{2}$ for three different 50-dimensional physical inputs to the SHG
device. The experimental noise eigenvectors (in green) are derived from
diagonalizing the sampled covariance matrix from 30 repeated measurements.
Each row represents a distinct physical input. Figure S19: Characterization of
the noise digital model. Left column: The experimental sampled covariance
matrix from 30 repeated measurements. Middle column: The predicted covariance
matrix computed with the predicted dominant noise eigenvectors $\vec{A}_{1}$
and $\vec{A}_{2}$ that are shown in Fig. S18. Right column: The difference
between the predicted and experimental covariance matrix. Each row represents
a distinct physical input, and the covariance matrix for a given row shares a
color bar that is shown on the right. The relatively low amplitudes visible in
the right column show the strong agreement between the digital model’s
predicted noise distribution and the experimentally-measured distribution.

###  2.5 Descriptions of PNN architectures

In this section, we describe the PNN architectures used for the different
machine learning tasks implemented with each physical system. In each
subsection, we also show the performance of the differentiable digital model
used, and (for the PNNs not thoroughly examined in the main paper) present
results illustrating how the physical transformations perform the
classification task.

Before mentioning specifics for each architecture, we mention several
commonalities among all PNNs we present in this work.

First, the same basic training procedure was used in all cases. We found that
the Adadelta optimizer [5] worked best, and speculate that its exclusion of
momentum and parameter-dependent adaptation makes it well-suited to PNNs,
whose parameters are generally not symmetric. We also usually utilized
learning rate scheduling, but this was of secondary importance.

Second, in all cases we compared trained PNNs to baseline reference models in
which all physical input-output transformations were replaced by identity
operations. This allowed us to confirm that, despite each PNN architecture
including some digital operations, the digital operations alone were not
responsible for the majority of the functionality of the trained PNN model.

Third, for convenience, we sometimes performed some digital pooling operations
to reduce the dimension of the measured physical output vectors. This is
because the physical input-output transformations we considered (for
convenience) had all equal input and output dimensions, and this operation
allowed us to modify the dimensionality of the input-output transformations,
allowing us to tackle machine learning problems with different dimensionality
in the input (feature) vector and output vector (for instance, the MNIST
dataset has a 196 input dimension and 10 dimensional output). Since this
operation is not trained and it can only contribute a minor computational
power to the PNN. Besides, any minor benefits it confers can easily be
performed physically as it is equivalent to changing the sampling rate of
output measurements (e.g., by using a slower digital-analog converter).
Finally, we note that for the oscillating plate and electronic PNNs,
downsampling operations of this kind were ultimately irrelevant, as both
devices were operated at input and output sampling rates limited by the sound
card and data acquisition device respectively.

Finally, in all cases it was necessary to include additional Lagrangian terms
in the loss in order to ensure that physical inputs and parameters sent to
physical systems were within the allowed ranges. This additional terms are
identical to the Lagrangian terms that have been discussed earlier in Section
1.4.

####  2.5.1 Figure 2: Femtosecond SHG-vowels PNN

The femtosecond SHG-vowel PNN was chosen to be the first example PNN of the
main manuscript for two reasons. First, the vowel classification task is
simple enough to facilitate a relatively simple PNN, and second, the SHG
physical transformation is relatively unusual, being a fully-nonlinear
transformation with nearly all-to-all coupling of inputs. Although the second
feature does not make SHG well-suited to the vowel-classification task (which
is linear), we felt it served as a clear, small-scale tutorial demonstration
of a multilayer PNN on a relatively high-dimensional task in which the
physical transformation used is highly nonlinear.

To make the SHG-vowel PNN intuitive and illustrative for the concept of PNNs
in general, we made several design choices. We used a portion of the pulse
shaper’s spectral modulations for trainable parameters, and a second, distinct
portion for encoding input data. In total, each SHG input-output
transformation maps 100 physical inputs, which are split between the input
data and trainable parameters, to a 50-dimensional output SHG spectrum.

In Fig. S20, we show some representative input and outputs of this
transformation, where physical inputs are split equally between the input and
parameters (50 dimensional input and parameters each). Consistent with the
analysis performed in Sec. 2.1.2, the transformation is complex, as evidenced
by the widely-varying outputs. Nevertheless, this transformation can be
effectively modeled via a data-driven model. By performing the neural
architecture prescribed in Sec. 2.4, we found a deep neural network that is
able to learn the transformation effectively. This digital model’s excellent
agreement with the experimental input-output transformation is shown in Fig.
S18.

Having characterized the performance of the differentiable digital model, we
now show how the SHG system performs vowel classification. To begin, we first
describe the vowel classification task that we perform. We use the vowel
dataset from Ref. [10], which can be downloaded at this https://static-
content.springer.com/esm/art%3A10.1038%2Fs41586-018-0632-y/MediaObjects/41586_2018_632_MOESM2_ESM.docx,
with an additional (untrainable) renormalization to account for the fact that
the SHG system takes in physical inputs between 0 and 1. For completeness, we
will paraphrase the description of the dataset that is in the supplementary
information of Ref. [10], before describing the renormalization we apply to
the dataset.

The vowel data considered is a subset of the Hillenbrand vowel database, which
can be found at https://homepages.wmich.edu/~hillenbr/voweldata.html. In this
dataset, each spoken vowel is characterized by their formant frequencies,
which refers to the dominant frequencies (in units of $\mathrm{H}\mathrm{z}$)
present in an audio signal. More specifically, each vowel is characterized by
a 12-dimensional vector, which contains the first three formants sampled at
the steady state (whole audio signal) and three temporal sections, centered at
20%, 50% and 80% of the total vowel duration respectively. The subset of
vowels we consider are the vowels spoken by 37 female speakers which have
complete formant information (some vowels in the full dataset have formants
that could not be measured). Since each speaker utters 7 different vowels, the
vowel dataset contains a total of 259 vowels ($N_{\text{data}}=257$).

In order to accommodate the finite input range of the SHG experiment, we
normalize them via the element-wise min-max feature scaling, where each
dimension of the input vector is renormalized by their corresponding maximum
and minimum values (over all 259 vowels). Mathematically, if the original
dataset is given by
$\set{(\bm{x}^{(k)},\bm{y}^{(k)})}_{k=1}^{N_{\text{data}}}$, the examples of
the transformed dataset
$\set{({\bm{x}^{\prime}}^{(k)},\bm{y}^{(k)})}_{k=1}^{N_{\text{data}}}$ is
given by
${x^{\prime}}^{(k)}_{i}=\left(x^{(k)}_{i}-x_{\text{min},i}\right)/\left(x_{\text{max},i}-x_{\text{min},i}\right)$,
where the element-wise minimum and maximum values are given by
$x_{\text{min},i}=\min\limits_{k}x^{(k)}_{i}$ and
$x_{\text{max},i}=\max\limits_{k}x^{(k)}_{i}$.

The full architecture for performing the vowel classification with the SHG
system is schematically shown in Fig. S21. In the following, we outline the
steps involved in performing an inference with the PNN. To show how the
intermediate quantities evolve in the PNN, we plot quantities that are
labelled by letter labels in Fig. S21 in subpanels of Fig. S22.

  1. 1.

The 12-D normalized vowel formant frequencies are overall rescaled, i.e.
$x_{i}\rightarrow ax_{i}+b$, where $a$ and $b$ are trainable parameters.

  2. 2.

The 12-D vector is then repeated 4 times. The 48-D output is the input to the
first SHG device, i.e.
$\bm{x}^{[1]}=[x_{1},x_{1},x_{1},x_{1},\dots,x_{12},x_{12},x_{12},x_{12}]^{\mathrm{T}}$
to arrive at a 48-D input for the first SHG system. This 48-D input is shown
in A of Fig. S21 and Fig. S22.

  3. 3.

For the first SHG system, this input is concatenated with a 52-D parameter and
the total 100-D physical input
$\bm{q}^{[1]}=[\bm{x}^{[1]},\bm{\theta}^{[1]}]^{\mathrm{T}}$ is sent to the
first SHG.

  4. 4.

The SHG performs the physical computation
$\bm{s}^{[1]}=f_{\text{SHG}}(\bm{q}^{[1]})$. The output of the first SHG
$\bm{s}^{[1]}$ is shown in B of Fig. S21 and Fig. S22.

  5. 5.

The 50-D output of the first SHG $\bm{s}^{[1]}$ is renormalized and an overall
trainable rescaling is applied, i.e.
$\bm{y}^{[1]}=a^{[1]}\bm{s}^{[1]}/\max(\bm{s}^{[1]})+b^{[1]}$, where $a^{[1]}$
and $b^{[1]}$ are trainable scalars.

  6. 6.

The 50-D output is then concatenated with a 50-D parameter and sent to the
next SHG device.

  7. 7.

Step 4 and 5 are repeated until the output of the fifth SHG device has been
renormalized and rescaled. Mathematically, they are given by,

$\displaystyle\bm{x}^{[l]}$ $\displaystyle=y^{[l-1]}$ (24)
$\displaystyle\bm{q}^{[l]}$
$\displaystyle=[\bm{x}^{[l]},\bm{\theta}^{[l]}]^{\mathrm{T}}$ (25)
$\displaystyle\bm{s}^{[l]}$ $\displaystyle=f_{\text{SHG}}(\bm{q}^{[l]})$ (26)
$\displaystyle\bm{y}^{[l]}$
$\displaystyle=a^{[l]}\bm{s}^{[1]}/\max(\bm{s}^{[l]})+b^{[l]}$ (27)

for $l=2,3,\ldots 5$. The output of each SHG device is shown in C-E of Fig.
S21 and Fig. S22.

  8. 8.

From the 50-D output, a 21-D section is cropped as follows
$\bm{y}^{[\text{c}]}=[y^{[5]}_{12},y^{[5]}_{13},\ldots
y^{[5]}_{32}]^{\mathrm{T}}$. The 21-D output is shown in F of Fig. S21 and
Fig. S22.

  9. 9.

Each 3-D section of this new 21-D output is summed to arrive at the final
output of the PNN, i.e.
$\bm{o}=[y_{1}^{[\text{c}]}+y_{2}^{[\text{c}]}+y_{3}^{[\text{c}]},y_{4}^{[\text{c}]}+y_{5}^{[\text{c}]}+y_{6}^{[\text{c}]},\ldots
y_{19}^{[\text{c}]}+y_{20}^{[\text{c}]}+y_{21}^{[\text{c}]}]^{\mathrm{T}}$.
This final output is shown in G of Fig. S21 and Fig. S22.

Analogous to Sec. 1.4, the PNN is trained by minimizing the following loss
function:

$\displaystyle
L=H(\bm{y},\text{Softmax}(\bm{o}))+\lambda_{\mathrm{\theta}}\sum_{l}L_{\text{bound}}(\bm{\theta}^{[l]},v_{\text{min}},v_{\text{max}})+\lambda_{\mathrm{x}}\sum_{l}L_{\text{bound}}(\bm{x}^{[l]},v_{\text{min}},v_{\text{max}})$
(28)

where $H$ is the cross-entropy function and
$L_{\text{bound}}(\bm{z},v_{\text{min}},v_{\text{max}})=\sum_{l}\frac{1}{N_{\bm{z}}}\sum_{i}^{N_{\bm{z}}}\max(0,z_{i}-v_{\text{max}})-\min(0,z_{i}-v_{\text{min}})$.
The last two terms are regularization terms to constraint the inputs and
parameters of each SHG device, whose parameters are
$\lambda_{\mathrm{\theta}}=0.5,\lambda_{\mathrm{x}}=2,v_{\text{min}}=0,v_{\text{max}}=1$.

Instead of vanilla stochastic gradient descent (as in (19)), we perform PAT
with the Adadelta optimizer. This is because the adaptive learning rate is
helpful for dealing with the different nature of parameters that exist in the
model. The initial learning rate is $1.0$ and it was reduced by half every 700
epochs.

Figure S20: Inputs and outputs of the physical transformation performed by the
femtosecond SHG device. Here, 100 dimensional physical inputs applied to the
pulse shaper are equally split between the input (machine learning data) and
parameters. The outputs of the SHG is 50 dimensional. As shown by the good
agreement between the black and green curves, the data-driven models emulate
the physical system well. Note that the physical inputs considered here is
part of the “test” set and has not been used to train the data-driven model.
Figure S21: The full PNN architecture used for the second harmonic generation
vowel classifier (for the case of 5-layers). Annotation letters (A-H) refer to
the plots in Fig S22.  Figure S22: The sequence of operations that make up the
trained SHG PNN vowel model. The lettered labels corresponds to intermediate
quantities that are labelled in Fig. S21.

####  2.5.2 Figure 4: Analog transistor-MNIST PNN

As described in the main article, the electrical circuit chosen to serve as an
electronic PNN demonstration was selected to exhibit highly nonlinear
transient dynamics by embedding a transistor within an RLC oscillator.

This circuit was, unexpectedly, particularly challenging to produce an
accurate differentiable digital model for. While we expect that incorporating
knowledge-based components or other physical insight could improve the
results, in order to keep our procedure general, we found that an effective
workaround was to simply reduce the dimensionality of the input-output
transformation. Empirically, we found that 70-dimensions (70 input samples,
and 70 output samples) was small enough to achieve reasonable results.

The digital model was found by performing a neural architecture search. The
best model found obtained a validation loss of 0.09 (with data normalized to
-1 to 1). It is a 4-layer neural network with swish activation functions, and
217, 427, 167, and 197 hidden units. The Adam optimizer was used with a
learning rate of 4.56$\times 10^{-4}$.

A related challenge with the analog electrical PNN was that our implementation
exhibits high noise, mainly due to lost samples during streaming over the USB
connection, and to jitter related to triggering at the limit of the device’s
temporal resolution. Even after implementing post-selection to account for
traces corrupted or lost during the fast read-write operation of the DAQ, the
mean deviation amplitude between the output voltage time series for two
identical inputs was over 5% of the maximum absolute output amplitude, and
nearly 20% of the absolute mean output amplitude.

Figure S23: The full PNN architecture used for the analog transistor circuit
MNIST digit classifier. For details see main text.

Supplementary Figure S23 shows the differentiable digital model’s predictions
compared to the circuit’s actual output. As a reference point, the figure
shows 40 output traces from the circuit, each taken for the same inputs. The
figures illustrates that a main limiting factor in the digital model’s
accuracy is in fact the large, complex noise in the circuit’s input-output
transformation.

These two considerations, the low dimensionality of the input-output
transformation, and the significant noise of the circuit’s operation,
influenced our PNN architecture design, shown in Supplementary Figure S24.

To make the 196-D MNIST images compatible with the 70-D input-output
transformations of the electronic circuit, we chose to sum blocks of the MNIST
image into 14-dimensional vectors (i.e., perform a pooling operation on the
196-dimensional vector in blocks of 14), then repeat this 14-dimensional
vector 5 times as the input to the first layer of the PNN architecture. Of
course, many other options are possible, but we found this strategy to be
effective.

Figure S24: The full PNN architecture used for the analog transistor circuit
MNIST digit classifier. Annotation letters (A,B,C,D,…) refer to the plots in
Fig S25 and the description in the main text. Figure S25: The sequence of
operations that make up the trained electronic PNN MNIST model. For details of
each section, see main text. In part O, the different-colored curves
correspond to the output from the 7 different 3-layer PNNs.

To obtain robust performance despite the circuit’s noise and the digital
model’s relatively poor accuracy, we chose to implement two additional design
components. First, we operate the PNN as a ‘committee’: 7, 3-layer PNNs each
produce a prediction, and the full PNN produces its prediction by taking a
weighted average of these. Second, we added trainable skip connections between
layers of the PNN. These connections, inspired by residual neural networks,
were added to allow the network to act further as an ensemble of sub-networks
[11], thus adding to the redundancy of the architectures (and, in principle,
helping it to achieve good classification performance even when any one sub-
part stochastically fails) [11].

Each inference of the full architecture involves the following steps. By
letter labels (A,B,C, …), we refer each step to its corresponding location in
Supplementary Figure S24 and to the subpanels of Supplementary Figure S25.

  1. 1.

196-D MNIST images (A in Supplementary Figures S24 and S25, downsampled from
784-D) are element-wise rescaled, i.e. $x_{i}\rightarrow x_{i}a_{i}+b_{i}$,
where $a_{i}$ and $b_{i}$ are trainable parameters .

  2. 2.

Pool the 196-D vector from 196 to 14-D (i.e., sum each region of 14).

  3. 3.

Concatenate this 14-D vector 5 times to produce a 70-D input for the first
physical transformation (B in Supplementary Figures S24 and S25)

  4. 4.

Run the physical transformation once, producing an output 70-D vector
$\vec{y}$ (B-C in Supplementary Figures S24 and S25).

  5. 5.

Add the output $\vec{y}$ to the initial input via trained skip connection,
i.e. $\vec{y}\rightarrow\vec{y}+a\vec{x}$, where $a$ is a trainable skip
weight and $\vec{x}$ was the input to the physical system (D-E in
Supplementary Figures S24 and S25).

  6. 6.

Apply trainable, element-wise rescaling to produce the input for the second
physical transformation, i.e., $x_{i}=y_{i}a_{i}+b_{i}$ (E-F in Supplementary
Figures S24 and S25).

  7. 7.

Run the physical transformation once, producing an output 70-D vector
$\vec{y}$ (F-G in Supplementary Figures S24 and S25).

  8. 8.

Add the output $\vec{y}$ to the initial input via a trained skip connection,
i.e. $\vec{y}\rightarrow\vec{y}+a\vec{x}$, where $a$ is a trainable skip
weight and $\vec{x}$ was the input to the physical system (H-I in
Supplementary Figures S24 and S25).

  9. 9.

Apply trainable, element-wise rescaling to produce the input for the second
physical transformation, i.e., $x_{i}=y_{i}a_{i}+b_{i}$ (I-J in Supplementary
Figures S24 and S25).

  10. 10.

Run the physical transformation once, producing an output 70-D vector
$\vec{y}$ (J-K) in Supplementary Figures S24 and S25).

  11. 11.

Add the output $\vec{y}$ to the initial input via trained skip connection,
i.e. $\vec{y}\rightarrow\vec{y}+a\vec{x}$, where $a$ is a trainable skip
weight and $\vec{x}$ was the input to the physical system (L-M in
Supplementary Figures S24 and S25).

  12. 12.

Sum 50 elements of $\vec{y}$ in sections of 5 to produce a 10-dimensional
vector (M-N in Supplementary Figures S24 and S25).

  13. 13.

Repeat the above steps for 6 other sub-PNNs, each with its own independent
trainable parameters, producing its own independent output 10-dimensional
vector (Outputs of each independent sub-PNN shown in O in Supplementary Figure
S25).

  14. 14.

Take a weighted sum of the 7, 10-dimensional vectors,
$\vec{y}_{\text{out}}=\sum_{i=1}^{7}w_{i}\vec{y}_{i}+b_{i}$ (P in
Supplementary Figures S24 and S25).

To train the PNN, we used PAT with the Adadelta optimizer with learning rate
of 0.75, which was reduced to 0.375 after 30 epochs, then to 0.28 after 45
epochs. We found Adadelta to be particularly effective compared to other
optimizers, and hypothesize that momentum may inhibit the suppression of
gradient errors in PAT, at least in these very small-scale proof-of-concept
PNNs. Although we found it resulted in only moderate improvements in small-
scale tests, we retrained the digital model for the electronic circuit every
two training epochs, using a combination of random 70-dimensional inputs and
local data (i.e., the inputs and outputs obtained from the current training
data). This was done in part to help compensate for variations in the
circuit’s response over time, and also to nominally help improve the digital
model’s accuracy within the local region of parameter space.

####  2.5.3 Figure 4: Oscillating plate MNIST PNN

Figure S26: Multiple microphone recordings for 196 uniformly distributed input
voltages at an input frequency of $192$kHz and the digital models’
predictions. The linearity of the input-output transformation allowed
excellent agreement between digital model and experiment as evident from the
almost indistinguishable traces. Figure S27: The full PNN architecture used
for the oscillating plate MNIST digit classifier. Annotation letters
(A,B,C,D,E) refer to the plots in Fig S28 and the description in the main
text. Figure S28: The sequence of operations that make up the oscillating
plate PNN MNIST model. For details of each section, see main text.

The linearity of the oscillating plate’s input-output transformation allowed
us to accurately learn even very high-dimensional input-output
transformations. There was no need to downsample the original 784-D MNIST
images as we could learn a 78-4D to 784-D dimensional input-output map.
Instead, we could unroll and encode a whole image into a time-multiplexed
voltage signal. Outputs of the same dimension were read out from the voltage
over time signal of the microphone.

The strong performance of a linear single-layer perceptron on the MNIST
dataset and the realization that the oscillating plate PNN performs a matrix
multiplication as in Eq. 23 inspired us to use the linear oscillating plate in
a similar fashion. By cascading a few layers of element-wise-rescaling and
passing the signal through the speaker we roughly emulate the operation of a
single-layer perceptron.

The final layer of the oscillating plate consists of reading out a short time
window from the microphone output voltage. The exact position and length of
the time window was a trainable hyperparameter of our network architecture,
however, its position was generally closer to the end of the output as
information from the input can only interact with later inputs due to
causality. The window is divided into 10 bins. MNIST digits were classified
according to which bin had the strongest average signal.

Each inference of the full architecture involves the following steps:

  1. 1.

784-D MNIST images are element-wise rescaled, i.e. $x_{i}\rightarrow
x_{i}a_{i}+b_{i}$, where $a_{i}$ and $b_{i}$ are trainable parameters.

  2. 2.

Run the physical transformation once, producing an output 784-D vector
$\vec{y}$.

  3. 3.

Apply trainable, element-wise rescaling to produce the input for the second
784-D physical transformation, i.e. $x_{i}\rightarrow x_{i}a_{i}+b_{i}$.

  4. 4.

Run the physical transformation once, producing an 784-D vector $\vec{y}$.

  5. 5.

Apply trainable, element-wise rescaling to produce the input for the third
784-D physical transformation, i.e. $x_{i}\rightarrow x_{i}a_{i}+b_{i}$.

  6. 6.

Run the physical transformation once, producing an 784-D vector $\vec{y}$.

  7. 7.

Take outputs 724 to 774 and average 5 consecutive points to produce a 10-D
vector $\vec{y_{\text{out}}}$.

Due to the simple PNN architecture and the very accurate linear digital twins,
we could achieve high accuracy on MNIST by only using _in silico_ training and
transferring the trained parameters onto the experiment. Using PAT still
resulted in performance gains, although not as dramatic as for the nonlinear
physical systems. Usually, _in silico_ performance in experiment only fell a
few percentage points short of PAT, and after just a few epochs ( 3) of
transfer learning with PAT, the full performance could be retrieved. We expect
therefore that, at least for very shallow networks, PAT is less crucial for
linear physical systems, and we expect that in silico training can be
generally applied to greater effect in such systems. While this may be useful
in some contexts, we note that sequences of linear systems cannot realize the
computations performed by deep neural networks, since they ultimately amount
to a single matrix-vector multiplication.

####  2.5.4 Figure 4: Femtosecond SHG-MNIST PNN

Figure S29: The full PNN architecture used for the SHG MNIST digit classifier.
Annotation letters (A,B,C,D,E) refer to the plots in Fig S30 and the
description in the main text. Figure S30: The sequence of operations that make
up the SHG PNN MNIST model. For details of each section, see main text. In
parts E and F, the grey section shows the region clipped. In part H, the
partially-transparent bars show the location of the correct classification.
The different color curves (blue, orange, green) correspond to 3 different
input images.

For the hybrid physical-digital MNIST PNN we present in Fig. 4, we chose to
use linear input layers to pre-process MNIST images prior to SHG
transformations. The structure of the input layer is in two segments, such
that the input dimension is reduced to 50, then transformed back to 196, which
is the input-output dimension of the SHG transformation. This seemingly
unusual choice was made because we observed that the SHG tended to produce
very little change in its output spectrum for high-dimensional random input
vectors. Thus, we view the two-step input layer as projecting the input MNIST
image into a low-dimensional space of SHG ‘input modes’, which ideally
produces highly variable output spectra. We found that, in order to maximize
the sensitivity of the SHG’s output spectra to the shape of the input
modulation, it was optimal to use an effective input dimension of about 10-30
“input modes” (i.e., to choose the input layer to be comprised of two
segments, with a middle dimension of 10-30).

In order to minimize the size of the SHG-MNIST model considered in the main
article however, we ultimately increased this input layer dimension to a
relatively large value (50) so that good performance could be achieved with
relatively few uses of the SHG process. Our desire to keep the model size
small (in terms of number of physical system uses) was driven by the
relatively slow (200 Hz) speed of the SHG transformation - if the SHG
transformation were used many times, it would have been difficult for us to
test different configurations. The result of these choices is, however, that
the SHG plays less of a role in the hybrid network’s calculations than it can
with a more optimal choice of input dimension.

Despite this deviations from a more optimal model, the hybrid model presented
in Fig. 4 nonetheless leverages the SHG’s nonlinear transformations to
effectively gain $\sim$ 7% improvement on the MNIST digit classification task
compared to when the SHG transformations are replaced by identity
transformations. While this seems like a small improvement, we note that the
difficulty of the MNIST digit classification task is highly nonlinear in
accuracy: a model reaching nearly 90% needs only $\sim 10^{3}$ operations,
while one reaching 96% or higher may require 100-1000 times as many
operations.

Similar to the electronic PNN, we found it was helpful to include skip
connections in the network. This is especially helpful because the SHG
transformation does not support an identity operation, and is therefore prone
to losing information.

The PNN architecture used for MNIST digit classification with the SHG physical
system in Fig. 4 of the main paper is depicted schematically in Fig. S29. The
flow of inputs through the network is shown in Fig, S30, where 3 different
inputs are shown by different colored-curves. In all three cases, the digits
are classified correctly.

Each inference of the full architecture involves the following steps:

  1. 1.

196-D MNIST images are operated on by a trained 196x50 matrix, producing a
50-dimensional vector. Without any nonlinear activation, a second trained
50x196 matrix is applied to produce the input to the first SHG transformation
(Fig. S30A). This is performed for 4 separate, independent channels (shown in
different columns of Fig. S30).

  2. 2.

The SHG apparatus is executed with the 196-dimensional input, producing a
196-dimensional output vector, which is depicted in Fig. S30B.

  3. 3.

The output of the SHG transformation is operated on by a trained element-wise
rescaling, i.e. $x_{i}\rightarrow x_{i}a_{i}+b_{i}$, resulting in Fig. S30C.

  4. 4.

To this is added the input to the SHG transformation, multiplied by a trained
(1-D) skip weight. This results in Fig. S30D.

  5. 5.

The SHG apparatus is executed with the 196-dimensional input, producing a
196-dimensional output vector, which is depicted in Fig. S30E.

  6. 6.

The SHG transformation’s output is clipped from index 70 to 140, and added to
the SHG input, which is multiplied by a trained skip weight. This produces a
70-dimensional vector depicted in Fig. S30F.

  7. 7.

The 4 independent channels’ output 70-dimensional vectors are concatenated
into a 280-dimensional vector, shown in Fig. S30G.

  8. 8.

The 280-dimensional vector is summed in bins of length 28 to produce a
10-dimensional vector $y_{\text{out}}$.

  9. 9.

Perform softmax with a trainable temperature $T$ on $y_{\text{out}}$ (for
training). Perform argmax to obtain the predicted digit (for inference).

##  3 Simulation analysis of candidate PNN systems

In this section, we analyze the simulations of several candidate PNN systems.
The purpose of this section is three-fold:

  1. 1.

To identify, based on the form of the numerical simulations, what types of
mathematical operations candidate physical systems approximate,

  2. 2.

To determine plausible self-simulation advantages, which set an upper-bound on
the computational benefits each PNN system could provide, and

  3. 3.

To analyze the scaling of the latter with respect to physical parameters.

Physical processes closely resemble the structure of DNNs that operate on
natural data: real physical processes typically feature low-order
nonlinearities, hierarchical structure, local correlations, sparsity, and
symmetry [12]. DNNs that perform well on natural data either exploit these
features in their design (as in convolutional DNNs, which intrinsically
account for translational symmetry and locality, or in networks trained with
data augmentation to improve robustness), or they acquire it as the result of
extended training with natural data. Although DNNs are usually nominally
deterministic algorithms, noise like that encountered in real physical systems
used for PNNs is often added to improve their performance. DNNs are often
trained with data that is distorted or has had noise added to improve the
generalization of their models or to reduce sensitivity to adversarial
examples. During training, adding noise (such as drop-out) is also common for
these purposes. In other applications, such as generative models, noise can
play a more fundamental role. Noise within PNNs may be useful for similar
reasons, provided it is not too large. In sum, all the reasons DNNs work well
for computations involving natural data and natural processes are
correspondingly strong rationale for why PNNs should be able to perform many
of the same tasks modern DNNs perform.

Physical systems that are good candidates for PNNs implement operations that
are well-suited to DNNs, even though they may not always be capable of
universal computation. Such candidate PNN processes implement controllable
matrix-vector operations and convolutions, as well as nonlinearities and
tensor contractions. In the vast majority of physical systems, arbitrary
control over these physically-implemented computations will not be easy.
Instead, the physical processes will provide these operations approximately
and controllable within some subspace of possible operations (typically
subject to sparsity, locality and other symmetries) rather than fully-
arbitrary mathematical operations. As a result, although these physical
systems may provide significant computational benefits via operations similar
to those implemented by modern DNNs, they are unlikely to be useful in
isolation, and will sometimes require cooperation from a universal digital
computer. In this case, it becomes essential to keep track of how many
operations the PNN is saving the digital computer - this and related
considerations are addressed in the last section of this document.

We find that, generally, self-simulation advantages are higher when the
physical system is intuitively more complex. Concretely, this corresponds to
it being designed such that:

  1. 1.

Nonlinearity is present, especially higher-order nonlinearities and those that
couple multiple degrees of freedom.

  2. 2.

All degrees of freedom are coupled in heterogeneous ways to many other degrees
of freedom.

  3. 3.

The dimensionality is high (e.g. 3-D versus 2-D). However, this will often
correspond to a reduced nonlinearity.

We find that self-simulation advantages significantly over $10^{6}$ are common
in candidate systems. This is comparable to the increase in speed of
supercomputers from 1995 to the present (from roughly 1011 to 1017.5 GFLOPs
[13]. Of course, the actual potential of PNNs depends both on how well the
self-simulation advantage estimates computational advantages for other, more
useful tasks, how well PNN models can be designed to leverage advantageous
physical computations, and how the physical devices used as PNNs might advance
in the coming years. Nonetheless, such large self-simulation advantages (SSAs)
confirm that, if these non-trivial steps can be completed successfully, non-
traditional physical systems do have enormous potential for advancing real-
world machine learning applications within the next decades.

###  3.1 General case

The findings summarized above can be appreciated intuitively by considering a
coupled nonlinear differential equation model of a physical dynamical system
with $N$ distinct degrees-of-freedom (modes) $x_{i}$, which interact with up
to $n$th-order nonlinearities, and are driven by $N$ driving forces $f_{i}(t)$

$\dot{x}_{i}(t)=\sum_{j=1}^{N}a_{ij}(t)x_{j}(t)+\sum_{j=1}^{N}\sum_{k=1}^{N}b_{ijk}(t)x_{j}(t)x_{k}(t)+\sum_{j=1}^{N}\sum_{k=1}^{N}\sum_{l=1}^{N}c_{ijkl}(t)x_{j}(t)x_{k}(t)x_{l}(t)+...+f_{i}(t).$
(29)

For such a system, if the coupling of the $N$ degrees of freedom is
heterogeneous and non-negligible over $M$ neighbours for all orders of
nonlinearity, then the simulation complexity is dominated by the highest-order
nonlinear term as $\Delta T\Delta fNM^{n}$, where $\Delta f$ is the bandwidth
of the system (the inverse of the necessary time step) and $\Delta T$ is the
simulated time. Note that, because sparsity, locality and other symmetries can
be exploited to make each term faster to numerically evaluate, it is important
to consider the coupling range $M$ of heterogeneous, non-negligible coupling.
In most systems, the coupling range depends on the order of the nonlinearity,
so a more general expectation is

$\text{Simulation Complexity}\sim\Delta T\Delta fNQ_{\text{max}},$ (30)

where $Q_{\text{max}}=\text{max}[M_{i}^{i}]$ is the maximum coupling
complexity over the different terms, where $M_{i}$ is the coupling range of
heterogeneous coupling for the $i$th order nonlinear term.

The power required to sustain such a physical system will usually be linear in
the number of modes $N$, and the time taken is of course $\Delta T$, so the
self-simulation advantage will typically scale with $Q_{\text{max}}$. This
however assumes that $n$ and $M$ are independent of the power supplied to the
physical system, which will generally not be true: in most systems, as the
energy localized in the system is increased, more nonlinearity and more
complex, long-range coupling are necessary to describe the physics. As the
dimension of the dynamics increases, the ability for localized degrees of
freedom radiating a fixed amount of energy to couple with neighbours increases
exponentially. This is intuitive: higher-dimensional physics is harder to
simulate.

Equation 29 also identifies, very generically, some of the physical data
transformations realized by physical systems. If the physical parameters of
the system that determine $a_{ij}$, $b_{ijk}$, $c_{ijkl}$ can be adjusted and
trained, then the evolution of the $x_{i}(t)$ implement matrix-vector
operations (first term) and higher-order tensor contractions (subsequent
terms). Depending on the locality and symmetry of $a_{ij}$, $b_{ijk}$,
$c_{ijkl}$, these operations will often be similar to convolutions. In
general, these coefficients could also be controlled in time, $a_{ij}(t)$,
$b_{ijk}(t)$, $c_{ijkl}(t)$. While this naturally opens up significantly more
complex sequences of operations in a single physical system, it may not be
compatible with in-place parameters (parameters which are stored in-place in
the physical device, and therefore do not need to be retrieved from a digital
memory each time an inference is performed). For this reason (see the last
section of this document), we expect the most useful, scalable
parameterization of most physical systems to be the use of trainable, but
time-independent physical parameters $a_{ij}$, $b_{ijk}$, $c_{ijkl}$.

Alternatively, if some of the $x_{i}(t)$ are controllable parameters (i.e.,
are treated as forced source terms, rather than terms obeying Eqn, 29), then
one may recognize the nonlinear terms as realizing tensor operations of order
$n$, where $n$ is the order of nonlinearity (for example, the quadratic term
$\sum_{j=1}^{N}\sum_{k=1}^{N}b_{ijk}x_{j}x_{k}$ could result (in part) in a
controlled matrix-vector operation if some of the $x_{k}(t)$ are controlled
parameters). If the drive terms $f_{i}(t)$ are used as controllable
parameters, the outcome will probably be similar to controlling $x_{i}(t)$,
depending on the memory of the system’s dynamics, and the strength of the
drive. For the same reasons mentioned in the last paragraph, such time-
dependent physical parameters may pose scaling challenges if they require fast
reading from a digital memory, and digital-to-analog conversion.

This general model also suggests that the greatest physical advantages can be
gained with operations that are typically not used in modern deep neural
networks, namely high-dimensional tensor operations, since these are
associated with higher-order nonlinear coupling. Of course, this depends on
utilizing these operations efficiently to perform machine learning, which has
not been widely studied. Understanding how to use operations of this kind
efficiently for machine learning could thus help to advance the computational
benefits obtained with PNNs.

In the rest of this section as follows, we will analyze the physical costs
(energy and time) to evolve a selection of realistic PNN modules based on
different physical systems, then compare these physical costs to those
incurred by performing a simulation with a state-of-the-art digital electronic
device. For these comparisons, we will consider the Nvidia DGX Superpod
supercomputer [14] which achieves 39 pJ per floating point operation and
2.8$\times 10^{15}$ operations/s, and the Nvidia Titan RTX GPU [15], which
achieves 9.8 pJ/operation and 32$\times 10^{12}$ operations/s for half-
precision operations.

###  3.2 Single-mode nonlinear optical pulse propagation: scaling

For ultrafast pulse propagation in nonlinear dispersive media, a standard
simulation algorithm is the split-step Fourier method [16]. We will use the
computational complexity of this algorithm as our basis for assessing self-
simulation advantage.

For an optical propagation length $L$, the split-step algorithm divides the
propagation into small steps whose length must be much shorter than the rate
at which nonlinear effects occur, $L_{\text{nl}}$. The pulse is described by a
vector of $N_{\omega}$ values which contain the complex amplitudes of the
electric field envelope. The number of spectral points is determined by the
maximum frequency range of the pulse and any new frequency content the
propagation creates, and by the time window needed to describe the pulse
throughout the propagation, $N_{\omega}=T\Delta\omega$, where $\Delta\omega$
is the bandwidth of the simulation and $T$ is the required time window of the
simulation.

In the split-step algorithm, linear effects like dispersion are efficiently
applied in the Fourier domain. The complexity of the split-step algorithm is
dominated by the Fourier transform operations required, so:

$N_{\text{ops}}\sim N_{z}N_{\omega}\log(N_{\omega}),$ (31)

where $N_{z}$ is the number of longitudinal steps required. The minimum step
size is proportional to the rate of nonlinear effects, $N_{z}=L\times\gamma
f(P)$, where $f(P)$ is the dependence of the nonlinear rate on the peak power
$P$ of the pulse and $\gamma$ is the nonlinear constant, which depends on the
material, the type of nonlinearity, and the transverse confinement of the
single spatial mode the pulse propagates in. For Kerr (cubic) nonlinearity,
$f(P)=P$, whereas for quadratic nonlinearity $f(P)=\sqrt{P}$, generally
$f(P)=P^{(n-1)/2}$ where $n$ is the order of the nonlinearity.

The peak power $P$ is related to the minimum temporal feature of the pulse,
$\sim 1/(\Delta\omega)$, where $\Delta\omega$ is the bandwidth of the
simulation. $N_{\omega}=T\Delta\omega$, where $T$ is the required time window
of the simulation. In other words, $N_{\omega}$ is the maximum time-bandwidth
product that occurs during the propagation, which generally increases for more
complex, dispersive and nonlinear pulse propagation. For the sake of
simplifying the discussion, we will let $P_{0}=E_{p}N_{\omega}$, though we
note that this will often overestimate the nonlinearity.

The maximum rate at which physical computations can be performed is determined
by the longest of the following: duration between input light pulses, duration
between electro-optic modulations, and the maximum duration of the optical
pulses during propagation, with the last required to ensure that each physical
computation is independent (which is usually, but not always required for
classification tasks). We will assume this corresponds to about 33 ps (so a
rate of 30 GHz). Pulsed sources of this repetition rate are possible [17], and
electro-optic modulators with bandwidths in excess of 40 GHz are standard.
Provided the length of the nonlinear optical medium, and the dispersion and
pulse spectral broadening are chosen so that the maximum pulse duration never
exceeds 33 ps, 30 GHz operation should therefore be feasible.

The physical advantage for computation (PAC) in time is then:

$\text{PAC}_{t}=\frac{\text{Time for simulation}}{\text{Time for physical
evolution}}=\frac{t_{c}N_{\text{ops}}}{t_{p}},$ (32)

where $t_{c}$ is the time per operation for the reference conventional
computer, and $t_{p}$ is the time the physical evolution takes (being limited
by one of the above considerations).

Similarly, the physical advantage for computation in energy is:

$\text{PAC}_{e}=\frac{\text{Energy for simulation}}{\text{Energy for physical
evolution}}=\frac{e_{c}N_{\text{ops}}}{E_{\text{phys}}},$ (33)

where $e_{c}$ is the energy per operation of the reference computer and
$E_{\text{phys}}$ is the total energy required to instantiate and measure the
pulse propagation. Instantiating the pulse propagation and measuring its
result requires significantly more energy than just that needed for the
optical pulse, as an example will show.

####  3.2.1 Quantitative example

In this section, we calculate the energy and speed self-simulation advantage
of a realistic, scaled device based on ultrafast nonlinear pulse propagation.
Our intention here is not to design a device that would maximize the self-
simulation advantage within the most optimistic conditions, but to consider a
realistic (though not trivial) experimental demonstration possible with
existing technology. Throughout, we will note room for improvement in near-
future devices, and at the end analyze both the realistic present-case and
near-future potential.

We consider a device as follows. An ultrashort pulse at 1550-nm is split into
20 wavelength bins, where it can be modulated by electro-optic amplitude
modulators and then recombined. Through these electro-optic modulators, input
vectors can be applied at $\sim$ 10-100 GHz via spectral modulations. For
physical parameters, the pulse is then shaped by a phase pulse shaper based on
a spatial light modulator or similar device. This device allows $\sim$ 103 or
more independent controllable parameters for the subsequent nonlinear pulse
propagation, which can be stored in-place over many uses of the PNN without
incurring energy costs due to memory reads or digital-analog conversion. The
modulated pulse then propagates in a nanophotonic waveguide based on
periodically-poled thin film lithium niobate. To maximize the complexity of
the dynamics, we assume the waveguide is poled to create an artificial Kerr
effect.

In Ref. [18], 1 pJ pulses were sufficient to observe significant
supercontinuum pulse propagation in a cm-scale waveguide. However, the pulses
used were simple, transform-limited pulses. To ensure very complex pulse
propagation with the structured, data-encoded pulses used for computation, we
assume that roughly 50 times as much pulse energy should be used (spread
across a longer, more complex pulse). To account for various optical losses in
a realistic device (transmission around 1%), and the wall-plug efficiency of
the initial mode-locked laser (around 10%), we assume a total of 50 nJ of
energy is required per optical pulse. This energy could be reduced by several
orders of magnitude in a fully-integrated device, since the optical losses
would be drastically reduced. Assuming a 0.5 cm waveguide, the pulse duration
should remain below 33 ps throughout the waveguide (for a dispersion of -15
fs2/mm, such a waveguide would correspond to at most 10-20 dispersion lengths,
so a maximum duration of 10-20 ps is reasonable to expect given the expected
pulse bandwidths described below).

To inject data to the physical system, it must be converted from a digital
computer’s digital format to a format (typically analog) suitable for the
physical system. The most energy-efficient ADCs are capable of approximately 1
pJ/conversion [19] at 24 GHz. For the sake of this order-of-magnitude
estimate, we will assume ADC and DAC operations each require 1 pJ/conversion.
Integrated electro-optic modulators that operate with $\sim$ pJ/modulation are
typical [20, 21, 22], while modulators capable of few-bit-depth modulation
appear possible [22, 23, 24]. For our purposes here, we will assume each fast
optical modulation requires an additional 1 pJ on top of the DAC cost.
Although on-chip photodetectors with ultralow capacitance may enable
amplifier-free detection of light pulses as weak as 1 fJ [25, 26], we will
assume that detectors add an additional 1 pJ on top of the ADC cost for a
modern experiment. Finally, the spatial light modulator requires roughly 10 W
to maintain, resulting in 330 pJ per 33 ps inference. Of course, in future
devices this would be significantly reduced, presumably to significantly below
1 W.

Adding these physical cost per inference is:

$330\text{ pJ}+20\times 1\text{ pJ}+20\times 1\text{ pJ}+20\times 1\text{
pJ}+50\text{ nJ}=50.39\text{ nJ}$ (34)

for the present-day experiment, and

$33\text{ pJ}+20\times 0.1\text{ pJ}+20\times 0.1\text{ pJ}+20\times 10\text{
fJ}+0.5\text{ nJ}\approx 537\text{ pJ}$ (35)

for an improved, near-future device. For both devices, the time-per-
computation would be limited by the propagation time, 33 ps, since fast
electro optic modulators and DAC/ADC would permit input/output rates exceeding
1/(33 ps) = 30 GHz.

We now consider the energy and time costs for simulating this system with
state-of-the-art conventional hardware. For the 50 pJ pulses we have assumed,
if the waveguide is designed to support weak anomalous dispersion, peak powers
of roughly 1 kW ($\sim$ 10 pJ localized into a 10 fs region) are probable,
corresponding to a nonlinear length scale of
$10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ for the device in Ref. [18] which
has an effective self-phase modulation coefficient of 100/Wm. To adequately
simulate this would require a longitudinal step size of around
$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ or smaller. To accommodate spectral
broadening up to around 400 THz, a spectral range $\Delta\omega$ of 800 THz is
required for accurate simulations. To accommodate the complex shaped pulses
and dispersion of the broadband light, a temporal window size of 10-50 ps
would be expected, since the initial pulse duration would be roughly 1-5 ps
depending on the spectral modulations applied (100-1000 phase degrees of
freedom across a $\sim$ 50 THz input bandwidth), and the 5 mm waveguide
corresponds to roughly 5-8 characteristic dispersion lengths for the broadened
pulse (assuming a dispersion of -15 fs2/mm).

Thus, using the Fourier split-step method, the number of operations for
simulations would be:

$N_{z}N_{\omega}\log_{2}(N_{\omega})=(5\text{
mm}/$1\text{\,}\mathrm{\SIUnitSymbolMicro m}$)\times 800\text{ THz}\times
25\text{ ps}\times\log_{2}(800\text{ THz}\times 25\text{ ps})\approx 5\times
10^{8}\text{operations.}$ (36)

Using the Nvidia Titan RTX, this would cost 5 mJ and take
$16\text{\,}\mathrm{\SIUnitSymbolMicro s}$. This corresponds to roughly 105
times more energy and 105 times longer than the present-day experiment, and
roughly 107 more energy and 105 times longer than the near-term device.

Using the Nvidia DGX SuperPOD, the simulation could be performed more quickly,
179 ns, but with greater energy cost, 20 mJ. This corresponds to roughly 105
times more energy and 104 times longer than the present-day experiment, and
roughly 108 more energy and 104 times longer than the near-term device.

Although the self-simulation advantage possible with nanophotonic nonlinear
optics is large, the simulation analysis presented here highlights some
important caveats. First, the scaling of the self-simulation advantage for
single-mode nonlinear optical pulse propagation is in general only marginally
superlinear with respect to the optical degrees of freedom, $N_{\omega}$,
since the logarithmic factor is relatively weak. Since larger $N_{\omega}$ can
give rise to shorter pulses, in some regimes $N_{z}$ could scale
optimistically as $N_{\omega}$, suggesting a best-case scaling of
$N_{\omega}^{2}\log_{2}(N_{\omega})$. It is not obvious if pulse propagation
that supports this complexity scaling will also be useful for practical
computations, however. A more crucial caveat is that the operations performed
by single-mode nonlinear optical pulse propagation are, by default, not well-
matched to the standard operations performed in modern deep neural networks,