Validating an algebraic approach to characterizing resonator networks

Resonator networks are ubiquitous in natural and engineered systems, such as solid-state materials, electrical circuits, quantum processors, and even neural tissue. To understand and manipulate these networks it is essential to characterize their building blocks, which include the mechanical analogs of mass, elasticity, damping, and coupling of each resonator element. While these mechanical parameters are typically obtained from response spectra using least-squares fitting, this approach requires a priori knowledge of all parameters and is susceptible to large error due to convergence to local minima. Here we validate an alternative algebraic means to characterize resonator networks with no or minimal a priori knowledge. Our approach recasts the equations of motion of the network into a linear homogeneous algebraic equation and solves the equation with a set of discrete measured network response vectors. For validation, we employ our approach on noisy simulated data from a single resonator and a coupled resonator pair, and we characterize the accuracy of the recovered parameters using high-dimension factorial simulations. Generally, we find that the error is inversely proportional to the signal-to-noise ratio, that measurements at two frequencies are sufficient to recover all parameters, and that sampling near the resonant peaks is optimal. Our simple, powerful tool will enable future efforts to ascertain network properties and control resonator networks in diverse physical domains.


I. Introduction
Resonator networks are a ubiquitous [1] and diverse class of systems, found in both natural [2] and engineered contexts.They can range in scale from astronomical systems [3] to biological and neural networks [4][5][6], and from small-scale micromechanical lattices to large integrated circuits [7] and solid-state materials.These many-body systems exhibit rich collective behaviors, such as brain memory, quantum [8,9] and classical computation, and optical properties of solids, making them an important subject of study.To understand and engineer this behavior, it is necessary to characterize the network building blocks (e.g., neurons, micromechanical resonators, ions, etc.) and their connectivity.Despite many realizations, a useful and simple model for these networks is as a collection of coupled mechanical mass-spring resonators (Fig. 1a).Using this model, the building blocks are defined locally by the elasticity, mass, and damping of each resonator, while the connectivity is captured by coupling springs.For a linear response, the resonator network is governed by the equation of motion where , , and  are the mass, damping, and elasticity matrices, respectively, and  ⃗ is the external force.Typically, the mechanical elements (, ,  , and  ⃗ ) of these systems are characterized by analyzing amplitude and phase spectra with non-linear least squares (NLLS), where the mechanical elements are fitting parameters.NLLS fitting requires both solving the coupled differential equations (i.e., determining a closed-form solution for  ⃗()) and a priori knowledge in the form of initial guesses for each of the network's mechanical parameters.These initial guesses are crucial, since the NLLS solution may converge to a local minimum determined by the initial guesses, and poor initial guesses can lead to inaccurate estimates of the true state of the network [10].Moreover, as the number of parameters increases, the optimization problem becomes more complex, and the number of potential local minima generally increases, making NLLS more inaccurate [10].
Recently, we developed an alternative, non-regressive algebraic approach to circumvent the issues of NLLS [10], a method we call Network Mapping and Analysis of Parameters, or NetMAP [11].In this approach, we transform Eq. (1) into where () = − 2  + i + ,  ⃗ is the amplitude of  ⃗ , and  ⃗ () is the complex response vector of the network which includes the amplitude and phase for each resonator in the network measured at a frequency  (Fig. 1b).We then measure and construct  ⃗ () at two or more values of , and use these measurements to rearrange Eq. ( 2) into an augmented homogeneous equation where  is a matrix of known, measured quantities determined by  ⃗ ().The vector  ⃗ consists of the desired, unknown elements of , , , and  ⃗ (Fig. 1b, and see Methods), and the solution space for  ⃗ is the null-space of .In contrast to NLLS, NetMAP does not require a priori knowledge (i.e.neither exact nor approximate knowledge of the elements of , , or ) nor iterative computation, but instead solves for  ⃗ directly with as few as two response vector measurements ( ⃗ ()).In Ref. [11] we used NetMAP to determine  ⃗ for small clusters of graphene NEMS resonators and we found excellent agreement with expected parameter values and broader spectral response.We calculate the relationship between the complex amplitude  ⃗ () of the masses and the underlying parameters.c) Steps of the validation process: 1) Modeling a mass-spring network, such as a monomer, shown here, or dimer.2) Setting the parameters  ⃗  to simulate.3) Solving the equations of motion (EOM) and simulating spectra with noise.4) Using simulated measurements of amplitude and phase at discrete frequency points to construct a matrix.5) Using Singular Value Decomposition (SVD) to recover the parameters  ⃗ ̂, and scaling the result as needed.Here  is shown without a hat to identify it as a known quantity while the others are scaled to . 6) Calculating the percent error  to compare the recovered parameters to the set parameters.7) Calculating the expected spectra for the recovered parameters, and calculating  2 to compare the simulated data and the expected curve.We find that  2 is correlated with error .
While NetMAP is a promising new technique, there are many unanswered questions regarding its accuracy in predicting the mechanical parameters of the network: How does the number  of response vectors, their signal-to-noise ratio (SNR), or the frequency at which they were measured affect the accuracy of NetMAP?How do the actual values of the mechanical parameters or the dimension of the solution space of Eq. ( 3) affect the accuracy?To answer these fundamental questions, we simulate noisy response vectors for single-resonator (monomer) and coupled-pair (dimer) network clusters using predetermined mechanical parameters and then we use NetMAP to predict these parameters.We find that the accuracy of NetMAP predictions improves with the SNR, as does measuring response vectors near spectral peaks.While measuring additional vectors moderately improves accuracy, we find a minimum of two response vectors is sufficient, thus requiring a remarkably small number of measurements.Moreover, the dependence of the accuracy on the null-space dimension is nuanced; the accuracy of a 1D null-space solution varies with the input parameters.

II. Methods
To assess the accuracy of NetMAP, we statistically compare the input parameters to the recovered parameters and the actual simulated spectra to the expected spectra from the model.The steps of the validation process are shown in Fig. 1c (see Supplemental Material i for details.)In the first step, we choose the size of mass and spring network to model ( cluster ).In this work, we modeled the monomer (a single resonator,  cluster = 1) and the dimer (two coupled resonators,  cluster = 2).We then set the network parameters  ⃗ in to fixed numerical values (step 2) and input them into the analytical solution of the equations of motion (EOM) to obtain an exact, noise-free complex response function for each resonator,   () (see Supplemental Material).The response vector  ⃗() solves Eq. ( 2) exactly and has  cluster complex components.
To simulate a random experiment (step 3), we add real and imaginary noise to each response component to obtain   (, ) =   () + Γ , (, ) + iΓ , (, ), where Γ(, ) is a pseudorandom number drawn from a zero-centered normal distribution of variance  2 and regenerated for each frequency .Next (step 4), we select a set of noisy response vectors  ⃗ (, ) at  discrete frequencies (e.g.  and   in Fig. 1c, so  = 2).Each  ⃗ () provides two vector equations corresponding to the real and imaginary parts of Eq. ( 2), for a total of 2 cluster linear equations, which we use to populate the matrix elements of .The matrix  has dimensions 2 cluster × , where  = dim ( ⃗).(See Supplemental Material for general closed form of  for the monomer and dimer.) To obtain the predicted parameters vector  ⃗ ̂ (step 5), we use the simulated  to solve  ⃗ = 0 ⃗⃗ with a computational Singular Value Decomposition (SVD) algorithm.An algebraic approach is not standard but is straightforward for describing resonator systems.NetMAP uses singular value decomposition (SVD), an algebraic approach to solving systems of equations with applications in medical imaging [12], antenna arrays [13], planar transmission lines [14], and movie recommendations [15].SVD is similar to eigensystem solvers, with singular values instead of eigenvalues and singular vectors instead of eigenvectors, while allowing the matrix  to be rectangular rather than square.In using SVD to solve  ⃗ = 0 ⃗⃗ , we seek a singular value of zero and its corresponding singular vector, which is a solution for the physical parameters  ⃗.The solution space for non-trivial  ⃗ ̂ will be at minimum onedimensional (1D), but may have higher dimension.The solution-space dimension provided by SVD is open to interpretation, but commonly determined by the number of singular values with  ≪ 1.We explore the accuracy of NetMAP with 1D, 2D, and 3D solution spaces.When we define higher-dimension solution spaces, we use the smallest  and their corresponding singular vectors  ⃗ ̂ in order of increasing value.In the 1D case, we scale  ⃗ ̂ so that the force  ̂ equals the input force .For higher-dimensional solution spaces, we supply additional parameter constraints (see Supplemental Material).
Finally, to test the NetMAP parameter predictions, we replicate steps 3-5 of the simulation for a total of 1000 runs to generate a sample distribution for  ⃗ ̂, and then use a one-sample statistical -test, to quantify the agreement between the  th predicted network parameter ̂ and the corresponding input parameter  ,in .For each  ⃗ ̂ from the sample, we also compute the fractional error for each parameter   = |̂ − ( in )  |/( in )  , and we calculate the correlation coefficient (  2 ) between the expected spectra   ̂() and the simulated noisy spectra   ().

III. Case studies A Monomer
To demonstrate NetMAP's algebraic approach, we first consider a test case with a lightly damped monomer.For a general monomer, the noisy spectrum is given by For this test case, we arbitrarily choose  ⃗ in with parameters: mass  = 4 kg, damping coefficient  = 0.01 N/(m/s), spring stiffness  = 16 N/m, and force amplitude  = 1 N.Moreover, we set the input error for this demonstration to  = 5 × 10 −3 m.Fig. 2a shows a simulated spectrum for the amplitude () and phase () of (, ); Fig. 2b shows the real and imaginary parts of the same (, ) plotted in the complex plane.To obtain  ⃗ ̂, we use (, ) evaluated at the two frequencies   and   shown in Fig. 2a,b, which we then use to construct .For a monomer sampled at two frequencies, the equation  ⃗ = 0 ⃗⃗ written explicitly is: where () = Re(()) and () = (().The simplest way to build  is to read the values from complex plane plots, as in Fig. 2b.By solving Eq. ( 4) with SVD, we obtain solutions for the physical parameters  ̂,  ̂, and  ̂.The observed -statistics for 1D, 2D, and 3D solution spaces is shown in Table 1.Given our large sample size (1000), the  0 approximates the number of standard error intervals the predicted value is from the expected.All -values exceed the  = 0.05 by at least a factor of 10, so we conclude all NetMAP predicted parameters for this trial agree with the set input values  ⃗ in .2c.The 1D result recovers all parameters within 0.04% of the input values in 95% of the trials (Fig. 2c).The 2D discrepancy results are similar or better than the 1D; for example, the 95% confidence range for the elasticity is ~2.5 × 10 −5 %.The 95% confidence interval for the mean Δ  /( in )  are much tighter; for the 1D mass, The fractional error,   = |̂ − ( in )  |/( in )  is a useful metric to assess the accuracy of NetMAP.However, an experimentalist without knowledge of input parameters cannot calculate the fractional error but they can calculate the coefficient of determination  2 .To elucidate the relationship between the average  2 and   , we plot the average error 〈〉 = 1 − ∑   versus 1 −  2 for each of the 1000 simulated trials and for the 1D and 2D null-spaces (Fig. 2d).For both dimensions, we observe a linear correlation ( 1 2 = 0.67,  2 2 = 0.4).For the 1D case, we observe 〈〉 decreases as √1 −  2 ( 1 2 ~0.74,  ≪ 0.001).The correlation between 1 −  2 and error provides a means for an experimentalist to assess the accuracy of recovered values without a priori information.The corresponding distributions for error 〈〉 are provided in Fig. 2e.While the 2D error is lower than the 1D (〈〉 1 = 0.0144 ± 0.0003% vs. 〈〉 2 = 0.0039 ± 0.0001%), we had to specify two input parameters ( in )  to solve for  ⃗ ̂.In general, for a D-dimensional null-space,  number of input parameters  ⃗ in are required to obtain a solution, which is a disadvantage in terms of a priori knowledge.Generally, 1D solutions are preferable and sufficient to recover the parameters with low error..The remaining () [colorful datapoints] are used only for validation, not for finding the recovered parameters.c) Box and whisker plot showing the spread in recovered parameters,   / ,in = (̂ −  ,in )/ ,in , for 1D and 2D solutions with 1000 runs of the same parameters.For 2D-SVD, additional a priori information is required: we fix  =   in order to select the solution from the 2D solution space.d) Error ⟨⟩ and 1 −  2 are correlated (1000 runs).We calculate ⟨⟩ from a priori information, but we may predict it approximately from the  2 value.e) Histogram and box plots of the percent error ⟨⟩ (1000 runs), showing a half-normal distribution.f) Expanding to a range of input noise  for this lightly damped monomer (80 runs per ), we find that the 2D-SVD solution is slightly more accurate than the 1D-SVD solution.The 3D solution is many orders of magnitude less accurate.The vertical grey line shows the input noise  = 0.005 m appearing in other subfigures (a,b,c,d,e) of this figure.The average percent error and the standard deviation are related by a power law.The average curve is calculated as a mean of the logarithm of the average errors and the shaded regions indicate 95% of the runs.
So far, we have presented results for a fixed input noise ( = 5 × 10 −3 m.).To determine how the level of noise affects the error, we sweep  through several orders of magnitude and run the simulation 80 times per noise value.The error 〈〉 vs.  is shown in Fig. 2f, with the  = 5 × 10 −3 m trial indicated with a vertical gray line.We also compute the signal-to-noise ratio defined as SNR =   and add it to the upper axis.We find that 〈〉 for all solution space dimensions varies approximately linearly with  (e.g. for 1D, 〈〉 =   with  = 1.0004 ± 0.0016).The 2D-SVD error is the lowest, but is followed closely by the 1D error.The 3D error is several orders of magnitude larger than either 1D or 2D. .e) Box and whisker plot showing the spread of recovered parameters as fractional discrepancy, Δ  / ,in = (̂ −  ,in )/ ,in , for 1D, 2D, and 3D solutions over multiple trials.For 2D-SVD, an additional parameter,  1 , is fixed at  1, in order to identify the solution within the 2D solution space.For 3D-SVD,  2 is also fixed.f) Histogram and box plots of the average error ⟨⟩.g) The average error ⟨⟩ is correlated with 1 −  2 .h) High SNR is key to minimizing the error.The error ⟨⟩ ̅̅̅̅ is proportional to the input noise .

A Dimer
We now similarly analyze a two-mass (dimer) resonator system (Fig. 3).For this system, we set the input simulation parameters as follows:  1 = 1 kg,  2 = 10 kg,  1 = 1 N/m,  2 = 10 N/m, coupling spring  12 = 1 N/m,  1 =  2 = 0.1 N/(m/s),  1 = 10 N, and noise  = 0.005 m.The force is applied to  1 .This case represents an anti-crossing scenario where the coupling splits the resonant frequency degeneracy.The noisy spectra of each resonator are shown in Fig. 3a-d, and we identify the two frequencies  and corresponding responses (black circles) needed to construct  and solve for  ⃗ ̂.The fractional discrepancy Δ  /( in )  for 1D, 2D, and 3D solutions is shown in Fig. 3e (see Supple-mental Material for -test results); we also include the discrepancy for the isolated resonance frequencies   = √   /  .We find that the damping of the undriven resonator,  2 , is the least accurate of the recovered parameters.Nonetheless, all parameters are recovered with an error standard deviation below 0.5%.The distributions for average fractional error 〈〉 are provided in Fig. 3f, where 〈〉 1 = 0.123 ± 0.002%, 〈〉 2 = 0.060 ± 0.001%, and 〈〉 3 = 0.060 ± 0.001%).Using  ⃗ ̂, we calculate and plot the recovered spectra  ̂1() and  ̂2(), shown as black dashed lines in Fig. 3a-d.As before, we plot the average error 〈〉 versus 1 −  2 (Fig. 3g) and the input noise  (Fig. 3h), and we observe a similar correlation with 1 −  2 and functional dependence 〈〉 ∝  for all null-space dimensions (see Supplemental Material).

IV. Optimizing frequency selection
We have thus far demonstrated NetMAP on a monomer and a dimer system using the minimum number ( = 2) of response vector measurements.Even in the presence of noise, the 1D null-space solution recovers all network parameters with high accuracy.While any choice of at least two frequencies for response vector measurements will allow the analysis to proceed, we expect some frequencies to yield better results than others, and that sampling response vectors at more frequencies will yield different results.To test these ideas, we varied the number and value of the sampling frequencies for both a monomer and a dimer system while characterizing the accuracy of NetMAP.

Monomer frequency optimization.
In order to assess the accuracy as a function of the number of measurement frequencies , we first consider a moderately damped ( = 16) monomer with  = 4 kg,  = 0.4 N/(m/s),  = 10 N/m,  = 1 N,  = 5 × 10 −3 m and we sweep  from  = 2 to 25.We start with two points near the resonance frequency  res = 1.58 rad/s and add additional points sequentially to the right (i.e. higher frequency) and left (lower frequency) of the peak, with Δ = 0.01 rad/s separation between measurement points (Fig. 3a, see Supplemental Material for the amplitude and phase spectra.).We replicate a given  for a total of 100 runs.The SNR for the peak datapoint is / = 3164.Fig. 3b shows the average error 〈〉 vs.  for 1D, 2D, and 3D null-space solutions, which we label as 1D-SVD, 2D-SVD, etc.The 1D error gradually decreases with , but with diminishing returns as the number of measurement points becomes large; the 1D mean errors for  ≥ 7 are statistically equivalent (see Supplemental Material for ANOVA post hoc analysis).The 2D error dependence on  is similar the 1D, but plateaus after  = 4.However, while the 2D solution has lower error than the 1D for low , the two solutions appear to converge for  ≥ 9; at  = 25 the mean difference between 1D and 2D error is 0.008 ± 0.002%.As with the first monomer test case, the 3D error is orders of magnitude larger than 1D or 2D error and does not decrease with  but instead oscillates with  as frequency points are incorporated to the left and right.Due to the high error of the 3D solution (sample distribution, 〈〉 3 = 8 ± 15%), it is not recommended for this system.
To investigate the effect of frequency choice on the error, we consider two response vectors with frequencies   and   and measure 〈〉 ̅̅̅̅ for 1D and 2D null-spaces as we vary the frequencies across resonance.We plot 〈〉 ̅̅̅̅ vs. (  ,   ) in Fig. 4c, and find that both the 1D and 2D solution are optimized near the peak resonance frequency ( res = 1.58 rad/s).The "cross" patterns indicate that if either   or   is on resonance, then the other frequency can take nearly any value and the error remains small.For the 1D case, however, there is a narrow diagonal line of high error (68.3% ± 1.2% standard deviation), indicating that   and   must differ to obtain a low-error 1D solution.This diagonal line is absent in the 2D case, so that replicated  ⃗ () with the same frequency are possible and yield low error.To observe more general trends, we plot the average error for different choices of   while Input (circled datapoints) and output (dashed black curve) spectrum for a simulated monomer with 25 circled points used for analysis.b) In order to compare the number  of frequency measurements, we plot the average error for m,k,b as a function of the number of frequency points used in the SVD analysis.The lines show the mean error ⟨⟩ ̅̅̅̅ of 100 simulated trials (shown as datapoints) for each number of frequency points.Violin plots show the distribution of ⟨⟩.For this monomer system, the 2D-SVD solution (orange) has the least error, then the 1D-SVD (blue), and the 3D-SVD (green) has the largest error.The average percent error of the parameters  ⃗ varies with the number of frequency points used in the SVD analysis.For this analysis, we add the frequency points in a specific order, starting with two points at each of the two amplitude peaks and sequentially adding points to the left and right of the peaks.The specific datapoints used for b are indicated in a. cdef) Comparing the choice of frequencies if there are  = 2 frequencies.c) The average error of the SVD solution varies with the choice of the two frequencies   and   .The 1D solution fails when the two frequencies are the same.Both the 1D and 2D solution have better results when the measurements are taken near the peak of resonance.d) The average error varies with the measured frequency   .Here the average is taken over all the results shown in (c), with   varying from 1.4 to 1.8 rad/s.e) For a double ( = 2 frequency) measurement of a monomer, the average error ⟨⟩ ̅̅̅̅ varies with the choice of frequencies measured.In this case, we fix   at the resonance peak and sweep   .We find that the optimum second frequency occurs when the phase at   is near   = − taking an average of all   values (Fig. 4d).This result shows that the optimal   to extract  ⃗ (  ) is near the resonant frequency,  res .Moreover, by fixing   =  res we determine which   will yield the lowest error.Fig. 4e shows the average error vs. the complex phase of the response vector  ⃗ (  ) for varying   .From this plot, the error of the 1D solution reaches a minimum near   = −3/4 and   = −/4 , and is greatest at   =   ≈ −/2 .In accord with Fig. 4c, the error of the 2D solution is minimum near   =   ≈ −/2.We see similar optimal phase choices for other monomer examples (see Supplemental Material).
For the monomer under study, the 2D null-space solution generally has the lowest error.However, increasing the dimension of the solution space assumes the singular value of the additional degree of freedom is sufficiently small, i.e.  ≪ 1.If this value is not small, the lower dimension null-space should have lower error because the singular vector corresponding to the large singular value does not solve the homogeneous Eq (3).To test this idea, we plot the average error for 1D, 2D, and 3D nullspaces against the value of the second smallest singular value ( 2 ) for each solution, as shown Fig. 4f (for this plot, we use the full set of swept frequency pairs from Fig. 4c.)We find that the 2D solution, which uses both the first and second singular vectors, is more accurate when  2 is small, while the 1D solution improves in accuracy as  2 grows.(See Supplemental Material for the 1D solution error as a function of both  1 and  2 .) The 3D solution additionally uses the singular vector corresponding to the third singular value ( 3 ≥  2 ).As expected, the 3D solution has consistently higher error than either 1D or 2D.
Our monomer case studies shed light on NetMAP best practices.For the above case study, the error generally decreases modestly as the number  of measured response vectors increase; before plateauing, the 1D error decreases by ~5 × with  = 7 and the 2D error decreases by ~1.5 × with  = 4.In either case,  = 2 suffices.The 1D solution is preferred to the 2D because it requires minimal prior knowledge and has reasonably low error (i.e.< 1%).However, the 2D error is consistently lower (mean of ~8.3 × across all ) if  1 ~2 ≪ 1.To obtain the lowest error with the 1D solution,  ⃗ () should be measured on resonance (  ≈ −/2) and /4 radians off resonance, where   = −3/4, −/4.The sampling frequencies from Fig. 2 were selected in this way.Finally, given the error is proportional to the noise, it is beneficial to perform noise-filtered, long-integration experimental measurements of (  ) and (  ).

Dimer frequency optimization
We now repeat the optimal-frequency analysis for a dimer system with the following input parameters:  5c.The error for 1D solutions starts at ⟨⟩ ̅̅̅̅ 1D = 0.50 ± 0.30% (standard deviation, SD) for  = 2, and improves with increasing  down to ⟨⟩ ̅̅̅̅ 1D = 0.125 ± 0.077% (SD) at  = 50.Analyzing the 1D case with ANOVA (see Supplemental Material), we see that solutions with  ≥ 16 have equivalent error and provide a modest ~2.5 × error improvement over the  = 2 solution.Compared to 1D solutions, the error for 2D solutions is higher, ranging from a minimum error of 2.37 ± 1.88% (SD) at  = 12 up to 215 ± 208% (SD) at  = 2.There is a pattern that repeats for every fourth additional frequency point: when a measurement is added near the stronger 0.774 rad/s peak, the 2D solution improves, and when a measurement is added near the 3.501 rad/s peak, the 2D solution becomes less accurate.Some of the measurements are numbered in Fig. 5 a,b to indicate the order in which they are included into Fig.5c.The 3D solution has the overall lowest error for 2 ≤  ≤ 11, with a minimum ⟨⟩ ̅̅̅̅ 3D = 0.052 ± 0.037% (SD) at  = 11, but then suffers an abrupt loss in accuracy for  ≥ 12 where the error jumps to 99 ± 79% (SD).To determine the optimal frequencies, we measure the average White plus signs indicate datapoints that are added as  increases from 2 to 50, and numbers indicate the order in which the points are included, with odd numbers for the larger peak and even numbers (not all shown) for the smaller peak illustrating how we alternate adding the frequency points.Only one resonant frequency,  res = 0.774 rad/s (blue points), appears in the spectrum of R1 (a), but a small second peak at  res = 3.501 rad/s (yellow points) appears in the spectrum of R2 (b).c) The average error varies with the number of datapoints  used in the analysis, with 1D-SVD in blue, 2D-SVD in orange, and 3D-SVD in green.For this analysis, we add the frequency points in a specific order, starting with two points at each of the two amplitude peaks and sequentially adding points to the left and right of the peaks.Simulated with 99 trials per each of 49 values of , totalling 4851 trials.d) If  = 2 datapoints are used, the 1D solution is most accurate when a measurement is taken at each of the two resonant peaks.e,f) The accuracy of the 2D and 3D solutions also varies with the pair of frequencies chosen.In general, it is less accurate for the two frequencies to be equal or nearly equal (dark diagonal lines).g) Replotting (d) in terms of R2 phase rather than frequency shows greater detail.Subfigures (defg) are simulated with 240000 trials.h) When   is fixed at the lower resonant peak,   = 0.774 rad/s, then the accuracy of the 1D solution is optimized when   is at the higher resonant peak, 3.533 rad/s (blue dip).The 2D solution is not accurate when each measurement is at one of the two resonant peaks (sharp orange peak).Shaded areas show 95% of simulated ⟨⟩ results and lines show ⟨⟩ ̅̅̅̅ .Simulated with 1600 trials per each of 199 values of   , totalling 318400 trials.
error with  = 2 and sweep both frequencies   and   across a range that covers both spectral peaks (Fig. 5d-f).Each solution space dimension has a high-error diagonal band corresponding to   =   .For the 1D solution, the lowest error occurs when one frequency is at the top of one peak and the other frequency is at the top of the other peak.The spectral peaks have the two highest SNR values that also sample both resonances.The 2D and 3D error patterns are more complex and more forgiving in terms of frequency choice; in both cases, it is sufficient to have one frequency near a spectral peak resonance, but the lowest error still occurs by sampling near each spectral peak.
To see the detail of the 1D case, we plot the Fig. 5d dataset as a function of complex phases of  2 (  ) and  2 (  ) in Fig. 5g (see Supplemental Material for phase spectra plots).The  2 () phase  2 () is one-to-one, unlike  1 (), so there is no ambiguity about which value of  2 corresponds to which value of .The phase plot shows consistently low error for quadrants II and IV, corresponding to recommended phase values  2 ∈ [0, ] for one measurement and  2 ∈ [, 2] for the other.In terms of the complex spectrum of  2 () (Fig. 5b), these phase values correspond to the upper and lower loops (spectral peaks).High error occurs in quadrants I and III, where the two response vectors are sampled on the same loop while the other is not sampled.Moreover, dark bands of high error occur when either phase is equal to 0 or , which correspond to response vectors near the origin of the complex plane with near zero SNR.We further examine the behavior of the 1D and 2D solutions near the spectral peaks by fixing   = 0.774 rad/s and sweeping   , as shown in Fig. 5h.For the 1D case, error decreases as   → 3.533 rad/s, which corresponds to the higher frequency spectral peak, where it reaches a minimum of 0.5%.The 2D case has similar behavior, but spikes to a maximum over 100% when   = 3.533 rad/s, indicating a failure in the 2D null-space accuracy when the response vectors are sampled directly on resonance.Thus, for a dimer with two resonant peaks and a 2D solution, measuring near-but not on-the spectral peaks is ideal.

V. Accuracy varies with input parameters
So far, we have presented case studies of monomer and dimer resonator systems with fixed input parameters,  ⃗ in .However, it is possible that the error behavior we observe depends on the choice of  ⃗ in .To probe the dependence of the error on the input parameters themselves, we run a full, replicated 2  factorial experiment for a general dimer system by varying the mechanical parameters (see Methods).We screened and ranked the mechanical parameters as predictors of average error for a dimer solved using 1D-SVD; including up to two-factor interactions, the order from most predictive to least is: ,  1 ,  2 ,  12 , and the  1 ⋅  2 interaction (see Supplemental Material for full effects model results).We previously discussed using 1 −  2 and  2 as indicators of average error (Fig. 3g and Fig. 4d).
To gain a deeper understanding of how the parameters  ⃗ of the resonant system affect the SVD accuracy, we plot the average error as a function of varied parameters in a few example dimer cases (Fig. 6).Each cartoon in Fig. 6 shows one parameter in orange that we vary while keeping all others fixed.We identify the two peak frequencies  res using a peak-finding function [16] and plot them below each cartoon to illustrate their variation with the parameter.To improve the accuracy of the NetMAP results, we choose 4 additional frequencies on each side of the two peak frequencies  res for a total of  = 10 frequencies for analysis.To demonstrate the impact of signal to noise ratio (SNR) on the error as we vary each parameter, we plot the  res datapoints with higher SNR in orange, highlighting how the error falls as SNR increases.Increasing the driving force (Fig. 6a) or the non-driven mass  2 (Fig. 6b) causes the SNR to increase and the error in the recovered parameters to decrease correspondingly.This matches our findings that error is inversely proportional to SNR, and we generally find that more accurate recovered parameters result for dimer systems with high force amplitude and larger  2 .We observe the error increases as we increase the coupling spring stiffness  12 from 0 to 20 N/m (Fig. 6c), and this is usually true for dimers.For this dimer system, a lower  1 corresponds to more accuracy (Fig. 6d), but this is not generalizable to all dimer systems.As  2 is varied in Fig. 6e, the resonance peaks show an anticrossing near  2 = 20 N/m, and the 1D solution loses accuracy near the anticrossing.This is a specific result for this particular dimer system.To see how an anticrossing affects a dimer system with different input parameters, we consider the system shown in Fig. 6f, with an anticrossing near  2 = 7 N/m, and find, in contrast, that the 1D solution has higher accuracy near the anticrossing, suggesting that anticrossings may increase or decrease the error and the variation in error is likely associated with the SNR.Furthermore, to illustrate the effect on accuracy when the peak-finder fails to identify one of the resonance peaks, we allow the peakfinder to miss the higher  res peak for  2 = 10 to 12 N/m, and find that the 1D solution loses accuracy while the 2D solution gains accuracy in that range.Thus, Fig. 6 shows a detailed view of how the parameter values affect the accuracy of the algebraic solution.
The frequency optimization and factorial studies provide broader NetMAP best practices for the dimer system.Driving with a higher amplitude force improves accuracy.For our case study in Fig. 5, the error decreases modestly as the number of frequency points increases, though  = 2 measured response vectors are sufficient.As it was for the monomer case, the 1D solution is preferred to the 2D or 3D because it requires minimal prior knowledge of the parameters, while still having reasonably low error, and because the 2D and 3D solutions have erratic accuracy, depending heavily on the particular frequency points sampled in our example in Fig. 5.To obtain the lowest error with the 1D solution,  ⃗ () should be measured at each of the two resonance peaks of the dimer spectrum.The two sampling frequencies in Fig. 3 were selected in this way.As with the monomer case, the error is proportional to the noise, so it is beneficial to perform noise-filtered, long-integration experimental measurements of (  ) and (  ).

VI. Discussion
For an experimentalist, NetMAP offers a means to analyze amplitude and phase data in order to reveal the physical parameters for each individual oscillator and the coupling between the oscillators.This phase-sensitive data is standard for lock-in amplifier measurements.However, for NetMAP to be useful, the experimentalist must measure the amplitude and phase for each oscillator in the network, which in some cases may be challenging.The SNR range we describe here (e.g. 10 3 , or 30 dB SNR) is attainable to experimentalists [11] working with a wide variety of systems, including electronic or micromechanical oscillator networks.
In cases where NetMAP returns an inadequate  2 value, it may be advantageous to combine NetMAP and NLLS fitting by using the results from NetMAP as the initial guesses for an iterative solution, and thereby improving  2 .Since  2 correlates with the error in the parameters, this is expected to improve the accuracy of the results.
We have shown how well NetMAP performs for a monomer (single mass) and dimer (double mass) system, but further investigation is needed to evaluate its performance in larger systems.Our calculations (see Supplemental Material) show that, even as the number of oscillators increases, only two discrete frequencies are needed for the required response vectors  ⃗ .For a larger system, the size of the response vectors would scale with the number of oscillators, increasing the size of the matrix , which may require solution spaces of higher dimension.Future work will explore NetMAP performance in larger systems.
Γ , (, ): A pseudorandom number generated from a Gaussian distribution (standard deviation ) to add noise to the real part of   for the  th resonator at a given discrete frequency .
Γ , (, ): A pseudorandom number generated from a Gaussian distribution (standard deviation ) to add noise to the imaginary part of   for the  th resonator at a given discrete frequency .
angular frequency.The resonator network response frequency is assumed to be equal to the driving frequency for the steady state solution.  ,   or  1 ,  2 , … ,   : input frequencies, the set of frequencies selected for analysis.
: number of trials (number of repeats + 1) = each resonator, from 1 to  cluster (e.g.  ,   ,  , 2 ,  , 2 ) i = imaginary constant  ⃗: parameters vector.It includes all unknown parameters and possibly 1 or more known parameters for scaling the solution.The order of the parameters in the parameters vector must be consistent with the order of the columns of .
⃗ in : input parameters.This is the a priori information.
⃗ ̂: the parameters output by NetMAP (Network Mapping and Analysis of Parameters).Components of this vector are  ̂1,  ̂1,  ̂1, et cetera.
: number of trials for a factorial experiment R1: resonator 1.The resonator we are driving with an oscillating force.
R2: resonator 2 in a two-mass (dimer) system.The resonator we are not directly driving.

Extended description of NetMAP
Resonator network theory A system of equations for a resonator network can be organized into the matrix form: where  ⃗ is the displacements for each resonator,  ⃗ is the oscillating forces driving each resonator, and , , and  are matrices with all information about the inertia, damping, and elasticity of the network.If the force vectors are all oscillating with driving frequency , then we guess an oscillating steady-state solution  ⃗ =  ⃗  i where  ⃗ is a constant complex amplitude vector with components   =    i  , and obtain where  ⃗ is the amplitude of the force vector,  ⃗ =  ⃗  i .More compactly, the equation of motion system can be written where the symmetric matrix () = − 2  + i +  contains complete information about the inertia, elasticity, and damping of the network.

Solving the equations of motion
In order to simulate spectra to validate NetMAP, our approach is to set the input values  ⃗ in , then calculate the spectra () with added random Gaussian noise, thereby simulating an experiment (steps 2 and 3 of Fig 1c).To calculate the input spectra, we implement Cramer's rule in SymPy [1] to symbolically solve Eq. (S2) given the input values  ⃗ in .This yields the expected response vector  ⃗(), where each vector component   corresponds to the th resonator.To simulate noisy measurements, we add noise from a random Gaussian distribution to each vector component:

𝑍 𝑖 (𝜔) = 𝑧 𝑖 (𝜔) + Γ 𝑥,𝑖 (𝜎, 𝜔) + iΓ 𝑦,𝑖 (𝜎, 𝜔). (S4)
We add noise to the real and imaginary parts of   rather than to the amplitude and phase because doing so better matches experimental results, especially at low amplitudes, where the phase becomes highly uncertain.We also use the same solution to the equations of motion to calculate the output spectra  ̂() for step 7 of Fig 1c, but we use the output parameters  ⃗ ̂ and do not add simulated noise.
While the method we present here would be appropriate for any resonator system, we consider as a particular case a micromechanical system where we measure the amplitude and phase of each resonant mass using scanning interference measurements (SIM).[2,3] Determining the number of frequencies required for NetMAP

Clusters
Consider the case of a linear chain of resonators (Figure S1a) where we apply an oscillating force of amplitude  to one resonator in the chain.For a linear chain of resonators (Figure S1a), with coupling springs of stiffness   only between nearest neighbors, the matrix  is tridiagonal (Figure S1b, left), with diagonal elements (yellow squares) of the form −   2 + i   +  −1 +   +   and off-diagonal elements (blue squares) of the form −  , corresponding to Eq. (S2).
In an experimental system that is driven at one resonator (Figure S1b, right, light green square), a finite number of resonators will respond, indicated by the shaded squares for  ⃗ (), which selects a relevant submatrix  ̃() (Figure S1b, purple square).Our response vector approach determines elements of  ̃() from the observations of  ⃗ () using Eq.(S3).However, the number of unknown elements in  ̃() and  is larger than the number of equations provided by a measurement of  ⃗ () at one frequency, thereby requiring multiple measurements of  ⃗ () at distinct values of .
The number of required frequencies is determined by the number of unknowns in the network or network cluster,  unknowns .Assuming a unique mass and damping for each resonator and setting the edge spring constants to zero, we write  unknowns in terms of the number of resonators in the network cluster: The full set of unknown parameters include   ,   ,   , all relevant   's, and the amplitude .The cluster size ( cluster ) is determined by the number of resonators with detectable oscillation, and will always be equal to or less than the total number of resonators in the network.Moreover, the response vector  ⃗ () will have  cluster non-zero components (although some might be zero or undetectable when sampled at some frequencies) and  ̃() will be a symmetric, tridiagonal  cluster ×  cluster matrix.Each measurement of  ⃗ () at a fixed value of  solves and thus provides two equations per resonator (due to real and imaginary parts): Therefore, to determine the unknowns, we only need to acquire response vectors at a minimum of two frequencies- ⃗ (  ) and  ⃗ (  ).These response vectors then provide a system of 4 ×  cluster linear equations.
In cases where the size of the resonator network is finite, we may have the cluster encompassing the entire network, such that  cluster =  unknowns .In some experimental cases, motion of some resonators at a given  will be small and undetectable by the apparatus and require acquisition of additional response vectors  ⃗ ().When a resonator's motion is undetectable at some or all driving frequencies, then the number of linear equations is reduced; if the undetectable resonator is in the interior of the cluster, we lose six equations, while if at the edge of the cluster we lose four equations.Therefore, in these cases it is important to ensure that all resonators have a measurable amplitude at the two driving frequencies, or else to measure  ⃗ at additional driving frequencies.
Reorganizing the system of linear equations.(S8) We combine and reorganize Equations (S7) into the following single linear homogenous equation where  is a real-valued matrix with matrix elements that depend on the known measured quantities   (),   (), and  (step 4 of Fig 1c).See below for the form of  for the monomer and dimer case.When constructing  in practice and calculating the error of  ⃗ () , the vector  ⃗-which we call the parameters vector-is an  unknowns -dimensional vector composed of the unknown mechanical parameters of the cluster:   ,   ,   , all relevant   's, and the force amplitude vector  ⃗ .Including  ⃗ () at additional driving frequencies will expand the dimensionality of  and  ⃗, which in general will over-determine the system of equations.Finally, we use singular value decomposition (SVD) to find the solution space of , which determines the output parameters vector  ⃗ ̂.
Here, we simulate cases where the force amplitude vector is  ⃗ = 〈 1 , 0〉, only driving the first resonator, but other force amplitudes 〈 1 ,  2 〉 are possible.Sampling response vectors from at least  = 2 drive frequencies, we then use equations of the form Eq. (S7) to construct  and the system of linear equations  ⃗ = 0 ⃗⃗ .We apply a NumPy SVD solver [4,5] in Python to factorize and identify the solution space of , and thereby generate a solution to recover the parameters  ⃗ ̂, where the hat symbol indicates that these are values recovered by the SVD analysis, not simulated measurements (step 5 of Fig 1c).
We obtain from the NumPy SVD solver the singular values  1 ,  2 , … ,   , sorted from smallest to largest, each associated with a singular vector  ⃗ 1 ,  ⃗ 2 , …  ⃗  , where  = dim ( ⃗  ) is the number of parameters.With noise, the null value of  ⃗ = 0 ⃗⃗ is not precisely zero, so we assume the singular vector associated with the smallest singular value  1 corresponds to the solution space.We may also consider multiple small singular values  1 ,  2 , … ,   to be degenerate, even if they are not precisely equal to zero, and, in that case, we must find the solution within a -dimensional solution space, for some positive integer  ≤ .

Scaling the singular vector: the question of solution space dimension
Since NumPy SVD solver provides normalized singular vectors, each of length 1, we must therefore scale the SVD output to identify the solution within the -dimensional solution space and recover the numeric parameters.If we assume the solution space is 1D, we use  = 1 known values to scale the output vector.Here, we will require the force amplitude  to equal the input value  in , and thus we obtain a unique scaled solution  ⃗ ̂1D =  ⃗ 1,un , where the normalized (i.e.unscaled) singular vector  ⃗ 1,un corresponds to the smallest singular value  1 and the scaling coefficient  allows  to equal the input value.Here we use an input parameter, whereas in our experimental work, [2] we estimate a parameter.
The smallest singular value  1 will only equal zero exactly in the limit that the noise  approaches zero, and  2 is larger.Nevertheless,  2 may approximate zero, which would correspond to a 2D solution space.If we interpret the solution space as 2D, then  = 2 known values is sufficient to define a unique solution  ⃗ ̂2D =  ⃗ 1,un +  ⃗ 2,un from within the 2D space, and, while any two choices of parameters are options, here we will fix both  and  1 to equal the input values and allow the other parameters to be scaled accordingly.
If we assume the solution space is 3D, we require  = 3 known values to define a unique solution  ⃗ ̂3D =  ⃗ 1,un +  ⃗ 2,un +  ⃗ 3,un ; for these three known values, here we use ,  1 , and  2 for the dimer or , , and  for the monomer, leaving the damping  as the only unknown parameter in the monomer case.
Ultimately, we obtain three possible solutions for the parameters,  ⃗ ̂1D ,  ⃗ ̂2D , or  ⃗ ̂3D , where the 1D solution is the easiest to obtain, requiring the least additional information.For an experimentalist, the choice of solution space dimension (the nullity ) may be informed by running simulations and by considering the size of the singular values.If the second smallest singular value  2 is low, then the 1D solution space may be insufficient and a 2D solution space may be more accurate, as shown in Fig 4d .In such a case, an experimentalist seeking to calculate the physical parameters using NetMAP has three choices: (1) assume a 2D solution space (or higher dimensional) and find a solution in that space, (2) use the 1D solution space solution, recognizing that the error will be larger, or (3) take additional spectral data to improve the accuracy of the 1D solution.In any of these situations, simulations may inform the choice.
In order to validate the accuracy of NetMAP, we use a priori information from the simulations (step  3e).The percent error for each parameter, the absolute value of the fractional discrepancy: describes the accuracy of the output parameters.Because we take an absolute value,   has a halfnormal distribution.When we repeat the simulated experiment multiple times with different random noise for each repetition, then we obtain a logarithmic average  ̅  over multiple trials of the error.We calculate the average error over parameters as where  is the total number of parameters and  is the number of parameters that have been fixed such that their error is zero, which should not contribute to the arithmetic mean.We consider the mean error ⟨⟩ ̅̅̅̅ averaging over parameters and multiple trials, to be the figure of merit that describes the accuracy of the SVD approach in describing a resonator system, with a lower error ⟨⟩ ̅̅̅̅ indicating a high accuracy.The input parameters are a priori information, only available when conducting simulations, and generally not available when conducting experiments, so we seek to predict the average error 〈̅ 〉 from values available to an experimentalist, including the correlation value and the signal to noise ratio (SNR).We compute the correlation R-value: to quantify the agreement between the input and output spectra (step 7 in Fig 1c).To compute the -value, we solve the equations of motion and calculate the output spectra  ̂() for each resonator  from the recovered parameters  ⃗ ̂, and plot both the noisy input spectra  ⃗ () (colorful datapoints) and the output spectra  ⃗ ̂() (black dashed line) for   ≈ 100 frequencies, many more frequencies than we use for NetMAP in the examples here.Since we may plot either the amplitude and phase or the real and imaginary parts of the complex amplitude   , we have a choice of which plots to use for calculating the R-value.For a monomer, we choose both Re() vs  and Im() vs , calculate  2 for each and use the arithmetic mean as our value for  2 in Fig 2d .For a dimer, we choose all four of Re( 1 ) vs , Re( 2 ) vs , Im( 1 ) vs , Im( 2 ) vs , calculate  2 for each and use the arithmetic mean as our value for  2 in Fig 3g .Since we compare parameters output from the SVD analysis and measured spectra, the R-values are accessible to experimentalists, whereas directly calculating the error ⟨⟩ ̅̅̅̅ required a priori knowledge of the parameters.By providing a prediction for the error ⟨⟩ ̅̅̅̅ , we will enable an experimentalist to estimate the accuracy of the parameters  ⃗ ̂ output from the SVD analysis without needing to know a priori the true parameters  ⃗.

Factorial methods
Table S1 Factorial parameters for monomer investigation, 2  ,  = 6,  = 30  To understand how the percent error of the algebraic response vector characterization varies with factors of our system, we ran full, replicated two-level factorial experiments-e.g. 2  experiments, [6] where  = .The input factors of these experiments include spring constants, damping, mass, the applied force, the standard deviation of the input noise of simulated spectra, and the number of frequencies sampled.We varied these factors over an order of magnitude.For the monomer study,  = 6 parameters, and for the dimer study,  = 8.Factors and levels for monomers and dimers are shown in Table S1 and Table S2, respectively.The response variables for this experiment included the average error for one-, two-, and three-dimensional null spaces.In the case of two-and three-dimensional null spaces, we scaled the singular vectors as described above (Methods section 5).We replicated monomer experiments 29 times, for a total of 30 full simulated experiments.We replicated dimer experiments 9 times, for a total of 10 full simulated experiments.The results of the simulations were analyzed using ANOVA factor screening methods, and the model assumptions were checked and verified for each.We analyzed the logarithm of the error ⟨⟩ ̅̅̅̅ because the raw data did not satisfy the equal variances assumption, and the logarithm of a half-normal distribution better approximates a normal distribution.The logarithm stabilized the variance appropriately.We used  = 0.05 significance cutoff for effects and interactions in the reduced model.

Simplifying assumptions
There are several simplifying assumptions we make for this investigation regarding how measurements are made and what information is available to an experimenter.We assume as an approximation that the resonant system is not over-driven, heated, or subject to temperature fluctuations.We assume that masses, spring stiffnesses, and damping coefficients are constant, regardless of the driving frequency or force.We assume that the experimenter has a correct model for the number of resonant masses and springs and their connections, and that the system has reached a steady-state oscillation, where the energy dissipated by the damping balances the energy input by the force.For this investigation, we consider the simplest two systems (Figure S2): a single damped driven resonant mass on a spring, which we call a monomer, and a resonant dimer system consisting of a pair of two masses, each connected to a wall by a spring and coupled together with a third spring.
For both the monomer and dimer presented here, we assume the experimenter drives one mass with a sinusoidal force  =  cos () and measures the amplitude  and phase  of each resonant mass, and the phase is measured with respect to the phase of the driving frequency.
The monomer system (Figure S2a) is fully described by the following parameters: mass , spring stiffness , damping coefficient , and driving force ().We assume that the driving force amplitude  and frequency  are known.If, however, an experimenter does not know the driving amplitude  in Newtons, knowledge of any one of the other parameters, , , or , is equivalently useful, or the experimenter can use NetMAP to obtain the force-normalized parameters, /, /, and / without requiring any additional information for a 1D solution space.In our related work, [2] we scale the parameters with respect to the spring constant, calculating /, /, and /.Higher dimensional solution spaces require more a priori information about the parameters.
The dimer system (Figure S2b) is fully described by the following parameters: masses  1 and  2 , with damping coefficients  1 and  2 , each attached to a wall by a spring of stiffnesses  1 and  2 and coupled to each other by a spring of stiffness  1 , also called  12 .We assume that the first mass of the dimer is driven and that there is no external driving force acting on the second mass, though our algebraic approach can be modified to address an additional driving force.

Monomers Simulating a monomer
For a monomer (Figure S2a) driven by a sinusoidal force, the equation of motion in complex form is the differential equation where  is the complex position of the mass, ̇ is its complex velocity, ̈ is its complex acceleration, and the real part of the equation is the physical description.We guess a steady-state solution to this differential equation of the form  =   i , where  is the complex amplitude of the position of the oscillating mass, and, through substitution for  and division by the oscillating factor   , obtain the time-independent equation of motion Separating into the real and imaginary parts, we find the equation of motions requires both (−ω 2  + ) Re() −  = 0 and −  Im() = 0 (S15) We assume an experimenter would measure the complex amplitude  = e i of the single oscillator in response to the applied force e i .

Recovering monomer parameters using SVD
We next calculate the unknown parameters , , and  from the known values , , and .Normally, the response (ω) is measured at a series of driving frequencies  and fit with an iterative least-squares approach that optimizes 1 −  2 , where  2 is the coefficient of determination.We describe an alternative non-iterative approach using singular value decomposition.We require at least  = 2 measurements of the complex amplitude, (ω 1 ) and (ω 2 ), and with singular value decomposition (SVD), the method supports any number of measurements above two, as follows.For each measurement (ω  ), we know that both Eq 3a and 3b must hold, and we organize these equations as follows: ⃗ = 0 ⃗⃗ where the measurement matrix  is calculated from the  spectrum measurements.For the monomer, each frequency contributes two rows to , so  = 2 frequencies provide a 4 × 4 square matrix.
We have flexibility in the number of measured frequencies by expanding to a rectangular matrix, where each additional frequency adds an additional two rows to the matrix.We assume an experimenter's measurement of each element of the measurement matrix  is subject to Gaussian noise affecting  but that the noise in frequency ω  is negligible.In order to algebraically solve the differential equation, we apply singular value decomposition (SVD) to the matrix to obtain its singular values and corresponding 4-dimensional singular vectors.For a square normal matrix with a normal eigenbasis, we may also call these the eigenvalues and eigenvectors.The singular value corresponding to zero has a singular vector corresponding to the parameters vector  ⃗.We use the known force amplitude  to normalize the parameters vector, thereby obtaining the physical parameters  ̂,  ̂, and  ̂, where the hat symbol indicates that these are calculated from our model of the measured data.
Details for the monomer case shown in Figure 2 of the main text , and  ̂= 15.99647 N/m.The percent errors for each of these is −0.022%, −0.0040%, and 0.022%, respectively.Each of these is within 0.022% of the correct values for m, b, and k.We also see that the recovered value  We set the input noise  = 0.0005 m.We use up to 25 frequencies for SVD analysis, namely [1.5795569The spectrum with these 25 frequencies identified is shown in Figure S3, left.For  = 25 frequencies, the  matrix has 50 rows.Its smallest singular value,  1 = 0.00059708 , corresponds to singular vector  ⃗ ̂ = ( ̂,  ̂,  ̂, ) =α( −0.36955 kg, −0.036957 N/(m/s), −0.92387 N/m, −0.092383 N), where α= in / −0.092383 N = 1 / −0.092383 = −10.825 is a normalization constant obtained from our knowledge of the force amplitude f for a 1D-SVD analysis.Dividing by α allows us to scale the singular vector to yield the modeled parameters vector.Therefore, we obtain  ̂= 4.00017 kg,  ̂= 0.400043 N/(m/s) and  ̂ = 10.00047N/m.The percent errors for each of these is 0.0041835%, 0.010689%, and 0.0046745%, respectively.Each of these is within 0.010689% of the correct values for m, b, and k.We also see that the recovered value Heatmaps in Figure 4c of the main text show how accurately  = 2, 1D-SVD and 2D-SVD recover the solution for this monomer.Figure S3, right shows a similar heatmap for  = 2, 3D-SVD.For a monomer, 3D-SVD only recovers damping ; all other parameters must be known a priori in order to find the solution in the 3D solution space, so the average error is showing the error in damping ̅  (Figure S3, right).Figure S4 shows how the average error for the monomer parameters decreases for 1D-SVD as the number of frequency points increases.A connecting letters report (Table S3) shows that there are diminishing returns to increasing the number of frequency points:  = 25 is not significantly better than  = 7.  Figure S5.The 1D solution is much less accurate when  2 is small.For  2 near zero, the error of the 1D solution rises to 10%, and for larger  2 , the 1D solution error is as low as 0.01%.This is the same dataset shown in Fig, 4f of the main text.Table S3 Connecting letters report for the data shown in Figure S4, corresponding to Fig. 4b.

Dimers Simulating the dimer
We consider how the algebraic approach presented here applies to a dimer system of damped driven oscillators: two masses  1 and  2 , with damping  1 and  2 , respectively, each coupled to a wall by a spring of stiffness  1 and  2 , respectively, and coupled to each other with a spring of stiffness  12 (Figure S2b).The equations of motion of the damped driven dimer in complex form are the coupled differential equations where  1 is the amplitude of the driving force applied to the first oscillator and  2 is the amplitude of the oscillating driving force applied to the second oscillator.The two forces are assumed to have the same frequency and phase, such that the force amplitudes  1 and  2 are each real constants, or one of the two is assumed to be zero.We guess a steady-state solution of the form  1 =  Recovering the dimer parameters with SVD Having now simulated an experiment, we consider how we can recover the physical parameters from the noisy spectra  1 (ω) and  2 (ω) and knowledge of the driving forces  1 =  1  i and  2 =  2   .We reorganize Eq (B2) into a different matrix equation than we used for Cramer's rule by creating a vector  ⃗ of parameters and a matrix  of measurements such that: where   = Re (  (  )) and   = Im (  (  )) provide information about measurements of each resonator at each of the discrete frequency points.
The measurement matrix  has dim( ⃗)=9 columns and 4 rows for the dimer, where  is the number of frequency points used for the analysis.In general we have  cluster complex equations describing the system.Taking the real and imaginary parts, we rewrite these as 2 cluster real equations and thus the real measurement matrix  has 2 cluster rows.For example, using measurements from two frequencies to describe a dimer would require an 8 × 9 matrix whereas using 30 frequencies would require a 120 × 9 matrix .We then analyze the measurement matrix  using SVD to identify the singular vectors and identify 1 or more singular vectors corresponding to the singular value of zero.The dimension of the solution space may be open to interpretation because the measurements are noisy, so the smallest singular values do not precisely equal zero.If we interpret the solution space to be 1-dimensional then there is one singular vector corresponding to the null singular value.We therefore scale the singular vector using our knowledge of one element of the parameters vector  ⃗.In particular, here we assume we know force amplitude  1 pushing the first mass.Hence the analysis recovers all but one of the elements of the parameters vector  ⃗ = ( 1 ,  2 ,  1 ,  2 ,  1 ,  2 ,  12 ,  1 ,  2 ).The resulting output parameters ̂ are not exactly equal to the input parameters  ,in due to noise in the measurements.If we do not add noise, the recovered parameters are exactly equal to the input parameters, within computer precision.Therefore, we add noise in order to validate NetMAP with noise.

Statistical comparison tests for dimer case
Table S4 shows the statistical comparison tests for the dimer case shown in Figure 3 of the main text.
Table S4.Observed -statistics and associated -values for network parameters of the dimer shown in Fig. 3 of the main text.Error is inversely proportional to the SNR Since the error is proportional to the input noise , the error is also inversely proportional to the SNR, as shown in Figure S6.

Null-space Dimension
Figure S6 For the dimer shown in Fig. 3 of the main text, error is inversely proportional to the SNR.This is another view of the data shown in Fig. 3h of the main text, where input noise is varied.The average SNR for resonator 1 (R1) and resonator 2 (R2) also varies, and the error of the recovered parameters is proportional to the input noise.The lines here are noisy but linear.
Details for the dimer case shown in Figure 5 of the main text Here we provide details relating to the dimer case shown in Figure 5 of the main text.Figure S7 shows the R1 and R2 spectra for this dimer.The second resonance peak in (a), near 3.501 rad/s, is imperceptible on the amplitude spectrum but does show an effect in the phase spectrum.Figure 5c of the main text shows how the accuracy of the 1D solution improves as additional response vectors are incorporated into the analysis.The improvement shows diminishing benefits as the number of response vectors becomes large.A connecting letters report (see Table S5) shows which of the errors are significantly different from each other.Levels sharing a letter are not significantly different from each other.

Experimental considerations for measuring absolute phase
Measuring the absolute phase requires some experimental consideration because the time elapsing between the driving force and the response usually is measured by a lock-in amplifier, and the electronics and optics will introduce delays before the driving force arrives at the resonator system and additional delays before the measurement can be read, creating a phase delay θ =  delay , which is added onto the absolute phase .We assume that this time delay can be measured and subtracted to obtain an absolute phase, as follows [2].An experimenter may obtain the time delay by measuring the spectra at lower frequencies, where the force and response are expected to be in phase such that the absolute phase  is zero.Then the total phase θ +  =  delay + 0, such that a linear fit of the phase versus frequency provides the slope  delay .The phase delay  is subtracted for each frequency  to obtain the absolute phase .

Additional Cases
Heavily damped monomer and thus construct the real rectangular 6 × 4 matrix (step 4) 0.0025 m.We measure at the two resonance frequencies, 0.758 and 4.02 rad/s.Then the mean SNR for resonator 1 measurements is 239 and for resonator 2 measurements is 202.For the strongly coupled dimer shown in Figure S10, we find that the 3D solution is the most accurate, while the 1D and 2D solutions are similar to each other.
Force applied to both resonators of a dimer As an additional example, we consider applying a force  ⃗ = ⟨ cos  ,  cos ⟩ to a dimer.That is, applying the same oscillating force to both masses.We set the input values to:  For dimer system in Figure S11 with force applied to both resonators, we find that the 3D solution is slightly more accurate than the 2D solution, and the 1D solution is the least accurate.The amplitude  2 of the second resonator is higher because it is also driven, and therefore the signal to noise ratio for the second resonator is better than cases where only one of the two resonators is driven.This explains why the solutions are overall more accurate than for dimer cases with only one driven oscillator.

Factorial experiments varying every parameter between two levels
Plots summarizing the factorial experiments are provided here for monomers (Figure S12) and dimers (Figure S13 and Figure S14).

Figure S12
The 1D-SVD error for a monomer is inversely proportional to the SNR, and varies with , , and .Each plot shows the complete set of data from the 2   experiments for the monomer described in Table S1.The color coding in each varies to show how the relationship between average SNR and average error varies with each of the six varied parameters, 2  ,  = 6,  = 30.
An experimentalist using NetMAP will wish to know the accuracy of the parameters vector without knowing the input parameters.A simplified linear model using  2 and mean SNR (Figure S13) can estimate the average error for a wide range of situations.

Figure 1 .
Figure 1.a) A resonator network as a chain of mass and spring oscillators.b)We calculate the relationship between the complex amplitude  ⃗ () of the masses and the underlying parameters.c) Steps of the validation process: 1) Modeling a mass-spring network, such as a monomer, shown here, or dimer.2) Setting the parameters  ⃗  to simulate.3) Solving the equations of motion (EOM) and simulating spectra with noise.4) Using simulated measurements of amplitude and phase at discrete frequency points to construct a matrix.5) Using Singular Value Decomposition (SVD) to recover the parameters  ⃗ ̂, and scaling the result as needed.Here  is shown without a hat to identify it as a known quantity while the others are scaled to . 6) Calculating the percent error  to compare the recovered parameters to the set parameters.7) Calculating the expected spectra for the recovered parameters, and calculating  2 to compare the simulated data and the expected curve.We find that  2 is correlated with error .
As an additional test of NetMAP accuracy, we compare the predicted  ̂() to the noisy simulated (, ) with correlation analysis. ̂() is shown as a black-dashed curve in Fig.2a,b.Using   = 100 simulated spectral data points, we compute the correlation coefficients   2 and   2 , where  = Re() and  = Im().For the data in Fig. 2b, the deviation of the  2 values from unity are 1 −   2 = 8.6 × 10 −8 and 1 −   2 = 8.0 × 10 −8 , or an average of 1 −  2 = 8.3 × 10 −8 .These  2 values indicate that the NetMAP prediction accounts for essentially all the variation of simulated spectra, despite calculating  ̂() from response vectors at just  = 2 frequencies.

Figure 2 .
Figure 2. a) Amplitude () and phase () spectrum of the lightly damped monomer, with input parameters: mass  = 4 kg, damping coefficient  = 0.01 N m/s , spring stiffness  = 16 N/m, and force amplitude  = 1 N.A subtle grey curve shows the exact spectrum while the datapoints in color show simulated spectrum measurements with standard deviation  = 0.005 m.The black dashed line shows the output spectra  ̂() and  ̂().The color scale corresponds to (b).b) The complex amplitude  =   is plotted in the imaginary plane, where the datapoints are phasors,  is the distance from the origin, and phase  is the polar angle.These are the same data as (a).Two measurements [(  ) and (  ), where   = 2.0000 rad/s and   = 2.0013 rad/s, black circles] are selected for input to NetMAP: one at the frequency of maximum amplitude and the other at a frequency corresponding to  = − 3 4

Figure 3 .
Figure 3.For an example dimer system with selected parameters, we demonstrate sampling at  = 2 frequencies, each at a resonance peak (black circles in (a,b,c,d)).a) Amplitude  and phase  spectra showing the motion of resonator 1 (R1).b) Spectra of resonator 2 (R2).c) Complex spectrum of R1, showing the same dataset as (a).d) Complex spectrum of R2, showing the same dataset as (b).e) Box and whisker plot showing the spread of recovered parameters as fractional discrepancy, Δ  / ,in = (̂ −  ,in )/ ,in , for 1D, 2D, and 3D solutions over multiple trials.For 2D-SVD, an additional parameter,  1 , is fixed at  1, in order to identify the solution within the 2D solution space.For 3D-SVD,  2 is also fixed.f) Histogram and box plots of the average error ⟨⟩.g) The average error ⟨⟩ is correlated with 1 −  2 .h) High SNR is key to minimizing the error.The error ⟨⟩ ̅̅̅̅ is proportional to the input noise .

Figure 4 .
Figure 4. Comparing various choices of frequency measurements for analyzing a monomer with  = 4 kg,  = 0.4 N m/s ,  = 10 N/m,  = 1 N, resonance frequency  res = 1.58 rad/s, input noise  = 0.0005 m. a)Input (circled datapoints) and output (dashed black curve) spectrum for a simulated monomer with 25 circled points used for analysis.b) In order to compare the number  of frequency measurements, we plot the average error for m,k,b as a function of the number of frequency points used in the SVD analysis.The lines show the mean error ⟨⟩ ̅̅̅̅ of 100 simulated trials (shown as datapoints) for each number of frequency points.Violin plots show the distribution of ⟨⟩.For this monomer system, the 2D-SVD solution (orange) has the least error, then the 1D-SVD (blue), and the 3D-SVD (green) has the largest error.The average percent error of the parameters  ⃗ varies with the number of frequency points used in the SVD analysis.For this analysis, we add the frequency points in a specific order, starting with two points at each of the two amplitude peaks and sequentially adding points to the left and right of the peaks.The specific datapoints used for b are indicated in a. cdef) Comparing the choice of frequencies if there are  = 2 frequencies.c) The average error of the SVD solution varies with the choice of the two frequencies   and   .The 1D solution fails when the two frequencies are the same.Both the 1D and 2D solution have better results when the measurements are taken near the peak of resonance.d) The average error varies with the measured frequency   .Here the average is taken over all the results shown in (c), with   varying from 1.4 to 1.8 rad/s.e) For a double ( = 2 frequency) measurement of a monomer, the average error ⟨⟩ ̅̅̅̅ varies with the choice of frequencies measured.In this case, we fix   at the resonance peak and sweep   .We find that the optimum second frequency occurs when the phase at   is near   = −

𝜋 4 or 4 .
= − 3 The colored bands correspond to 95% of the solutions.f) The second smallest singular value  2 predicts the average error of the SVD solutions across all combinations of frequencies from trials shown in (c).
coupling spring  12 = 5 N/m,  1 = 0.5 N/(m/s),  2 = 0.1 N/(m/s),  1 = 1 N, and noise  = 5 × 10 −5 m.Oscillating force is applied to  1 .The simulated spectra are shown in Fig. 5a,b (see Supplemental Material for corresponding amplitude and phase spectra); the spectral peaks inferred from the  2 () spectrum are at  res = 0.774 rad/s and  res = 3.501 rad/s.The higher frequency peak of  2 () has a markedly smaller maximum amplitude.The average errors for up to  = 50 sampled response vectors are shown in Fig.

Figure 5 .,
Figure 5. Analysis of optimal frequencies for NetMAP with an example two-mass (dimer) resonator system, with  1,in = 8 kg,  2,in = 1 kg,  1,in = 2 N m ,  12,in = 5 N m ,  2,in = 7 N m ,  1,in = 0.5 N m/s ,  2,in = 0.1 N m/s ,  1 = 1 N, and  = 5 × 10 −5 m. a,b) Complex spectrum of R1 (a) and R2 (b), showing the simulated spectrum with noise (colorful datapoints) and the spectra of the recovered parameters  ̂() (black dashed line, for 1 trial).White plus signs indicate datapoints that are added as  increases from 2 to 50, and numbers indicate the order in which the points are included, with odd numbers for the larger peak and even numbers (not all shown) for the smaller peak illustrating how we alternate adding the frequency points.Only one resonant frequency,  res = 0.774 rad/s (blue points), appears in the spectrum of R1 (a), but a small second peak at  res = 3.501 rad/s (yellow points) appears in the spectrum of R2 (b).c) The average error varies with the number of datapoints  used in the analysis, with 1D-SVD in blue, 2D-SVD in orange, and 3D-SVD in green.For this analysis, we add the frequency points in a specific order, starting with two points at each of the two amplitude peaks and sequentially adding points to the left and right of the peaks.Simulated with 99 trials per each of 49 values of , totalling 4851 trials.d) If  = 2 datapoints are used, the 1D solution is most accurate when a measurement is taken at each of the two resonant peaks.e,f) The accuracy of the 2D and 3D solutions also varies with the pair of frequencies chosen.In general, it is less accurate for the two frequencies to be equal or nearly equal (dark diagonal lines).g) Replotting (d) in terms of R2 phase rather than frequency shows greater detail.Subfigures (defg) are simulated with 240000 trials.h) When   is fixed at the lower resonant peak,   = 0.774 rad/s, then the accuracy of the 1D solution is optimized when   is at the higher resonant peak, 3.533 rad/s (blue dip).The 2D solution is not accurate when each measurement is at one of the two resonant peaks (sharp orange peak).Shaded areas show 95% of simulated ⟨⟩ results and lines show ⟨⟩ ̅̅̅̅ .Simulated with 1600 trials per each of 199 values of   , totalling 318400 trials.

Figure 6 ., 1 N
Figure 6.Varying dimer parameters affects the resonance frequencies and average error.In (a -e), the input parameters are  1 = 1 kg,  2 = 10 kg,  1 = 1 N/m,  12 = 1 N/m,  2 = 10 N/m, damping  1 =  2 = 0.1 N m/s , and force amplitude  1 = 10 N, except the single parameter that is varied: force amplitude in (a), mass 2 in (b), coupling spring in (c), spring  1 in (d), and spring  2 in (e).In (f),  2 is varied and the input parameters are  1 = 5 kg,  2 = 3 kg,  1 = 12 N/m,  12 = 1 N/m,  1 = 1 N m/s ,  2 = 0.5 N m/s , and  1 = 10 N. The resonator cartoons above each plot indicate the strength of parameters with line thickness, and the varied parameter is indicated in orange in the cartoon.The average error of the recovered parameters is plotted as a function of the varied parameter for three solution space dimensions, with 1D-SVD (blue), 2D-SVD (orange), and 3D-SVD (green).Error for individual simulations ⟨⟩ is shown as datapoints and the average error across the simulations ⟨⟩ ̅̅̅̅ is shown as lines.

2 :
: number of frequencies for calculating -value    2 : coefficient of determination comparing the amplitude spectrum  ̂() to the measured amplitude   ().  coefficient of determination comparing  ̂() to   ()    2 : coefficient of determination comparing  ̂() to   ()    2 : coefficient of determination comparing  ̂() to   ()  2 : The average of    2 and    2 , taking the average over all resonators  : monomer spring constant : number of parameters in factorial.Used in the context: 2  .Force  ⃗ : force vector, with components for each oscillating mass.The capital notation indicates that it is time-dependent (oscillating).dim( ⃗ ) =  cluster  1 or : the oscillating force pushing mass 1.  ⃗ , force vector, just the amplitude (with the time-dependent part divided out). ⃗ =  ⃗ cos   1 or f, the force amplitude pushing mass 1.  1 =  1 cos  Results of singular value decomposition  1 ,  2 , … ,   the singular values of , sorted from smallest to largest  ⃗ 1 ,  ⃗ 2 , …  ⃗  , the associated singular vectors  1 ,  2 , … ,   the subset of singular values that correspond to the solution space  ⃗ 1,un the unscaled singular vector associated with the smallest singular value  1 .It is normalized to length 1.  ⃗ 2,un the unscaled singular vector associated with the second smallest singular value  2 .It is normalized to length 1.  ⃗ ̂1D ,  ⃗ ̂2D ,  ⃗ ̂3D : the recovered parameters, calculated using  = 1,  = 2, and  = 3, respectively.Error The error is a figure of merit to show how accurate recovered parameters are. the percent error in mass.The discrepancy, Δ = ̂ −  ,in .For example, Δ =  ̂−  in .̅  : the logarithmic average of the error, taken over many trials, where  indexes the parameter (for box and whisker plots) average across parameters for an individual simulation trial 〈̅ 〉: logarithmic average across trials first, then arithmetic average across parameters ⟨⟩ ̅̅̅̅ : arithmetic average across parameters first, then logarithmic average across trials.

Figure
Figure S1 a) An oscillator network made up of a linear chain of masses and springs.b) Pictorial visualization of equation of motion showing the cluster response vector  ⃗ () and submatrix ℳ ̃().

6 of− 1 ,
Fig 1c) to calculate the discrepancy Δ  = ̂ −  ,in between the input and output parameters, where  indexes the parameters.The fractional discrepancy, has a normal distribution over multiple trials and is shown as box and whisker plots (Fig 2c and Fig.

Figure
Figure S2 Cartoon of a) monomer and b) dimer system with two coupled masses.

2 . 4 N,
2D Solution space for the monomer exampleIn our example, the vector   ⃗ 2,un +   ⃗ 1,un = [-1.62549828e-02,9.98136999e-03, -5.83247471e-02, 9.98115410e-01] +  [-2.42090968e-01, -6.05227330e-04, -9.68363873e-01, -6.05227821e-02] spans the 2D solution space, but we need to determine both  and  to obtain the physical parameters.For example, if we assume the experimenter independently knows both   = 1 N and  in = 4 kg, then it is straightforward to obtain the coefficients  = −16.5227and  = − 6.60069e-7 to then obtain  ̂= 0.009999992 and k̂ = 15.999999993.Generating the curves from these two recovered parameters yields   2 = 1 − 1 × 10 −11 and   2 = 1 − 4 × 10 −7 .In this case, the 1D solution and 2D solution were almost the same.We note that  ≫ , thus selecting a 2D solution almost entirely in the 1D solution space.The 1D solution is preferred because it requires independent knowledge of only 1 parameter, and we find that it is often reasonably accurate, even when the 2D or 3D solution is more accurate.Details for the monomer case shown in Figure4of the main text For the monomer in Figure 4, the maximum number of frequencies is  = 25.We set the input values to:  = 4 kg,  = 0. = 10 N/m,  = 1 N. Then the resonance frequency  res = 1.58 rad/s and rad/s is more accurate than the individually recovered values for mass and spring stiffness; this is often true for various sharply peaked systems.The percent error for √  m ̂ compared to √  in  in is 0.00025%.

Figure
Figure S3.(left) The spectrum also shown in Fig 4a, with corresponding color scale.Vertical lines and circles show the  = 25 frequencies used for analysis. is the phase.(right) Heatmap showing the average error for  = 3,  = 2 analysis of the same monomer system analyzed in Figure 4.
FigureS4shows how the average error for the monomer parameters decreases for 1D-SVD as the number of frequency points increases.A connecting letters report (TableS3) shows that there are diminishing returns to increasing the number of frequency points:  = 25 is not significantly better than  = 7. FigureS5expands on Fig.4f, showing the variation of the error for 1D-SVD for both of the smallest two singular values.Since the error varies more with  2 than with  1 , Figure4fin the main text shows the variation only with  2 .

Figure S4 .
Figure S4.Average error ⟨⟩ for , ,  as a function of the number of frequency points used in the SVD analysis of a monomer.These are the same datapoints as Fig 4b for 1D-SVD in the main text.

Figure S7 .
Figure S7.Spectra corresponding to Fig. 5 of the main text.a) Amplitude  1 (), b) Phase  1 (), c) Amplitude  2 (), d) Phase  2 ().The black circles and vertical grey lines indicate the two input frequencies   = 0.774 rad/s and   = 3.501 rad/s, each at a resonant peak, used for the SVD analysis ( = 2.The subtle grey curves show the exact spectra calculated from input parameters using Cramer's rule.The colorful datapoints show simulated spectra with added noise  = 5 × 10 −5 m.The dashed black line shows the output spectra  ̂ and  ̂ calculated from the 1D-SVD recovered parameters  ⃗ ̂ for one trial, showing agreement with the exact spectra in grey and with the noisy simulated spectra in rainbow colors. As an additional example, we consider a heavily damped monomer (step 1).We simulate a spectrum for a mass-and-spring where we set  = 4 kg,  = 8 N/(m/s),  = 9 N/m,  = 1 N, and Gaussian noise with standard deviation σ = 0.0005 m (step 2).With these input parameters, the quality fac-5 rad/s greatly overestimates the resonant frequency for the heavily damped case.)We mimic an experiment by simulating the spectrum with noise.We can use any frequency points for the analysis with SVD, and to demonstrate this we choose  = 3 frequencies at random, ω 1 =0.183389 rad/s, ω 2 =0.545455 rad/s, and ω 2 = 1.481318 rad/s, as shown in FigureS8.At these three measurement frequencies, we obtain amplitude (ω 1 ) =0.111738 m, (ω 2 ) = 0.112556 m, (ω 3 ) =0.083927 m and phase ϕ(ω 1 ) =-0.155537 rad, ϕ(ω 2 ) = -0.507897rad, and ϕ(ω 3 ) =-1.542800 rad (step 3).The signal to noise ratio for each measurement is and the mean SNR is 205.We use a parameters vector  ⃗ = (, , , ) 1 = 1 kg,  1 = 0.1 N s/m,  1 = 1 N/m,  1 =  2 = 10 N,  2 = 10 kg,  2 = 0.1 N s/m,  2 = 10 N/m,  12 = 1 N/m,  = 0.005 m, and we measure at the resonance frequencies, 1.00 and 1.45 rad/s.Then the mean SNR for resonator 1 measurements is 16130 and for resonator 2 measurements is 10630.The results are shown in Figure S11.

Figure
Figure S11 An example of a dimer where the same oscillating force is applied to both resonators

Figure S13 .
Figure S13.The results of the factorial simulations for dimers are shown on the left, showing that the error ⟨⟩ (vertical axis) varies with 1 −  2 (horizontal axis) and SNR (color scale).On the right, a linear fit to the simulated data offers a simplified model for estimating the error from the SNR and  2 values.

Figure
Figure S14 Reduced model effects summary of factorial experiment on the dimer system.The stronger effects and interactions have larger LogWorth values.Plots were generated in JMP software v 16.1.

Table 1 :
Observed -statistics and associated -values for network parameters for the monomer shown in Fig.2.Results for the 1D, 2D, and 3D solution spaces are shown.
A useful measure of the accuracy of NetMAP is the fractional discrepancy Δ  /( in )  = (̂ − ( in )  )/( in )  .The Δ  /( in )  distributions for the 1D and 2D null-spaces are shown in Fig.
Table S2 Factorial parameters for dimer investigation, 2  ,  = 8,  = 10 %.This high accuracy likely arises because we choose frequency   at the peak amplitude and the lightly damped monomer has a sharply peaked resonance, so √/ is well defined.