Reprogrammable, in-materia matrix-vector multiplication with floppy modes

Abstract

Matrix-vector multiplications are a fundamental building block of artificial intelligence; this essential role has motivated their implementation in a variety of physical substrates, from memristor crossbar arrays to photonic integrated circuits. Yet their realization in soft-matter intelligent systems remains elusive. Here, we experimentally demonstrate a reprogrammable elastic metamaterial that computes matrix-vector multiplications using floppy modes—deformations with near-zero stored elastic energy. Floppy modes allow us to program complex deformations without being hindered by the natural stiffness of the material; but their practical application is challenging, as their existence depends on global topological properties of the system. To overcome this challenge, we introduce a continuously parameterized unit cell design with well-defined compatibility characteristics. This unit cell is then combined to form arbitrary matrix-vector multiplications that can even be reprogrammed after fabrication. Our results demonstrate that floppy modes can act as key enablers for embodied intelligence, smart MEMS devices and in-sensor edge computing.

I Introduction

Linear transformations, encoded in matrix-vector products, are ubiquitous in modern computing; with applications spanning from deep learning to 3D graphics [1] and data compression [2]. In fact, universality theorems prove that any function can be approximated by a network of programmable linear transformations, interconnected by fixed nonlinear activation functions [3]. This widespread applicability has motivated the development of specialized hardware to compute matrix-vector products, ranging from the well-established Graphics Processing Unit (GPU), to emergent physical computing approaches such as photonic integrated circuits [4] and memristor crossbar arrays [5].

Within the paradigm of physical computing, mechanics has gathered considerable excitement for a range of special-purpose but highly relevant applications [6]: Mechanics offers the capability to directly process signals such as speech [7] or gait [8] without transduction into the electrical domain—facilitating zero-power, in-sensor information processing; it forms the basis for embodied intelligence in soft robotic structures—enabling such robots to mimic the agility and dexterity of living organisms [9, 10, 11], and, when scaled down to microscopic dimensions [12, 13], can perform computations with ultra-low energy consumption [14, 15]—enabling intelligent devices that can operate under extreme power constraints. The overarching challenge in achieving mechanical computation lies in mapping relevant information-processing operations, such as matrix-vector multiplications, into the natural responses of mechanical systems.

Figure 1: A metamaterial (red) computes a matrix-vector multiplication by decomposing it into elementary operations (blue), corresponding to individual matrix elements. Each elementary operation (purple) is implemented by a metamaterial unit cell. The input vector is presented as a displacement applied on the input control rods (orange, $x_{i}$) and the result is deterined by measuring the displacement of the output rods (blue, $y_{i}$). Each internal output is connected to an input, ensuring that the resulting metamaterial is frustration-free irrespective of the particular matrix being implemented.

In this work, we experimentally demonstrate a soft mechanical metamaterial that computes a matrix-vector multiplication using floppy modes (zero modes). The paper is structured as follows: We start by decomposing the matrix-vector product into a network of elementary operations, each of which will be computed by a single metamaterial unit cell (Fig. 1). Then, we implement this unit cell in a planar network of point masses subject to ideal constraints. The resulting unit cell presents two zero modes that branch and cross—as required by the matrix-vector multiplication, and preserves deformation compatibility when tiled. Afterwards, we investigate, using Finite Element Method (FEM) simulations, how a real metamaterial geometry deviates from the ideal point-mass model model. We use automatic differentiation [16] to adjust the design and preserve accuracy in the presence of nonzero bending stiffness and finite compliance, and discuss the limitations of the metamaterial arising from such real-world effects. Finally, we experimentally realize the metamaterial using water-jet cut rubber, and characterize its displacement using optical flow. We experimentally characterize the metamaterial unit cell, a 2x2 matrix-vector product, and a tunable unit cell that can be reversibly switched between multiple matrix coefficients. We observe that the metamaterial can compute the desired matrix-vector product with remarkable accuracy, even in the presence of fabrication tolerances and material non-idealities. These results demonstrate that floppy modes can be harnessed to implement intelligent responses in elastic systems.

II Design

Among the deformation responses of mechanical systems, floppy modes, or zero-modes, are unique in the fact that they are associated to a vanishingly low amount of stored elastic energy, i.e., they require small forces to actuate. Thus, they are a prime candidate for computation in passive elastic systems where deformations must be driven by the enegy available in the signal. Floppy modes originate in a rank deficiency of the rigidity matrix, caused by either an insufficient number of constraints [17, 18, 19, 20], or by the topological properties of the system [21]. Their global nature and topological connections point towards exciting research questions, but also pose a challenging design problem for applications that require a large number of on-demand floppy modes, such as matrix-vector multiplication (a matrix-vector product requires a number of modes equal to the dimension of the input space). Several approaches towards on-demand floppy modes have been introduced: Combinatorial methods draw from a small set of unit cells with well-defined compatibility characteristics [22, 23, 17]. They can be highly efficient, but are limited by the small number of unit cells and large design space. Techniques based on iterative addition of constraints can accurately produce designs with a few desirable modes, but introduce additional unwanted modes as a by-product [24]. Finally, machine learning methods can identify the number of modes corresponding to a particular unit cell [25], but they do not yet possess the capability to solve inverse design problems. While progress has been made in the particular case of line modes [26, 27, 28, 29, 23], branches and crossings—necessary for the realization of mechanical information-processing circuits—still pose a challenge for these methods. We will tackle the challenge of designing a material with a large set of complex floppy modes by introducing a unit cell with well-defined compatibility characteristics, that can be combined to achieve the desired response.

In this section, we introduce a metamaterial design that performs a matrix-vector multiplication $\vec{y}=A\vec{x}$. The input vector $\vec{x}$ is presented to the system by prescribing the displacement of a set of input degrees of freedom, and the output result $\vec{y}$ is read out by measuring the displacement of a different set of degrees of freedom. Since the system should be able to respond to any (small) input vector with no resistance, the space of floppy modes should span all possible input displacement combinations. At the same time, we require the output displacement to be entirely prescribed by the input vector, meaning that the system should not present additional floppy modes. For large input spaces, this involves designing a large number of modes—beyond the limit of direct optimization methods. We can overcome this limitation by decomposing the matrix-vector multiplication as a planar network of computation tiles (Fig 1), whose outputs $u_{1}$ and $u_{2}$ are related to the inputs $v_{1}$, $v_{2}$ as $v_{1}=A_{ij}u_{1}+u_{2}$ and $v_{2}=u_{1}$. Each tile can be understood as individually implementing a small matrix-vector multiplication, with the form

$$\begin{pmatrix}v_{1}\\ v_{2}\end{pmatrix}=\begin{pmatrix}A_{ij}&1\\ 1&0\end{pmatrix}\begin{pmatrix}u_{1}\\ u_{2}\end{pmatrix}.$$ (1)

Each individual computational tile is physically implemented by a metamaterial unit cell. Since unit cells are only required to perform a 2$x$2 matrix-vector multiplication, they require only two floppy modes. This simplified requirement can be met with existing floppy mode design methods—here, we solve it by hand-designing the unit cell.

Figure 2: (a) single degree of freedom, the angle of the support beams control the coupling between x and y (b) the unit made of 9 DOFs. The angle at site 3 and 6 (0 and 1) control the coupling between $u_{1}$ ($u_{2}$) and $v_{1}$. (c) the 2 floppy modes of the system

By construction, the tiling in Fig. 1 exhibits the correct number of global floppy modes: Each individual row added to the system will add $2n_{col}$ floppy modes, but will also incorporate $2n_{col}-1$ constraints ($n_{col}$ constraining the input $u_{2}$ with the output $v_{1}$ of the upper row, and $n_{col}-1$ intra-row constraints). This will increase the number of floppy modes by one, reflecting the fact that the extended system contains an extra input. Adding an additional column will add $2n_{row}$ floppy modes, but also $2n_{row}$ constraints ($n_{row}$ constraining the input $u_{1}$ of each site with the output $v_{2}$ of the previous column, $n_{row}-1$ intra-column constraints plus an extra constraint setting the input $u_{2}$ of the topmost site to zero, as shown in Figure 1). Since the constraints always connect an output (by definition prescribed) with an input (by definition free), and the directed computation graph presents no cycles, the constraints will neither neither be redundant nor introduce frustration. Such designer compatibility relation is reminiscent of traditional combinatorial design approaches for floppy mode systems [22]. However, here the site is continuously parameterized by $A_{ij}$, thus, the same design can encode a variety of matrices by designing the individual building blocks. This is a key to our ability to design and re-program the metamaterial.

Figure 3: (a) Evolution of the loss function during the optimization of the design from 1x1 lattices (black) to 6x6 (yellow). For large matrix sizes, the loss function cannot go below a saturation threshold. (b) Minimum aspect ratio required to accurately approximate a set of unitary random matrices. In the red curve, only $\theta_{3}$ is varied, while for the red curve, 5 angles are allowed to vary. (c) Input/output relation for a single unit cell, as a function of the angle of site 3, from a low (black) to high (yellow) amount of bending stiffness (black). We see that as the the aspect ratio is reduced, both the response at each given angle, and the maximum attainable matrix coefficient are lowered. (d) 3D model of a unit cell

The unit cell design consists of nine sites (Fig. 2a,b) that can move on a plane. Before introducing any constraints, the sites are free to move along both $x$ and $y$ axes, and thus the constraint-free unit cell has 18 floppy modes. As a consequence of the Maxwell-Calladine theorem, restricting the unit cell to the two prescriptive modes requires incorporating 16 linearly-independent constraints. To produce an initial guess for the metamaterial design, we start by modelling the constraints as ideal rods, for which bending stiffness is zero; and we consider the displacements infinitesimal, so as not to alter the length of the rods. We incorporate two classes of constraints: Local constraints (Fig. 2a)—represented by a pair of rods connected to a fixed boundary conditions, which restrict the site motion to a specific line along the $x$ and $y$ axes—introducing a rank-1 condition of the form $(\cos\theta,\sin\theta)^{T}(x,y)=0$; and pairwise constraints, which are realized by rods connecting two sites, resulting in a rank-1 condition of the form $(\cos\theta,\sin\theta)^{T}[(x_{1},y_{1})-(x_{2},y_{2})]=0$ between the displacements of the two sites.

The proposed unit cell design is introduced in Fig. 2b. The two floppy modes (Fig. 2c), spanning the two-dimensional input space, can be unterstood as follows: When the input $u_{1}$ is displaced, the rank-1 conditions introduced by the couplings between sites 2-3, 3-4 and 4-5 cause the output degree of freedom $v_{2}$ (horizontal displacement of site 5) to follow the input degree of freedom (satisfying the condition $v_{2}=u_{1}$ required by Eq. 1). At the same time, the condition introduced by the local constraints in site 3 couple its vertical displacement to its horizontal displacement, with a proportionality factor given by $\tan\theta_{3}$. This vertical displacement is then transferred to site 6, and converted to a horizontal displacement, scaled by the factor $\tan^{-1}\theta_{6}$, and transferred to the horizontal displacement of site 7. Similarly, the displacement $u_{2}$ is scaled by $\tan^{-1}\theta_{0}$ and $\tan\theta_{1}$, and then transferred through sites 1-4 and 4-7 to the vertical displacement of site 7. The diagonal coupling between sites 7 and 8 ensures that the displacement along the diagonal direction of site 8 matches the diagonal displacement of site 7. This diagonal displacement the projection on $1/\sqrt{2}$($dx_{7}$, $dy_{7}$), where $dx_{7}$ is scaled displacement of input $u_{1}$, and $dy_{7}$ is the scaled displacement of input $u_{2}$. This projection amounts to adding the horizontal and vertical displacements together. A local constraint on site 8 prevents the output site to move in the tangent direction to the 7-8 couplings. The output vertical displacement amounts to the projection of the diagonal displacement into the y-axis, introducing a total scaling of $1/2$. By setting the angles $\theta_{0}$ and $\theta_{1}$ to compensate for the factor of $\frac{1}{2}$, and setting the angles $\theta_{3}$ and $\theta_{6}$ to implement a coefficient $2A_{ij}$, the resulting system satisfies the required input-output characteristics.

The final system retains two floppy after the addition of the 16 constraints (7 local constraints plus 9 coupling constraints). These two resulting floppy modes are plotted in Fig. 2c. Since any prescribed input displacement can be expressed as a combination of these modes, in the ideal rod approximation, this (idealized) matrix-vector product requires no energy to be actuated. The floppy modes in our system are reminiscent to line-modes [26, 27, 28, 29, 23] observed in two-dimensional soft matter systems. However, in contrast with prior works, our unit cell design allows for line modes to branch and cross, as required for the realization of the matrix-vector product.

III Numerical validation

Until now, we have considered a design composed of ideal rods (zero bending stiffness) under infinitesimal displacements. However, these assumptions do not accurately describe the dynamics of experimental samples. To capture the effect of the finite bending stiffness of the rods, we construct a three-dimensional finite element simulation of the metamaterial (see Appendix A). Using automatic differentiation, it is also possible to compensate for some of these errors by fine-tuning the geometry.

Attempting to compensate for real-beam effects by adjusting the beam angles illustrates the limitations of the system (Fig. 3a): For small matrices, we can approximate the objective transformation with arbitrarily high accuracy. However, as the matrix size grows beyond $4x4$, the optimization algorithm cannot decrease the loss function below a given floor value. The origin of this error can be found in the relation between beam angle $\theta_{3}$ and matrix coefficients (Fig. 3c). While systems composed of ideal rods (infinite aspect ratio/zero bending stiffness) can approximate any matrix coefficient, rods with finite aspect ratio can only achieve a maximum matrix coefficient. This limitation arises from the ratio between longitudinal and bending stiffness in finite aspect ratio rods: Increasing rod angles has the effect of stiffening the input of the unit cell. Beyond a specific angle, the decrease in input displacement due to the finite longitudinal stiffness of the preceding rods exceeds the increase in output displacement due to the increase in angle. At this point, the matrix coefficient has reached saturation. We quantify this effect by determining the maximum matrix size that can be represented using a given aspect ratio. To do so, we consider 5 random matrices of each size, and perform an optimization to determine the minimum aspect ratio at which it is possible to achieve arbitrary accuracy. The results are plotted in Fig. 3b)—highlighting a linear relation between aspect ratio and maximum matrix size. Given that MEMS resonators with aspect ratios up to 1000:1 [30] are possible, we anticipate that our results can lead to matrix-vector products with dimensions above 64x64.

IV Experiments

Figure 4: (a) Complete experimental input/output transfer relation for a single unit cell, the measurements, in blue, are gathered from an optical flow algorithm. The fitted linear response (solid red lines) is within 20$\%$ of the target matrix value. The dashed olive line corresponds to a sigmoidal fit of the output displacement. (b) Picture of the experimental setup during characterization. The blue components are the output of the linear actuators. Videos of the experiment are provided as supplementary information. The red lines represent the deformation of the unit, it is computed using an edge detection algorithm and some edited to remove the edges of holes around the sample.

Figure 5: (a) input/output relationship of a 2x2 lattice. The matrix of the transformation is within 15% of the target matrix (b) Picture of the 2x2 lattice, with the deformation plotted in red. The first output is zero because the components cancel each other out, while they add for the second output.

We fabricate and optimize the design on a rubber sheet with a thickness of $6\,\text{mm}$, using water-jet cutting. The input displacement is prescribed using stepper motors, and the output displacement is measured with a camera and processed using an optical flow algorithm. The experiments on a single unit cell reveal that, in the linear regime (where input displacements remain below a few mm), the displacement transfer matrix coefficients between input and output are within 20% of the target value from the ideal unit. Beyond the linear regime (around 2 mm—for a unit of 250x275mm), the output displacement saturates. Interestingly, the saturation function resembles a sigmoid function. In neural networks, sigmoid functions are one of the most commonly used nonlinear activation functions—thus, the saturation characteristics of the material present an exciting prospect for the realization of full machine learning models in elastic systems. In addition to the sigmoidal-like saturation, the experimental results present two deviations from the ideal response: First, we observe that the response presents some hysteresis—which we attribute to the visoelastic response of the polymer material and non-idealities in the actuator-sample couplings. Second, the linear region is asymmetric—saturation occurs earlier under compressive loads than under tensile loads. Experiments on a metamaterial consisting of two unit cells reveal a similar response—with $12\%$ error between the objective matrix and the measured linear transformation, and similar hysteresis and saturation characteristics.

Figure 6: (a) Picture of the sample during characterization, with the lines showing a compliant mechanism in the low-stiffness (green) and high-stiffness (red) configurations. In the low-stiffness configuration, the beams are buckled, while in the the high-stiffness configuration, the beams are straight. (b) Experimental results demonstrating the switching between three different matrix coefficients, $A_{ij}=+.25$ (green), $A_{ij}=0$ (purple, red), $A_{ij}=-.33$ (blue). Switching between green and blue curves does not alter the number of floppy modes. (c, d) Variable stiffness mechanism in the clamped position (c), where the tension in the upper beams hold the shaft in place and (d) in the ’free’ position the buckling beams allow the shaft to move

We now proceed to demonstrate a programmable metamaterial that can be re-configured to implement different matrix-vector products. Reconfigurable floppy modes are an timely area of research, with existing results focusing on controlling the number of modes [31, 32]. The requirements for a programmable matrix-vector product are different, as the number of modes shall remain constant while the mode shape is altered. To allow for multiple possible values for the transformation coefficients, we incorporate additional rods in the site 3, the site that defines $A_{ij}$, the coordinate transformation. Increasing the number of constraining rods causes the system to be over-constrained—thus eliminating the zero modes. To achieve programmability, only a sub-set of the rods is active in each configuration. We accomplish this by incorporating bi-stable elements into the metamaterial [33]. In particular, we use a variable stiffness compliant mechanism [34], connected to the base of each of the support rods (Fig. 6a). When the compliant mechanism is in the low-stiffness configuration, the rods are free to move, and the constraint does not act into the system. By switching the compliant mechanisms between different configurations, the unit cell coefficient $A_{ij}$ can be switched between three different values (Fig. 6b)—although only two of them present the correct number of floppy modes. The programmable unit cell shows an increased hysteresis due to slack in the fabricated compliant mechanism. There are several avenues to increase the resolution in the matrix coefficients, including using more than four switchable support rods, or cascading multiple unit cells by cascading programmable units.

V Discussion

The present work demonstrated a mechanical metamaterial capable of re-programmable matrix-vector multiplication. The results can be interpreted in the light of the in-materia computing paradigm, where non-traditional material responses—here, floppy modes—are harnessed to perform computations—in this case, matrix-vector products. The key factor limiting the matrix size is the aspect ratio of the supporting beams, with matrix sizes above 64x64 attainable within present technology. Since, relevant machine learning such as phoneme recognition can be efficiently encoded in classifiers acting on a as little as 39 features [35], the proposed passive matrix-vector multiplier can have an enabling role towards passive speech-recognizing MEMS [7]. We have also shown a reprogrammable unit cell where every weight $A_{ij}$ can be adjusted from a pair of values while preserving the number of zero modes. This level of programmability—that can be expanded using more complex designs—has the potential to be applied to relevant problems, as machine learning models are known to be resilient to extreme weight quantization, down to 3 values per weight [36].

VI Appendix

VI.1 Simulation and optimization

We simulate and optimize the metamaterial using the finite element package FEniCSx [57, 56, 58, 59]. To do so, we partition the structure into sub-components 7. Each sub-component is obtained by extruding a two-dimensional geometry, composed of lines and circle arcs. The lines and arcs are generated using Python functions, from a set of high-level parameters, such as centers, angles, radiuses and thicknesses. We use JAX to map the high-level parameters into the elementary features (arcs and lines) to allow for automatic differentiation. Automatic differentiation is used to assemble the geometry from a graph of interconnected components. The graph specifies the topological characteristics of the system—i.e., which component is connected to which. Then, by computing the Jacobian of the distances between component boundary contact point with respect to the high-level parameters, $J_{contact}$, a distance minimization algorithm is capable of automatically constructing the full metamaterial geometry—enforcing continuity between components.

Figure 7: (a) Division of the planar geometry into sub-components (each component is plotted in a different color). The solid dots represent the contact points—that are automatically matched between each component to ensure compatibility. (b) finite element meshing of a single component, red dots represent boundary degrees of freedom, which match between components, which are deterministically positioned and thus guaranteed to match between components. (c) Gradient of the external mesh node coordinates with respect to a high-level parameter (in this case, the radius) (d) Gradient of the misfit $\eta$ with respect to the mesh node coordinates.

To simulate the system, we extract the stiffness matrix for every component using second-order Lagrange elements and a linear elastic energy function. Since we are only interested in the degrees of freedom located on boundaries; we use the Guyan reduction [45]. This computation allows us to determine the exact, statically condensed component stiffness matrix, that contains only a few hundred rows and columns. The determination of the statically condensed stiffness matrix, requiring a matrix inversion, is performed on the GPU—we found out that, given the relatively small size of each component (a few thousands of degrees of freedom), GPU-based dense linear algebra achieved the highest performance. The matrices for each component can be calculated fully in parallel.

To be able to assemble a full system matrix from the reduced component matrices, we make sure that all connecting boundaries are meshed with identical, regular grids. Then, we assemble the statically condensed stiffness matrix from every component into a global stiffness matrix. This global stiffness matrix can also be decomposed in inner degrees of freedom (corresponding to components that are connected to each other (not accessible from outside the metamaterial), and external boundaries (that correspond to the system inputs and outputs). Since we are interested in the response of the metamaterial to displacements in the external boundaries, we apply a second Guyan reduction. This second Guyan reduction can be very expensive for large systems. Thus, it is performed recursively, assembling intermediate structures with a larger and larger number of components until the full system is assembled.

Once the global stiffness matrix has been statically condensed to remove internal boundaries, we divide its degrees of freedom into input and output displacements,

$$K_{global}(\vec{\theta})=\begin{pmatrix}K_{in,in}(\vec{\theta})&K_{in,out}(\vec{\theta})\\ K_{in,out}^{T}(\vec{\theta})&K_{out,out}(\vec{\theta}),\end{pmatrix}$$ (2)

where $\vec{\theta}$ is the vector of high-level geometric parameters.

The transfer function between input and output displacements can be calculated as $T(\vec{\theta})=-K_{out,out}^{-1}(\vec{\theta})K_{in,out}^{T}(\vec{\theta})$. Once this transfer matrix $T(\vec{\theta})$ has been calculated, we prescribe the inputs to move along a specified direction—given by the displacement direction of the actuators, and project the resulting output into the direction of measurement. The projection $A(\vec{\theta})=\hat{P}_{out}T(\vec{\theta})P_{in}$, where $P_{in}$ and $P_{out}$ are the input and output projectors, provides the effective implemented matrix $A$. In our case, $A$ should match the objective matrix of the design, $T$. The agreement between the implemented and objective matrices can be quantified by defining a loss function of the form

$$\eta(\theta)=\sum_{i=0}^{N_{in}}\sum_{j=0}^{N_{out}}(T-A(\theta))^{2}$$ (3)

The entire evaluation tree, from the component stiffness matrices to the loss function $\eta(\theta)$ can be automatically differentiated using JAX. Once the change in misfit with respect to the component stiffness matrices has been determined, the gradient with respect to the finite element mesh can be calculated in FEniCSx using the built-in automatic differentiation capabilities in the UFL library [51]. This produces a gradient of the misfit with respect to the mesh coordinates (Fig. 7d), indicating how the mesh should be deformed to improve the matrix-vector multiplication accuracy. This mesh gradient can be projected into a gradient of the boundary coordinates with respect to the high-level parameters (Fig. 7d), producing the gradient of the misfit with respect to the high-level parameters. It should be noted that not all parameter updates are valid; some changes to the geometry will result in disconnected components. To prevent this, we project the gradient into the nullspace of the contact point distances $J_{contact}$ with respect to the high-level parameters, and perform an additional boundary distance minimization optimization to preserve the well-connectedness of the metamaterial.