Introduction

Calculators of various kinds have emerged by forging numerical algorithms with the corresponding technological platforms. While the algorithms describe the mathematical paths on how solutions to problems can be found, the platforms are responsible for the transliteration of this abstract path into measurable quantities. The algorithms, the platforms, and their fusion define such systems’ features and limitations. Following the ever-growing demand for ultrafast, compact, low/near-zero-power, and integrable cyber-physical devices for mathematical computations, it is organic that significant research efforts focus on making these numerical systems as optimal and efficient as possible.

This quest led to the exploration and development of unconventional analog computing systems that exploit electromagnetic waves to deliver low-power, ultrafast, compact, parallelized computations1,2,3,4,5,6,7. The two main categories in this domain involve systems that utilize free-space (scattering) elements (e.g., lenses6), and waveguides (e.g. photonic systems8,9 and phased arrays10). Sufficient free-space propagation can act as dense matrix multiplication11. Realized with traditional optics, this results in bulky devices6,12, while metasurfaces can be more compact13,14,15. However, in both, there can be major bottlenecks regarding photonic and electronic integration. Waveguiding systems such as photonic integrated circuits (IC) offer more mature solutions for integrable and reconfigurable devices at the expense of the challenge of rendering all of the needed functionality in a 2D plane16. In all cases, their main challenge is reconfigurability, since its implementation requires some form of a-priori mathematical calculations. For instance, metasurfaces/complex media17,18 require optimization, and photonic meshes require operator decomposition19,20,21,22.

In terms of their mathematical abilities, the above examples demonstrate wave-based analog computing with functionalities such as integration/differentiation in space23,24,25,26 and time27, matrix-vector multiplication28, emulating equations through physical phenomena29,30, or acting as platforms for neural network functionalities3,6,9,31. The intersection with the metamaterial paradigm delivered a series of remarkable analog computing devices with matrix multiplication1,7,23 and ultimately equation-solving (matrix inversion) capabilities2. In most cases, the matrix computations (especially the matrix inversion2) were performed through stationary algorithms32, such as the Jacobi method, where the matrix (operator/kernel) does not change with the iteration count.

The fundamental and far-reaching question we address here is whether a wave-based analog metastructure can be reconfigurable simply and intuitively, without needing a-priori calculations. Most importantly, the resolution to this question endows one with the ability to implement stationary and non-stationary algorithms. We propose a device based on an RF waveguide architecture with reconfigurable components. Regarding stationary problems, we use this device to perform a matrix inversion of a random (statistically uncorrelated) large set of matrices. As for non-stationary problems, we demonstrate both root finding using Newton’s method and Inverse Design. All three examples are possible only due to the reconfigurability of the device and hint at the possibility of deeper explorations into the realm of advanced numerical algebra methods.

Results

A conceptual representation of the main idea is pictorially summarized in Fig. 1A. The main property of our proof-of-concept system distinctively different from all previous metastructure approaches (e.g.,2) is its reconfigurability; the metastructure has the ability to rapidly take on different matrices (operators or kernels) K. To facilitate this, we employ a wave-based direct complex matrix (DCM) architecture, which offers an intuitive and simple implementation of any desired matrix33. Using waves instead of currents and voltages, it is analogous to the crossbar architecture used in electronic analog computing systems34,35,36,37, and it can be seen as a generalized phased array feed10. In this device, a collection of n × n tunable phase shifting and amplifying elements (which can also act as attenuating elements) connect an input vector of n complex amplitudes on an array of transmission lines to a similar output vector through combiners. This architecture can be seen schematically in Fig. 1B and its corresponding experimental implementation in Fig. 1C.

Fig. 1: A reconfigurable wave-based analog computing metastructure.
figure 1

A A conceptual representation that describes the main objective, i.e., a reconfigurable device that can provide us with repeated matrix inversions of arbitrary matrices in order to achieve stationary and non-stationary algorithms. The device can implement any given kernel K = I − A and give the \({\left(I-K\right)}^{-1}={A}^{-1}\). B The central component of the design consists of a direct complex matrix (DCM)33 architecture of 5 × 5 elements. C The experimental realization of this design for the 45 MHz operating frequency. D The essential element of the DCM is the multiplier module, consisting of both a phase shifting and an amplification part (which can also function as an attenuation part); controlled with an onboard microcontroller unit. Finally, E and F depict the performance of the DCM machine in the open-loop (matrix-vector multiplication) and the closed-loop (matrix inversion) setups. In both, we compare the experimentally obtained matrix to one computed conventionally. We see good agreement.

Full size image

The key component of the metadevice is the multiplier module (Fig. 1D) so named because given an input signal characterized by its complex amplitude at 45 MHz, Vin, it will render a similar output Vout = zVin, where z is a complex multiplication factor. This module consists of two basic components: (a) a voltage-controlled phase shifter with over 360 degrees of potential phase rotation, and (b) a voltage-controlled amplifier with  â‰ˆ47 dB of dynamic range (−30 dB to +17 dB). Each multiplier contains an embedded microcontroller unit (MCU), which allows each device to be controlled externally through a suitable communication network and a computer (see supplementary material). While the design frequency is 45 MHz, the module could be implemented at RF (GHz) and photonic (THz) platforms, following the same principle of operation. The experimental DCM implementation consists of 25 multipliers to yield 5 × 5 complex matrices. The ingress and egress stages33 are implemented with five 1-to-5 power splitters (ingress stage) and five 5-to-1 signal combiners (egress stage). The multipliers are clustered into five groups, one for each matrix row. In Fig. 1B we depict the planar schematic of the DCM suitable for photonic implementation. However, for the RF implementation, we stacked and routed the components vertically (Fig. 1C), making the device compact for our particular wavelength and platform choice. Different stacking or integrated circuitry approaches can potentially be used to further reduce its overall footprint.

The metadevice can be operated in one of two configurations with dramatically different results. When the DCM is set in an open-loop configuration (Fig. 1E inset), it can be used to rapidly calculate parallelized matrix-vector multiplication. However, a closed-loop configuration can be created by connecting the outputs and the inputs with a feedback loop using properly designed couplers (Fig. 1F inset). When the DCM is in a closed-loop configuration, the metadevice can rapidly calculate parallelized matrix inversion (equation solving). This is a unique feature of metadevices/metastructures2,33 that incorporate feedback loops.

Inverting stationary matrices

First, we investigate the stationary analysis capabilities of our metadevice. For the assessment of the open- and closed-loop operation, we performed a series of randomized trials, one instance of which is presented in Fig. 1E, F. For each trial, a random passive matrix \(A\in {{\mathbb{C}}}^{5\times 5}\) was chosen and applied to the metadevice in both configurations. Each measurement was performed by exciting each input port in turn with all other inputs appropriately terminated and then observing the complex amplitudes on each of the output ports. For the open-loop configuration, this corresponds to performing five matrix-vector multiplications, or as A ⋅ I where each column of I is progressively applied (one column at a time) as separate vectors. The closed-loop configuration was measured similarly, but in this case, we are probing the steady state of the metadevice, which corresponds to (I − K)−1 ⋅ I = A−1 ⋅ I. While this measurement technique fully characterizes both configurations, in practice dense complex vectors will be input and read to achieve parallelized results.

The estimated relative error ∥Aexact − Ameas∥2/∥Aexact∥2 for both cases (Fig. 1E, F) revealed an error of about 0.001 and 0.005, respectively. Despite the component imperfections, misalignments, measurement noise, and other stochastic errors, the measured results are in excellent agreement with the theoretical values. The calibration procedure of the metadevice and the statistical analysis of the full trial set (100 values) are presented in the Supplementary Material.

A single multiplier module has a rise time of approximately 80 ns to achieve its desired complex value. This value is approximately 4T assuming one-period duration of T = 1/45 MHz ≈ 22.2 ns. In the open-loop configuration, the total response requires approximately 5T, including signal delays in connections and splitters. The duration of the closed-loop case is affected by the platform and the condition number of the inverted matrix33, but in principle is in the same order of magnitude. Possible photonic implementations may further reduce this time to the picosecond range3 and below38.

Root finding with Newton’s method

We now apply the implemented metadevice to two characteristic non-stationary problems that highlight its mathematical abilities: (i) root finding of a system of five equations with five unknowns using Newton’s iterative technique and (ii) implementing an inverse design problem using the Lagrangian multiplier formalism for constrained optimization. Both cases require that the kernel be reprogrammed in each iteration step. Note that our approach is not restricted to these two problems; instead, we choose these to highlight the potential of the introduced metadevice.

For the first case, we construct a simple nonlinear toy problem, and we apply Newton’s algorithm39 (Fig. 2A) to find one possible root. The vector problem statement reads

$${{\bf{f}}}({{\bf{z}}})={\left[{f}_{1}({{\bf{z}}}),{f}_{2}({{\bf{z}}}),{f}_{3}({{\bf{z}}}),{f}_{4}({{\bf{z}}}),{f}_{5}({{\bf{z}}})\right]}^{T}={{\bf{0}}}$$
(1)

where \({{\bf{f}}}\in {{\mathbb{C}}}^{5\times 1}\), \({{\bf{z}}}={\left[{z}_{1},{z}_{2},{z}_{3},{z}_{4},{z}_{5}\right]}^{T}\in {{\mathbb{C}}}^{5\times 1}\), and 0 is the zero vector. We construct the vector function to have the following polynomial form

$${f}_{1}({{\bf{z}}})=({z}_{1}-{r}_{1})({z}_{2}-4.2i)({z}_{3}+2)({z}_{4}-5i)({z}_{5}-3.5)$$
(2)
$${f}_{2}({{\bf{z}}})=({z}_{1}-3.9)({z}_{2}-{r}_{2})({z}_{3}+2.5i)({z}_{4}-3.2i)({z}_{5}-4.2)$$
(3)
$${f}_{3}({{\bf{z}}})=({z}_{1}+5.2i)({z}_{2}-4)({z}_{3}-{r}_{3})({z}_{4}-4i)({z}_{5}-7.1)$$
(4)
$${f}_{4}({{\bf{z}}})=({z}_{1}-3)({z}_{2}-7i)({z}_{3}+4)({z}_{4}-{r}_{4})({z}_{5}-5i)$$
(5)
$${f}_{5}({{\bf{z}}})=({z}_{1}-5.2i)({z}_{2}-4)({z}_{3}+4.75i)({z}_{4}-8)({z}_{5}-{r}_{5})$$
(6)

where \({{\bf{r}}}={\left[{r}_{1},{r}_{2},{r}_{3},{r}_{4},{r}_{5}\right]}^{T}\) are the vertices of a regular pentagon with 1/4 radius (see Fig. 2B) and the other factors represent additional extraneous roots far from the starting point.

Fig. 2: Experimental verification of Newton’s root finding method with the proposed metadevice.
figure 2

A A comparison between the experimental results and a numerical implementation of the algorithm. The faint solid lines are cases where additional stochastic noise has been added to the system. We observe that the experimental trajectory is well contained within these simulated noisy paths. B A comparison between the numerical (left column) and measured (right column) evaluation of the inverse Jacobian at different iterations.

Full size image

For the evaluation of Newton’s method, we need to calculate the Jacobian matrix, i.e., \({J}_{ij}=\frac{\partial {f}_{i}}{\partial {z}_{j}}\) or

$${J}_{f}({{\bf{z}}})=\left(\begin{array}{cccc}\frac{\partial {f}_{1}}{\partial {z}_{1}}&\frac{\partial {f}_{1}}{\partial {z}_{2}}&\cdots \,&\frac{\partial {f}_{1}}{\partial {z}_{5}}\\ \frac{\partial {f}_{2}}{\partial {z}_{1}}&\frac{\partial {f}_{2}}{\partial {z}_{2}}&\cdots \,&\frac{\partial {f}_{2}}{\partial {z}_{5}}\\ \vdots &\vdots &\ddots &\vdots \\ \frac{\partial {f}_{5}}{\partial {z}_{1}}&\frac{\partial {f}_{5}}{\partial {z}_{2}}&\cdots \,&\frac{\partial {f}_{5}}{\partial {z}_{5}}\end{array}\right)$$
(7)

therefore the root can be estimated by the following iterative process

$${{{\bf{z}}}}_{n+1}={{{\bf{z}}}}_{n}-\alpha {J}_{f}^{-1}({{{\bf{z}}}}_{n}){{\bf{f}}}({{{\bf{z}}}}_{n})$$
(8)

where α = 0.2 is a relaxation constant33.

Table 1 shows the required algorithm steps to implement the iterative scheme described by Eq. (8). Note that the Jacobian changes value in each iteration and that its inverse is required to be calculated anew. The results are then used to update the z. The method converges successfully after a few iterations.

Table 1 Root finding with Newton’s algorithm and DCM metadevice
Full size table

Matrix inversion represents a computationally intensive task, expedited by our metadevice. In Newton’s method, matrix inversion stands as the primary complexity, with an \({{\mathcal{O}}}({n}^{3})\) complexity assumption for an n × n matrix40. On the other hand, evaluating the Jacobian typically demands \({{\mathcal{O}}}(n)\) operations, facilitated by automatic differentiation techniques29. The efficient computation of the Jacobian at each step can be provided to the system either from a dedicated subsystem responsible for its calculation and feed into the main system.

A numerical version (using MATLAB) is compared with the experimental results illustrated in Fig. 2B. We observe that for both MATLAB and the experiment, the estimation vector converges close to the exact roots. Furthermore, the estimated vector reaches a stationary point as the iteration count increases. After 15 iterations, the relative error is ∥z − r∥2/∥r∥2 ≈ 0.0023. This is similar to the accuracy achieved for stationary trials, thus representing the accuracy floor of our system. A similar picture is also visible by comparing three specific iterations, as illustrated in Fig. 2C, where a comparison of the full Jacobian is presented.

The experimental results do not precisely follow the paths indicated by the numerical implementation realized using MATLAB. This can be explained by adding random noise to the Jacobian on each iteration step. The added noise \(N\in {{\mathbb{C}}}^{5\times 5}\) is a random complex matrix that follows a normal distribution inside a disk with radius \({r}_{N}=0.01{\lambda }_{\max }\), where \({\lambda }_{\max }\) is the maximum eigenvalue of the Jacobian. The noise creates many possible paths, all of which successfully converge, and we observe that our measured results comfortably lie within these families of curves. Note that some solution branches are more susceptible to this noise than others (e.g., r1 (blue) and r3 (red) curves in Fig. 2B) and this is due to the details of the toy problem solved.

Generally, the device’s numerical accuracy has a threshold that depends both on the implementation and on the measuring apparatus (vector network analyzer (VNA)). In our case, the closed-loop operation was on the order of a few microvolts, approaching the sensitivity floor of the VNA utilized (refer to Supplementary Material for details). Additionally, the implementation would influence the scalability of systems handling matrices with tens or hundreds of inputs. Thoughtful RF/Photonic IC designs might pave the way for incorporating larger kernels, albeit requiring more sensitive VNAs. It is crucial to emphasize phase and amplitude stability, as highlighted in ref. 33.

When higher-precision computations are required, this device can be a part of a mixed-precision computing system. In these systems, part of the calculations are done in a fast, low-precision estimation stage and then fed and further refined at a higher-precision stage, similar to the in-memory mixed-precision approaches in electronic platforms41.

Inverse design of a metastructure

For the second example, we chose the case of an inverse design problem (Fig. 3A). We assume that our design consists of a collection of m = 5 two-dimensional (2D) scatterers with circular cross-section at fixed known locations r = [r1, â€¦, r5], each with an unknown bounded permittivity \({{\mathbf{\varepsilon }}}=[{\varepsilon }_{1},\ldots,{\varepsilon }_{5}]\in {{\mathbb{C}}}^{5\times 1}\). The goal is to achieve a specific user-defined scattered field measured at a series of n = 4 detection (objective) points, o = [o1, â€¦, o4]. Note that in our case we assume a collection of cylindrical circular scatterers (2D) excited with a monochromatic incident field of λw wavelength. The x-propagating incident field (kx) is polarized in the z-direction (TE wave - Ez) with the ejωt convention.

Fig. 3: A metastructure that designs a metastructure: numerical and experimental results.
figure 3

A Schematic of the numerical test case. A set of five two-dimensional (2D) scatterers with unknown permittivities, to be determined via our analog metadevice, \(\varepsilon=\left[{\varepsilon }_{1},\ldots,{\varepsilon }_{5}\right]\) subject to a plane wave excitation. Pictured is the complex-valued scattered field. The scattered fields at the observation points [o1, â€¦, o4] are used as the benchmark values of this problem in which the objective fields are shown as the color in each torus. The algorithm tunes each permittivity in order to match the scattered field (center of each torus) to the objective fields. B The relative error \({{\mathcal{E}}}\) for the algorithm computed both numerically, and experimentally using the metadevice, under various noise and filtering scenarios. The experimental device gives a minimum relative error of 0.00172 at 87 iterations. C Objective permittivities (black rings) compared to those computed numerically (blue cross) and experimentally via the metadevice (yellow rhombus) using the described algorithm at iteration 87. D Evolution of scattered field vector up to iteration 87 in comparison to objective fields for experiment and simulation under various noise and filtering scenarios.

Full size image

The scatterers are coupled, making this a nonlinear problem modeled using the Lippmann–Schwinger42 scheme, solved with a standard discrete dipole approximation (DDA) methodology43. Each scatterer will respond to the local (self-excluded) electric field, which consists of the known incident electric field, \({e}_{{{\rm{inc}}}}\in \,{{\mathbb{C}}}^{5\times 1}\), and the scattered field from all other scatterers, \({e}_{{{\rm{sca}}}}\in \,{{\mathbb{C}}}^{5\times 1}\). The scatterers exhibit a complex polarization vector \(p=A({e}_{{{\rm{inc}}}}+{e}_{{{\rm{sca}}}})\in \,{{\mathbb{C}}}^{5\times 1}\) where A is the normalized polarizability diagonal matrix, i.e., A = diag(ε − εbackground). The field interaction between the scatterers is expressed via the Green matrix \(G\in \,{{\mathbb{C}}}^{5\times 5}\) (hollow symmetric matrix) such that esca = Gp. We may express the polarization vector p = A(einc + Gp), which indicates the mutual dependence of p. Therefore, the polarization vector can be calculated as \(p={({A}^{-1}-G)}^{-1}{e}_{{{\rm{inc}}}}\). Finally, we use the four objective points o to measure the scattered field vector \({e}_{{{\rm{meas}}}}={G}_{{{\rm{pr}}}}p\in \,{{\mathbb{C}}}^{4\times 1}\) where \({G}_{{{\rm{pr}}}}\in {{\mathbb{C}}}^{4\times 5}\) is the propagator Green function. The measured field is then compared to a (user-defined) objective \({e}_{{{\rm{obj}}}}\in \,{{\mathbb{C}}}^{4\times 1}\).

A typical constrained minimization problem (primal) can be written as39

$$\begin{array}{rcl}&&{\min }_{\begin{array}{c}x,y\end{array}} \; \; f(x,y)\\ &&{{\rm{s.t.}}} \quad g(x,y)\le 0\end{array}$$
(9)

where f(x, y) are the objectives and g(x, y) are the constraints. For such problems, the Lagrangian (dual) problem is expressed as

$${\max }_{\lambda }{\min }_{\begin{array}{c}x,y\end{array}}\quad {{\mathcal{L}}}(x,y,\lambda )=f(x,y)+\lambda g(x,y)$$
(10)

Note that x and y may be subject to further requirements such as domains and bounds.

For our particular example, we have that x = p, y = ε, and f(p, Îµ) = 1/2∥Gprp − eobj∥2 and g(p, Îµ) = 1/2∥(A(ε)−1 − G)p − einc∥2. In our formulation, the objective is the scattered field at the observation points. The constraints comprise the self-consistency of the polarization vector (physics). Also, the permittivity vector is subject to specific bounds, i.e., \(\varepsilon \in {\mathbb{R}}\) and ε ∈ [1, 5]. Note that g(p, Îµ) is nonlinear with respect to p and ε and therefore requires a non-stationary approach.

Following an initialization, our numerical evaluation of the above is implemented by a non-stationary algorithm that requires repeated application of the following three steps. First, we minimize with respect to ε by examining \({\nabla }_{\varepsilon }{{\mathcal{L}}}(p,\varepsilon,\lambda )=0\). At this step, we project the resulting permittivity vector to the desired domain and bounds. Second, we minimize with respect to p by examining \({\nabla }_{p}{{\mathcal{L}}}(p,\varepsilon,\lambda )=0\). At this stage, the required stationary matrix inversion is performed with our metadevice. Finally, we maximize for λ by using \({\nabla }_{\lambda }{{\mathcal{L}}}(p,\varepsilon,\lambda )=0\). These steps are repeated until convergence is achieved, i.e.,

$${{\mathcal{E}}}=\parallel \!\!{e}_{{{\rm{obj}}}}-{e}_{{{\rm{meas}}}}{\parallel }^{2}/\parallel \!\!{e}_{{{\rm{obj}}}}{\parallel }^{2} < \delta$$
(11)

(for more information, see SM).

As a numerical test case, the scatterers are assumed to be lossless with permittivity of \(\varepsilon=\left[{\varepsilon }_{1},{\varepsilon }_{2},{\varepsilon }_{3},{\varepsilon }_{4},{\varepsilon }_{5}\right]=[3.5,\,1.5,\,1.5,\,3.5,\,1.5]\). The objective scattered field at the detection points o, as depicted in (Fig. 3A), is eobj = [−0.0086 − 0.0078j, 0.0089 − 0.0132j,−0.0066 − 0.0120j, 0.0043 − 0.0004j]. The values were extracted from the DDA method and verified with a full-wave COMSOL simulation. Note that Fig. 3A depicts the complex (hue/saturation) of the electric scattered field (Ez), i.e., the difference between the total field and the incident excitation.

Figure 3B depicts a set of four cases for the same algorithm. In the first case (black line), the idealized (noiseless, no filtering) computer evaluation of the algorithm is given - we observe that after only 20 iterations the error drops below 10−3. The experimental results are presented in Fig. 3B as red dots. The measured results exhibit an optimal point (minimum error) after 87 iterations, with an error of 0.00172. As an analog device, there is an additional systematic/stochastic/experimental noise to the system which affects the fidelity of the matrix inversion. We apply a simple averaging filtering scheme on the polarization estimation, i.e., pnew = (1 − αF)p + αFpprevious, with αF = 0.25, as a way to partially mitigate this noise. The filter affects the convergence speed by increasing the iteration count but also significantly improves the accuracy/fidelity of the matrix inversion, hence the metadevice’s performance. This feature is illustrated in Fig. 3B, where the retrieved experimental results are compared to the idealized computer evaluation with the applied filter (blue line). We also performed a series of 100 randomized cases of the idealized filtered computer evaluation with added noise to the estimated/measured polarization vector (faint blue lines in Fig. 3B). The noise profile is similar to the one used in the first example (Newton’s method). The measured results are well contained within these error bounds. Note that iteration count is not equivalent to time. For a traditional computer evaluation, each iteration (with its required matrix inversion) could ultimately be slower than the convergence time of an optimized hardware implementation of the metadevice.

Due to systematic/measurement noise, the error begins to grow after the experimental accuracy floor is obtained - an indication that a termination criterion could be applied at this point. This result also agrees with the maximum accuracy we obtained in the previous non-stationary example. More sophisticated error-correcting and filtering schemes can push the accuracy below this threshold. For instance, αF could be adaptively tuned during the non-stationary evaluation to realize a mixed-precision computing system.

At the minimum error point (iteration 87), the extracted permittivity estimation is illustrated in (Fig. 3C). Notice that the values are very close to the numerical test case objectives and permittivities. Finally, Fig. 3D illustrates the path of the scattering vector, emeas, for these 87 iterations. Similar to the above example, the faint paths represent the added noise effects on the numerical evaluation.

Discussion

For both presented non-stationary examples, it is evident that our metadevice can act either as an ultrafast analog computing machine and mathematics calculator with waves, or in a broader sense as an electromagnetic emulator for inverse design44. It can be used for a plethora of realistic problems where the linear response of a system (i.e. matrix-vector multiplication) or the solution of a system of equations (stationary problems, matrix inversion) is required. Active control systems45, real-time signal processing46, and radar systems47 are a few indicative systems that require real-time analog processing. Moreover, the intuitive reconfigurability of this metadevice also enables the performance of constrained optimization tasks, like the ones required in non-stationary problems such as inverse design, where the desired response of complex media requires intensive optimization48. In short, this metastructure can design metastructures. Finally, an adaptation of the above proof-of-concept metadevice in RF-IC, photonic, or hybrid platforms can make it an excellent candidate for on-the-fly or computation-through-propagation ultrafast, parallelized calculations. Scaling such metastructures to larger matrices is an open challenge.