{'fraction_non_alphanumeric': 0.03145584083882242, 'fraction_numerical': 0.013935564816059506, 'mean_word_length': 4.802912116484659, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computing coherent light scattering on the millimetre scale using a recurrent neural network without training

Heterogeneous materials such as biological tissue scatter light in random, yet deterministic, ways. Wavefront shaping can reverse the eﬀects of scattering to enable deep-tissue microscopy. Such methods require either invasive access to the internal ﬁeld or the computational solv- ing of an inverse problem. However, calculating the coherent ﬁeld on a scale relevant to microscopy remains excessively demanding for consumer hardware. Here we show how a re- current neural network can mirror Maxwell’s equations without training. By harnessing public machine learning infrastructure, the light-ﬁeld throughout a 6mm 2 area or 110 3 µ m 3 volume can be calculated in 16 minutes. The elimination of the training phase cuts the calculation time and, importantly, it ensures a fully deterministic solution, free of training bias. We integrated our method with an open-source electromagnetic solver. This enables any researcher with an internet connection to calculate complex light-scattering in volumes that are larger by two orders of magnitude.

Heterogeneous materials such as biological tissue scatter light in random, yet deterministic, ways. Wavefront shaping can reverse the effects of scattering to enable deep-tissue microscopy. Such methods require either invasive access to the internal field or the computational solving of an inverse problem. However, calculating the coherent field on a scale relevant to microscopy remains excessively demanding for consumer hardware. Here we show how a recurrent neural network can mirror Maxwell's equations without training. By harnessing public machine learning infrastructure, the light-field throughout a 6 mm 2 area or 110 3 µm 3 volume can be calculated in 16 minutes. The elimination of the training phase cuts the calculation time and, importantly, it ensures a fully deterministic solution, free of training bias. We integrated our method with an open-source electromagnetic solver. This enables any researcher with an internet connection to calculate complex light-scattering in volumes that are larger by two orders of magnitude.
Machine learning has pushed automation into areas that were formerly considered the exclusive remit of humans. Recently it was demonstrated that artificial neural networks can also help us circumvent the curse of dimensionality in traditionally hard scientific computations 1 . Problems that scale poorly and were once thought to be intractable, are now within reach of industrial, yet publicly available, machine learning infrastructure.
Advances in optics and photonics increasingly rely on our ability to accurately compute how light propagates and scatters as described by Maxwell's Equations 2,3 . Coherent calculations provide essential information that can be impractical or impossible to obtain experimentally. Notwithstanding, finding a numerical solution to Maxwell's equations is challenging in its own right. A volume of relevance to microscopy can span hundreds of wavelength per dimension, while tens of samples per wavelength are demanded for most algorithms to keep error accumulation and numerical dispersion in check 4 . The problem's solution alone may thus take up a significant proportion of the computer memory with its billions of free parameters, not unlike deep learning models such as OpenAI's GTP-3 5 .
In what follows, we show how Maxwell's equations can be rephrased as a recurrent neural network that does not require training. This eliminates a large computational cost while ensuring that the solutions are fully deterministic and generalisable. The pre-defined neural network allows us to effectively harness the economy of scale of cloud-based machine learning infrastructure. We demonstrate its potential by calculating the lightscattering of a guide-star embedded deep within a millimetrescale heterogeneous structure and use the exit field for refocusing. Next, we show how the highly efficient cloud-based calculations make it possible to calculate the complete scattering and deposition matrices. Finally, we analyse and compare the performance on a desktop, a GPU workstation, and on Google Colaboratory; respectively improving the efficiency 20 and 100-fold.

Results
Maxwell's equations as a neural network We start this section by describing a two-layer forward neural network that requires slow training (Fig. 1a). Next, we incorporate an extra layer for preconditioning that accelerates training (Fig. 1b). This network is transformed into a recurrent neural network that does not require training (Fig. 1c).
A large neural network model is typically trained to produce target output values, t, for a large number of test cases. Training is the optimisation of parameters, p, to minimise the difference between the neural network's output and the training target. This can be written as where the function N (p) computes the network's outputs for the corresponding training inputs. Such parameters typically represent the weights and biases of thousands of neurons and millions of connections 7,8 . Maxwell's equations can be rewritten in a similar form.
For coherent light with an angular frequency, ω, Maxwell's equations for an inhomogeneous material can be written as the time-independent complex vector functions, ∇ × H(r) = j(r) + ∂D(r) ∂t = j(r) − i ωD(r).
Provided with a source of electromagnetic radiation, j(r), we aim to determine the electric vector field, E(r). The electric displacement vector field, D(r), the magnetic flux density, B(r), and the magnetising field, H(r), are related by the constituent relations. Without loss of generality, we consider D(r) = 0 r (r)E(r) and B(r) = µ 0 H(r), where only the relative permittivity, r (r) = n 2 (r), is spatially-variant. More general relations can account for anisotropy, chiral, and magnetic properties 9 . In a dielectric, Eq. (2) and Eq. (3) can be combined in the vector Helmholtz equation where k 0 = ω 0 µ 0 . In principle, this partial differential equation can be solved numerically by discretising space and determining the field, E, for a source, j, on a sufficiently dense sampling grid. Without the position-dependency, Eq. (4) can be written succinctly as The matrix, M , is generally too large to be represented directly in computer memory. However, the result of its multiplication with any vector can be computed efficiently. It calculation consists of two terms. The former is a convolution. The latter term can be seen to be a diagonal matrix with the values of k 2 0 r at the sample points. We calculate E = M −1 j by solving the minimisation problem for the parameters, p, The parallel with Eq. (1) is apparent. The neural network's response, N (p), is replaced by the electromagnetic response, M , and the training data, t, is replaced by the electromagnetic source current-density, j. However, there are notable differences. The neural network function, N (p), is non-linear, and in general, any solution with sufficiently low loss is adequate for machine learning applications. In contrast, Maxwell's equations are linear and typically only permit a unique solution. It is thus important that the optimisation of Eq. (6) deterministically converges to the true solution. Fig. 1a shows how Maxwell's equations can be equated with a particular neural network with parallel convolutional and multiplication layers. Although the neurons are represented as a flat layer, the activity of each input neuron corresponds to the electric field, E, at a point in 3D space and for each polarisation. As in the hidden layer of d'Arco et al., the activity of the corresponding neurons in the output layer is proportional to the 3D currentdensity source, j, of the electric field. Determining the scattered light field, E, thus corresponds to training, the neural network, i.e. optimising its parameters, so that it outputs the source currentdensity, j, of the waves. This can be done using standard training algorithms such as stochastic gradient descent (SGD) and adaptive moment estimation (Adam) 11 . While it converged to the true solution for small problems, we found that training such neural network was impractically slow. A 16λ × 16λ size system could take up to a minute to converge on an average desktop computer.
To make the neural network more responsive to training, we alter its layers. Eq. (5) can be left-multiplied by Γ −1 , the inverse of a non-singular preconditioner 12 , leading to the equivalent expression The modified potential, V := −1 − k 2 0 ( r − 0 ) / i ωµ 0 s , is determined from the material property distribution. The background permittivity 0 ∈ C can be chosen to minimise max r | r (r) − 0 | for optimal convergence. The complex scaling constant, s, must be chosen so that 1 + V < 1 and ℜ p, s −1 M p ≥ 0 for all p. This is possible for any gain-free system 12 .
The inverse preconditioner can also be written as Γ −1 = −s −1 V G, using the spatially-invariant Green's function, G := This can be implemented as a convolutional neural network layer that performs the inversion G := Eq. (7) can now be rewritten as V (1 +GV ) E = −Γ −1 j, and solved by training the neural network in Fig. 1b. Although, compared to that of Fig. 1a, this network has an extra layer, the use of preconditioning reduced the required training time by approximately 8-fold. More importantly, the preconditioning enables us to write the solution to Eq. (5) and Eq. (7) as the Neumann power series This maps directly to the recurrent neural network depicted in Fig. 1c. This neural network does not require training. When the current density distribution of a light source is provided as recurrent input, it infers the electric field distribution, E. Its accuracy can be readily verified by inserting the solution in the forward problem Helmholtz equation, represented by the network of Fig. 1a.
We implemented the recurrent method in PyTorch and integrated it into the electromagnetic solver MacroMax 13 . Its previously CPU-bound calculations can now be seamlessly executed on top-end machine-learning infrastructure. We used the publicly available Google Colaboratory to directly study coherent light scattering on scales relevant to microscopy.
Coherent optical scattering on the millimetre scale. The high efficiency and minimal memory requirements of the recurrent algorithm opens a door to calculations in large heterogeneous materials. To demonstrate this, we determine the field emitted from a visible light guide-star embedded deep within a heterogeneous scatterer. Next, we refocus on the source using phase-conjugation and analyse the effectiveness of the phase conjugation. Fig. 1f shows the emission over a 2.5 mm×2.5 mm-area on both sides of the 2.5 mm × 1.25 mm scattering material.
To refocus on a guide-star embedded within a multi-millimetre sized scattering medium we first record the wavefront on the left boundary of (Fig. 1f). Next, we conjugate phase of this wavefront and use it as a source. The wavefront can be seen to focus back onto the guide star in (Fig. 1g). It took only 16 min to calculate the field in the displayed 2.5 mm×2.5 mm system, running the algorithm on freely available Google Colaboratory 16 . It is equally possible to simulate 3D spaces with an equivalent number of sample points (i.e. in this case it would be a system of approximately 110 µm × 110 µm × 110 µm size). To emphasise the magnitude of the problem, it is worth noting that the data displayed in Fig. 1f-g alone contains 2.2 GB of information.
Calculating all possible internal fields. With the significantly improved calculation efficiency, we were able to solve problems that were previously deemed impractically large to address numerically. An example that fulfils this criterion is calculating all  14 . Vectors corresponding to the extremes of this graph were used as input modes for graphs a-b. (f) Complete scattering matrix with 1022 modes in total. Coherent backscattering can be observed to manifest in reflection matrices (located in bottom-left and top-right quadrants). This corresponds to the increased field amplitude on the diagonals of those regions 15 the internal fields for a scattering matrix. The scattering matrix is a mathematical construct that holds complete information on all the ways the light can enter and exit the system 14 . While the mathematical formalism for the scattering matrix was derived in the early days of quantum mechanics 17 , it did not see usage in adaptive optics until recently. Furthermore, instead of determining full scattering matrices most work focused on transmission and reflection matrices [18][19][20] , as these are more feasible to measure experimentally. Unfortunately, physical measurements of these matrices are limited for analysis purposes, as it is practically impossible to capture all light modes in an experimental setup 21 . But luckily, numerical field calculations can model complete and lossless data acquisition. While, recently new methods were proposed to calculate the entire scattering matrix rapidly 22 , calculating the internal fields of scattering matrix systems remains impractically slow.
We calculated scattering matrix using visible light for 128 µm× 128 µm heterogeneous system, with all the internal fields, in just 42 min on Google Colaboratory (Fig. 2). We use this scatter-ing matrix (Fig. 2f) to calculate open and closed channels of the system 14 , which are light modes that are either completely transmitted or blocked by the system, respectively (Fig. 2a-b). These modes correspond to the extremes in the transmission graph displayed in Fig. 2e, and their energy density projection as they propagate through the system is displayed in Fig. 2d.
While the scattering matrix shown in Fig. 2f is sufficient to control transmission. The controlled deposition of energy at the centre of the sample, as shown in Fig. 2c, requires full knowledge of the complete deposition matrix 23 of size 640, 000 × 1022. Focussing onto the infinity-shaped path inside the scattering system (Fig. 2c) was done by a pseudo-inversion of the matrix (see Methods).

Discussion
The timing comparisons of Fig. 1d and e show that electromagnetic field calculations are more efficient by two orders of magnitude and, conversely, that problems can be addressed with a Fig. 3 Evaluation of the influence of sampling on accuracy and the simulation size. The finite-difference frequency domain method is compared with the Fourier-space convolution used in the proposed method. (a) The root-mean square error as a function of sampling density for a problem size of 16λ × 16λ. The convolutional neural network is more accurate at 4 samples per wavelength than the finite difference method is at 64 samples per wavelength, with 2.6% versus 3.6% error respectively. (b) The maximum calculation size as a function of calculation time for MaxwellFDFD 24,25 , and the proposed algorithm on the same CPU. Sampling densities of 64 and 4 samples per wavelengths were used to ensure that the root-mean-square error is below 3.6% (see Fig. 3a).
size that is a 100-fold larger. The efficiency advantage is already visible for 2D calculations with 10 4 sample points, and it grows approximately linearly with the problem size.
Compared to other algorithms, such as FDFD, our neural network method is significantly more accurate for larger step sizes, it can compute the same problems much faster and it requires less memory to do so, enabling the calculation of larger systems, as shown in Fig. 3. To compute an electric field matrix, the method requires storage of 7 additional matrices of the same size, which can be reduced to 4 by a sacrifice of some computation speed.
In conclusion, we demonstrated how Maxwell's equations can be rephrased as a recurrent neural network that yields their solution without the need for training. This stands in sharp contrast to conventional machine learning approaches, where the connection weights are determined during a computationally expensive learning phase that is prone to introduce training bias. Our recurrent neural network has a single hidden layer which performs a shifted Green's function, while the other connection weights have a one-to-one correspondence to the problem's permittivity distribution. By incorporating preconditioning from the outset 26 , convergence is efficient and fully deterministic. This auspicious recurrent neural network structure enabled us to leverage powerful, yet publicly available, machine learning infrastructure on the cloud. We integrated this into the open source electromagnetic calculation library, MacroMax. This enables any researcher with an internet connection to solve 100-fold larger complex scattering problems.

Methods
Implementation of the solver. The electromagnetic solver was implemented using the PyTorch (1.10.1+cu111) machine learning library 27 , and integrated into the MacroMax library 9 . The widely used PyTorch framework was chosen for its built-in ability to handle complex numbers and perform fast Fourier transforms. Its ubiquity allows us to leverage the latest technological advances and cloud computing platforms. The training of the forward neural networks, depicted in Fig. 1a-b, was done using PyTorch's Adam optimiser. We found that Stochastic Gradient Descent converged slower for this network topology. A fast-Fouriertransform-based convolution layer was used to avoid finite differences and handle the Laplacian operation more efficiently than the built-in version. The recurrent neural network of Fig. 1c also implements the Green's function as such convolution. Unlike the preceding topologies, this network has fixed weights and the solution as output. As such, it does not require training, and can infer the electric field on mm-scale. We integrated this approach with the electromagnetic solver, MacroMax 13 , to make it more accessible and extend its use to e.g. birefringent and magnetic materials.
Scattering system parameters. Fig. 1f-g was calculating for a wavelength of 633 nm. The scattering layer consists of packed spheres with a radius of 30 µm ± 10 % and a refractive index of 1.33. It is 1.25 mm thick. The system is padded with 16-wavelength thick absorbing boundaries with a linearly increasing extinction coefficient from 0 to 0.2.
In Fig. 2a-c we use a 500 nm light source. Similarly, the system contains a scattering system made of spheres with radius of 2.0 µm ± 5 % and a refractive index of 1.33. In Fig. 2a-b these spheres compose a scatterer that is 64 µm thick. The layer in Fig. 2c is thicker (85 µm) but has a rectangular gap in the centre of the target region. For the purpose of the mode transmission diagram in Fig. 2e, the scattering system was padded with 25λ thickness boundaries, which have a linearly increasing absorption coefficient from 0 to 0.025. This was done to reduce reflection effects from the edges of the system, which interfered with the accurate calculation of displayed values. Otherwise, 16λ boundaries were used, with a linearly increasing absorption coefficient from 0 to 0.2.
Calculation of deposition and scattering matrices. We compute the internally and externally scattered fields for a complete basis of incident plane waves. To facilitate the study of transmission matrices, we assume a wide slab geometry for the scatterer, orthogonal to the z-axis. To minimise edge effects for a finite slab, we adopt periodic boundary conditions in the transverse dimension. Common fast Fourier components are used for inward and outward travelling waves so that free space propagation corresponds to the identity matrix. The code included with this manuscript lists the basis vectors in raster-scan order for the forward and backwards propagating waves, respectively. Although deemed out of scope, pairs of orthogonal propagating polarisations are listed for vector fields. In total, the scattering and deposition matrices used to produce Fig. 2 have 1022 columns.
The deposition matrix used for Fig. 2c has 640, 000 rows, one for each internal and external value of the field. While the scattering matrix is a square matrix with the same plane wave basis for the row space. The scattering matrix can be considered a 2 × 2-block matrix with four quadrants. The two quadrants on the diagonal correspond to the forward and backward transmission matrices, while the off-diagonal quadrants are the front and back-reflection matrices 14 .
To get the open and closed channels ( Fig. 2a-b), we calculated the singular value decomposition of the forward transmission matrix, (i.e. top-left quadrant of scattering matrix in Fig. 2f) using LAPACK's algorithm based on a divide-andconquer method.
A pseudo-inversion of the deposition matrix was done to refocus within the deposition region inside the scattering material. The singular value decomposition of the deposition matrix was calculated using ARPACK's implicitly restarted Lanczos algorithm. The large computer memory requirements of this matrix restricted us to 750 eigenvectors for use in the inversion, instead of full 1022.
Evaluation of time efficiency. In Fig. 1d-e, the median iteration time is used to avoid giving excessive weight to the initial iteration which includes the initialisation steps. While the initialisation overhead is negligible for typical problems, it hampered a direct comparison for the smallest problem sizes. Typical systems require around 1000 iterations to converge, depending on the range of permittivity and the size of the system. Vectorial calculations with a polarised electric field take approximately 3.9 times longer than scalar calculations with the same number of sample points.

Data availability
All data underlying the results was generated by the algorithm described in this manuscript. The complex-valued scattering matrix shown in Fig. 2f will be made available for download after the manuscript is accepted.

Code availability
The algorithm, as well as the data visualisation, is implemented in Python using PyTorch. The complete source code with examples is openly available as a Git repository 13 . The PyTorch implementation is integrated in the MacroMax electromagnetic calculation library, which is freely available on the Python Package Index 28 .