Integration of spectator qubits into quantum computer architectures for hardware tuneup and calibration
IIntegration of spectator qubits into quantum computer architectures for hardwaretuneup and calibration
Riddhi Swaroop Gupta, ∗ Luke C. G. Govia, and Michael J. Biercuk ARC Centre of Excellence for Engineered Quantum Systems,School of Physics, The University of Sydney, New South Wales 2006, Australia Raytheon BBN Technologies, 10 Moulton St., Cambridge, MA 02138, USA
Performing efficient quantum computer tuneup and calibration is essential for growth in systemcomplexity. In this work we explore the link between facilitating such capabilities and the underlyingarchitecture of the physical hardware. We focus on the specific challenge of measuring (“mapping”)spatially inhomogeneous quasi-static calibration errors using spectator qubits dedicated to the taskof sensing and calibration. We introduce a novel architectural concept for such spectator qubits:arranging them spatially according to prescriptions from optimal 2D approximation theory. We showthat this insight allows for efficient reconstruction of inhomogeneities in qubit calibration, focusingon the specific example of frequency errors which may arise from fabrication variances or ambientmagnetic fields. Our results demonstrate that optimal interpolation techniques display near optimalerror-scaling in cases where the measured characteristic (here the qubit frequency) varies smoothly,and we probe the limits of these benefits as a function of measurement uncertainty. For morecomplex spatial variations, we demonstrate that the NMQA formalism for adaptive measurementand noise filtering outperforms optimal interpolation techniques in isolation, and crucially, can becombined with insights from optimal interpolation theory to produce a general purpose protocol.
I. INTRODUCTION
The scale-up of near-term quantum devices presentsnew challenges associated with fabrication, tuneup, char-acterization, and control as the total number of qubitsincreases [1–5]. For instance, in larger devices, we si-multaneously see that fabrication tolerances may lead toperformance variation across systems [6], while the num-ber of measurements required in order to bring systemsonline grows rapidly. These concerns add to the generalissue of decoherence due to undesired coupling of qubitsto spatially inhomogeneous ambient fields. Accordingly,the procedure by which quantum devices must be char-acterized and tuned-up presents a challenge which growswith the size and complexity of each new generation ofdevice [7, 8].In response, there has been an emergence of interest inthe deployment of automated and nondestructive datainference and control approaches to this problem, ex-ploiting the presence of spatial correlations in the per-formance variations [8]. One interesting approach to thisproblem involves characterization of classical noise or de-vice variations in space using spectator qubits [9]; addi-tional dedicated devices are employed to gain informa-tion about hardware performance such that information-carrying data-qubits may operate undisturbed. Evenwithin this framework we are still faced with the chal-lenge of efficiently extracting information from poten-tially large numbers of these ancillary sensing devices.In such circumstances, classical adaptive filtering tech-niques can play an essential role in efficiently charac-terising, calibrating and controlling mesoscale quantum ∗ [email protected] devices in an attempt to extract improved performance[2, 4, 10]. In earlier work, the authors presented twodistinct methods with radically different applications forthe approximate reconstruction of unknown continuousphysical phenomena using discretized projective mea-surements on qubits. The first method [11] was an ef-ficient approach to quantum state tomography of single-mode continuous variable systems using bivariate opti-mal Lagrange interpolation at the Padua points [12–16].The second method was inspired by probabilistic roboticsand enabled adaptive scheduling of measurements on 2Dmulti-qubit arrays for noise spatial mapping, referred toas Noise-Mapping for Quantum Architectures (NMQA)[17, 18].In this work, we combine these apparently disparatetechniques in order to address the problem of how tomost efficiently extract useful information for mesoscaletuneup and calibration using spectator qubits. We modelthe underlying hardware performance variations as a“field” and thus seek the most efficient means to char-acterize the spatial field. To this end we are inspiredby the observation that ideal interpolation theory pro-vides insights on how to conduct optimal sampling onany continuously-varying, bivariate function (e.g. space,time, relative phase). We associate such optimal sam-pling points with the physical locations of sensor-qubitsin order to efficiently determine spatial variations in thefield at the location of the (unmeasured) data-qubits.In particular, we suggest sensor-qubits should be lo-cated on a grid defined by the so-called Padua points[13–16], which enables efficient Lagrange polynomial in-terpolation of the underlying field. We employ numericsimulations of field characterization using different algo-rithmic approaches to estimation and explore the influ-ence of the underlying sensor-qubit locations on the qual-ity of the estimates. We demonstrate that Lagrange- a r X i v : . [ qu a n t - ph ] A p r Padua reconstructions in 2D outperforms other interpo-lation or statistical estimation methods for the recon-struction of polynomial fields, while using ∼ × fewersensor-qubits. Our simulations show how these benefitsare lost in conditions when optimal interpolation theorycannot be applied, and demonstrate that the alternativeadaptive measurement approach NMQA [17, 18] achievesthe lowest expected error. We study the impact of quan-tum projection noise on field estimation quality, includ-ing the role of averaging over the discretized outcomes ofprojective measurements. Our results demonstrate thatNMQA is agnostic to the spatial arrangement of sensor-qubits on hardware, motivating the use of a Padua sensorgrid as a default architectural choice for embedded sen-sors.We outline our approach in Section II, where wepresent a classical spatial reconstruction and predictionproblem over a 2D arrangement of sensor- and data-qubits. Two distinct types of device configurations areconsidered where a sensor-qubit sub-lattice is nested in afixed data-qubit grid. For each device configuration, weintroduce the appropriate spatial reconstruction strategyand test performance. Results from numeric simulationsand comparative analyses are presented in Section III.We conclude with a summary and brief future outlook inSection IV. II. APPROACHES TO SPATIAL FIELDESTIMATION VIA SPECTATOR QUBITS
We consider a notional 2D multi-qubit array in whichqubits are partitioned into ‘sensor-’ and ‘data-’ qubits;these non-interacting qubits are subject to spatial inho-mogeneities in a device parameter ( e.g. qubit frequency),modelled as an external classical field with finite spatialcorrelations. For concreteness we implement a term cou-pling to the qubit Hamiltonian ∝ σ z , as would be the casefor qubit frequency variations or an inhomogeneous mag-netic field. To probe this Hamiltonian term we assumea single-shot Ramsey experiment with a fixed interroga-tion time can be performed on each qubit, though thegeneral framework is compatible with arbitrary measure-ment routines. The question arising in this work is theextent to which the spatial configuration of the sensor-qubits affects the accuracy of the field reconstruction atthe proximal data-qubits.In the limit of high-fidelity measurements on eachsensor-qubit, the task of mapping the field at the data-qubit locations transitions from an estimation problemto a bivariate interpolation or an approximation prob-lem. Let x (similarly, x (cid:48) ) correspond to the spatial coor-dinates of a sensor-qubit (data-qubit) in 2D, and let f ( x )(similarly, f ( x (cid:48) )) be the estimated field value measuredby repeated projective measurements. A key insight inthis work is the recognition that some interpolation tech-niques are associated with an optimal choice of point-setat which to sample, χ := { x i } | χ | i =1 , and that these points should be associated with the spatial locations of sensor-qubits. The canonical example in 1D is the Chebshevpoints. In 2D, the Padua-points are a recently discov-ered family of point-sets for optimally sampling a bivari-ate continuous function, for which Lagrange polynomialinterpolation yields the lowest-error reconstruction [13–16].The Padua points are families of point-sets that occuron the unit-square [ − , . The geometric interpretationof Padua points can be obtained by envisioning them asequally spaced points along a generating curve on theunit square, defined by γ κ ( t ) := ( − cos(( κ + 1) t ) , − cos( κt )) . (1)The intersection of this generating curve with itself, theedges of the square or its vertices yields the Padua pointset. The Padua points are equispaced along t and indexedby j, j (cid:48) as: t ( j,j (cid:48) ) := jκ + j (cid:48) ( κ + 1) κ ( κ + 1) π, j, j (cid:48) ≥ , j + j (cid:48) ≤ κ (2)Each Padua point-set is associated with a positive integer order , κ , which indicates that a Lagrange interpolant canbe constructed at the Padua points using a family ofpolynomials with degree n ≤ κ . Different κ ’s yield Paduapoint-sets, χ κ , that differ in both size and arrangementon the unit square, with the number of points in the setgiven by: | χ κ | = ( κ + 2)( κ + 1)2 . (3)The structure of the Padua grid is unique with respectto the data-qubit grid for κ >
2, while for κ = 1 ,
2, thePadua grid is degenerate with a square grid containingthe relevant number of sites.For any point in the unit square x (cid:48) ∈ [ − , , thePadua-Lagrange interpolant [13, 14] of a true field, de-noted f , is written: L ( f )( x (cid:48) ) := (cid:88) x ∈ χ κ f ( x ) l ( x, x (cid:48) ) . (4)Here the Lagrange basis polynomials, l ( x, x (cid:48) ), form an or-thonormal basis over the Padua points, such that for anytwo Padua points x , x ∈ χ κ , l ( x , x ) = 1 if and onlyif x = x , and l ( x , x ) = 0 otherwise [14, 15]. If eitherone or both of x , x is not a Padua point, l ( x , x ) ≥ x and x using matrices involving the so-called‘product-grids’ of the Chebychev polynomials in 1D, asdetailed in [5], and re-stated for completeness in the Ap-pendices .Locating sensor-qubits at Padua points may be con-trasted with methods in which the sensor-qubit array ismatched to the underlying geometry of the data-qubits,but incorporated as a nested, regular sub-lattice. In thiscircumstance, a polynomial basis need not be an optimalchoice for interpolation. We may instead perform inter-polation using radial basis functions (RBF) [19, 20]. TheRBF method uses a non-polynomial (e.g. sigmoidal orGaussian) function basis to conduct interpolation, andscales well to multi-variate interpolation problems.So far we have discussed deterministic analysis meth-ods such as interpolation, but we also consider nonlin-ear statistical estimation techniques that incorporate aprojective-measurement model. We choose an algorithmpresented in previous literature as NMQA [17, 18] whichestimates an unknown 2D field by adaptively schedul-ing measurements on sensor-qubits. NMQA can be im-plemented directly on any sensor-qubit arrangement, in-cluding both regular (rectangular) lattices or Padua loca-tions for sensor-qubits. However, not all of the availablesensor-qubits need to be measured, as NMQA’s controllerselects the next physical measurement adaptively basedon state estimation in the previous iteration step. It isthus distinct from interpolation methods as it focuses onbuilding a “map” of f iteratively as new measurementsare obtained using a spatial information sharing and pre-diction mechanism. NMQA procedures are characterizedin detail in Refs. [17, 18]. III. RESULTS AND DISCUSSION
In this section, we employ a numeric simulation tocharacterize the performance of these techniques for truefields exhibiting different forms of spatial variation, whileincorporating different spatial arrangement of sensors onhardware. Our studies consider a fixed 5 × f at the locations of the data-qubits onthe array. For any mapping algorithm we calculate errorvia the uniform error, i.e. the infinity norm, and refer tothis quantity as the “mapping error”. The novel data in-ference procedures discussed here are compared againsta baseline approach using a square grid of sensors andsimply assigning sensor results to nearest-neighbor data-qubits (“Standard” protocol). Full details of all hardwaregeometries used in comparative simulations are providedin the Appendices . A. Dependence of interpolation performance onsensor grid
We begin by examining the dependence of mapping er-ror on the choice of interpolation strategy and sensor gridin Fig. 1(b). Calculations take as input data the binaryoutcomes ‘0’ or ‘1’ of simulated projective single-qubitmeasurements and receive the same fixed total measure-ment budget. Measurements are ‘batch-processed’, andthe total measurement budget is fixed as a multiple of the
Qubits ( M app i ng e rr o r (r ad ) - order true polynomial field (arb. units)(a) Regular Padua PaduaMax regularRegular Padua (b)
FIG. 1. Impact of sensor-qubit layout on field mappingquality in a 2D device. (a) Fixed 5 × d × d grid for d = 4(left) or according to Padua points of order κ = 1 , n for three different interpolation strategies. Toprow of insets depicts an example true polynomial field of de-gree n on a fixed 5 × n with random coeffi-cients, and the output of the polynomial is scaled and shiftedto lie between [0 , π ] radians. Interpolation at Padua points oforder κ is performed using Lagrange polynomials with κ = n (black dashed triangles). Interpolation on a d × d regular gridis performed using RBF where the total number of sensorsis within ± d = 2 , , , , d = 9, sensor-qubits thatoverlap with data-qubits are removed, leading to a total of54 sensor-qubits. In the extreme case where the number ofsensor-qubits on regular grid are much greater than Paduagrid for all orders n <
9, RBF interpolation is performed on afixed, regular grid using d = 9 (green dashed markers). Error-bars represent 1 standard deviation within the distribution of50 random trial polynomials. For all interpolation schemes m = 50 measurements per sensor-qubit is used. grid size, so that every sensor-qubit is measured m times.Thus for m >
1, measurements at every sensor-qubit areused to produce an estimate for the expectation value ofthe measurement, and this floating-point number is thengiven to the algorithm.The main panel in Fig. 1(b) shows simulated map-ping error vs. the polynomial order n of different fields f ,as achieved using different interpolation methods. Eachdata point in Fig. 1(b) shows the average mapping errorand standard deviation for a set of 50 randomly sam-pled polynomial fields of a given order n . Examples oftrue fields for different n are depicted as a color-scale inthe top row of insets. The mapping error for Lagrange-Padua interpolation is plotted for κ = n and m = 50measurements per sensor-qubit (black dashed triangles).This dataset may be compared directly with a RBF inter-polation using the maximum number of sensors possibleirrespective of the order n of the field f : a 54 sensor-qubit regular grid (green dashed) with an approximately2 : 1 ratio of sensor- to data- qubits.Since the number of sensors in Lagrange-Padua inter-polation scales with n , we may also directly compare in-terpolation on a regular vs. Padua grid with approxi-mately the same number of sensor-qubits. The choices of d × d regular grid for d = 2 , , , , κ = 1 , , , , n < n ≥
5) correspondto a maximum size difference of ± ±
3) sensor-qubitsbetween regular vs. Padua grids.In Fig. 1(b) we observe that for all polynomial orders n ,the Lagrange-Padua interpolation approach outperformsboth alternative RBF interpolation strategies, demon-strating the impact of the sensor-qubit architecture oninterpolation quality. This observation holds irrespectiveof whether regular and Padua grids are approximatelythe same size or if RBF is provided a distinct advantagewith a high number of sensor-qubits.The inclusion of error in sensor-qubit measurementsdue to quantum projection noise represents a deviationfrom the core assumptions of ideal interpolation theory(for both RBF and Padua), where the value of the fieldat the sensors is assumed to be known exactly. Indeed,comparisons to optimal interpolation theory are only for-mally facilitated in the regime where ideal measurementsare obtained as m → ∞ . Nonetheless, Lagrange-Paduareconstructions give low-error reconstructions of polyno-mial fields for κ = n despite the fact that m is finite,and that all interpolation methods receive noisy mea-surements at the location of sensor-qubits. B. Comparison between interpolation and adaptivemeasurement strategies
Having confirmed the utility of the Lagrange-Paduaapproach, we now relax assumptions of ideal interpola-tion theory in Fig. 2. We plot average mapping error vs. the total number of sensor-qubits, but we change thetrue field so that a representation on the finite polynomialbasis increases in difficulty from (a) to (c), which havefixed underlying fields. For each data point, we simu-late 50 repetitions of the measurement outcomes for thesame underlying field, and report an average mappingerror and standard error.Panel (a) in Fig. 2 corresponds to a simple linear vari-ation, and in panel (b) we plot a common test-functionin interpolation literature known as the Franke func-tion [21]. This function consists of a weak linear back-ground favouring a Lagrange-Padua approach, but withthe addition of a superposition of two Gaussian func-tions, favouring the Gaussian basis in RBF interpolation.In panel (c), the true field is non-polynomial with a form ∝ cos(exp(2 x + y )) sin( y ) that does not permit any finite-polynomial or Gaussian representation.The deviation from fields exactly describable by poly-nomial functions motivates the inclusion of the alternateadaptive mapping strategy, NMQA. In NMQA, a single-qubit binary measurement is fed iteratively into the al-gorithm; NMQA revises the state estimation procedureafter each measurement, and the NMQA controller iter-ation adaptively chooses which sensor-qubit to measurenext. This adaptive measurement procedure means thatfor any fixed total measurement budget, sensor-qubitsunder NMQA are not all uniformly measured m times.We compare the different interpolation and adaptivemeasurement strategies across these three different fieldsusing simulations setting m = 50, but varying thenumber of sensor-qubits employed. Fig. 2(a) numeri-cally demonstrates that the Lagrange-Padua interpola-tion method for low order κ performs the best of ourchosen methods. This benefit arises because the under-lying field is a degree-1 polynomial. In this case, thesimple structure of the underlying field makes the useof additional sensor-qubits with κ > m, κ regime.Moving to other functional forms for f , our resultsfor the Franke function in Fig. 2(b) show that the per-formance of the Lagrange-Padua method improves withincreasing κ , but does not out-perform other methods,particularly for small or moderate values of κ , where itis not even well distinguished from the Standard assign-ment approach. RBF’s superior performance in the high-data limit is somewhat unsurprising since Gaussian com- Interpolation (Padua)Interpolation (Regular)Linear field Franke field Non-polynomial field Q ub i t s ( ) Q ub i t s ( ) Q ub i t s ( ) M app i ng e rr o r (r ad ) True Field (rad)
Total sensors qubits (num)(a) (b) (c)
Standard
FIG. 2. Performance of interpolation strategies and NMQA for reconstructing different true fields for both regular and Paduageometries. True fields are depicted as a colorscale on a fixed 5 × d × d grid, we plotNMQA (filled circles) vs. radial basis functions (blue squares). Standard assignment approach is shown as a threshold usinga fixed arrangement of 9 sensor-qubits and averaged over 4 non-unique orientations (grey solid line). Data shown for m = 50measurements per sensor-qubit, with d = 2 , , , , κ = 1 , , , , , ponents of the Franke function are much stronger thanthe background linear drift, and these Gaussian spatialvariations are naturally expressed in the functional basisof RBF. In this case, NMQA performs marginally betterin the low-data regime but is outperformed by RBF inthe large-measurement limit.Finally, we observe strikingly different behavior when f is not approximated well by the natural functional basisof any interpolant. The final panel Fig. 2(c) represents achallenging field reconstruction problem for all of our cho-sen methods. Interpolation using Lagrange polynomialson Padua grids vs. RBF on regular grids shows that bothinterpolation methods perform similarly and provide lit-tle if any benefit relative to the Standard assignment ap-proach. In this case we observe that NMQA substantiallyoutperforms any other reconstruction strategy, illustrat-ing the relevance of an adaptive measurement strategywhen one cannot use an interpolant in some ideal func-tional basis relative to the true field.A surprising feature of panels Fig. 2(a)-(c) is that theerror scaling behaviour of NMQA appears to dependonly on the total number of sensor-qubits on the unitsquare, and not on the spatial arrangement of sensor-qubits (Padua vs. regular grid), unlike interpolation. Thedifference between NMQA and interpolation strategiescan also examined by varying the number of single-qubitmeasurements per sensor-qubit, m . In Fig. 3, we plotmapping error vs. m for two different field configura- tions (rows) and three different orders for the sensor grid(columns). The limit m → ∞ implies that a true field isknown perfectly on the locations of the sensor-qubits forideal interpolation, indicating that errors in functionalevaluation are reduced from left to right along the x -axisin each panel of Fig. 3.In Fig. 3(a), κ = n = 1 is close to satisfying con-ditions for optimal interpolation on a linear field andLagrange-Padua method provides the lowest error recon-struction for all m . The opposite is true in (d)-(f), whereNMQA significantly outperforms Lagrange-Padua on anon-polynomial field. For the panels with κ > m (as projection noise decreases with m ).However, this behaviour disappears in the more challeng-ing field reconstruction settings of (e) and (f), whereNMQA’s benefits are more substantial. A comparisonof open circles in the top vs. bottom row of panels showsthat mapping error trajectories for NMQA do not changewith the type of true field. This observation empiricallyconfirms that, unlike interpolation, NMQA’s efficacy de-pends on branching random processes [17], and hence itsstructure precludes choosing a functional basis for mapreconstruction.If a priori information about the true field is available, M app i ng e rr o r (r ad ) Measurements per sensor, (msmts/qubit)
Interpolation (Padua) NMQA (Padua) L i nea r N on - po l y no m i a l (c)(f)(e)(d)(a) (b) FIG. 3. Expected mapping error in radians (rad) vs. mea-surements per sensor-qubit m for linear field ( n = 1) in toprow (a)-(c) and non-polynomial field in bottom row (d)-(f).Columns depict increasing order of the Padua grid from left toright κ = 1 , ,
10. In each panel, mapping error vs. m numberof measurements per sensor-qubit is plotted for NMQA (openred circles) vs. Lagrange polynomial interpolation (blackdashed triangles) over 50 trials with visible error bars de-picting 1 std. error. Padua sensor arrangements (red dots)and associated generating curve (red dashed) for each columnis shown as bottom left insets in (d)-(f). For a linear field in(a), κ = n = 1 is optimal Lagrange-Padua interpolation for m → ∞ and (b)-(c) represent an unnecessary increase in thenumber of sensor-qubits. In (d)-(f), a non-polynomial truefield permits no finite order polynomial representation. the results of Figs. 2 and 3 demonstrate that the ap-propriate configuration of Lagrange-Padua methods willyield the lowest error field reconstruction on unmeasureddata-qubits. If no such information is available, thenNMQA satisfies our physical intuition that reconstruc-tion performance improves with an increase in the num-ber of sensor-qubits, whereas this intuition is not alwaystrue for RBF or Lagrange interpolation if the interpolantis poorly specified with respect to the properties of thetrue field.One anticipates that if κ is sufficiently high relativeto the polynomial order of the true field, and func-tional evaluation errors approach zero for large m , thenLagrange-Padua will always reflect a low-error recon-struction relative to RBF and NMQA approaches. Sincethere is no finite polynomial representation of the test-function in Fig. 2(c), we have not been able to numer-ically ascertain the value of κ that results in Lagrange-Padua interpolation outperforming NMQA. IV. CONCLUSION
Under the spectator qubit paradigm, we investi-gated protocols to efficiently characterize spatial inho-mogeneities in qubit calibration or performance, mod-elled as a classical scalar field in 2D. By collecting mea-surements on a dedicated sub-lattice of sensor-qubits, weestimated the field values at the locations of proximal,unmeasured data-qubits in a 2D multi-qubit device, as-suming a dephasing-type Hamiltonian. Drawing from op-timal interpolation and statistical estimation theory, weused simulated data to compare the performance of re-construction methods in the presence of different fieldconfigurations, studying the impact of the underlyingsensor-qubit architecture of field characterization.Our results showed that in most circumstances it is ad-vantageous to arrange sensor-qubits at the Padua points,an optimal point-set for 2D interpolation. In circum-stances where the field is well approximated by poly-nomial functions, Lagrange-Padua interpolation outper-forms comparable interpolation strategies using regularsensor-qubit lattices, as well as adaptive measurementstrategies. If, however, fields are not well approximatedby polynomial functions, adaptive inference proceduressuch as NMQA perform best. Crucially, as we haveshown here, the performance of NMQA shows limitedsensitivity to the underlying sensor arrangement.For complicated spatial variations, our results indicatethat both finite-measurement effects and the presence ofnon-polynomial field structure reduce the utility of op-timal Lagrange-Padua or RBF interpolation strategiesin a spectator-qubit paradigm. In some cases a pri-ori information about spatial variation in a target fieldparameter may be available, but unless the interpolantcan be appropriately specified, increasing the number ofsensor-qubits can increase overall field-reconstruction er-ror. This observation runs counter to the physical intu-ition that providing more sensors in a field reconstructionproblem should improve the quality of the reconstruc-tion. By contrast, the noise-filtering properties of adap-tive strategies such as NMQA enable the algorithm toreduce mapping error as the total amount of informationincreases ( e.g. increasing sensor counts).In this application NMQA also natively accommodatesthe potential to incorporate temporal drifts which are notcompatible with interpolation strategies. The dynam-ical model in NMQA can be modified to include tem-poral dynamics that may either be specified a priori orlearned through data. In the adaptive control [22] andprobabilistic robotics literature [23], particularly in thecontext of spatio-temporal mapping applications for self-navigating vehicles [24], a time-varying environment is of-ten sampled rapidly such that a formal dynamical modelneed not be specified and estimated state information is‘forgotten’ by the algorithm over some time-scale. Con-temporary inference proposals use hidden Markov auto-regressive moving average models [25] or Gaussian pro-cesses [26] to introduce temporal correlations to an oth-erwise static pattern estimation problem [27–29]. A fullspatio-temporal analysis of NMQA remains an excitingsubject for future work.A hybrid of both NMQA and Lagrange-Padua tech-niques yield a general protocol for spatial field character-isation. For higher accuracy, one may first use NMQAfor a coarse grained characterisation of a target field, fol-lowed by a Lagrange-Padua method on local regions thatare approximately polynomial (with interpolation orderand number of measurements informed by the results ofNMQA). For quasi-static linear or quadratic spatial vari-ation, only first or second order polynomial interpolantsare needed, in which case the Padua point-set is an exactsubset of a regular square grid, as explained in
Appen-dices . Thus, NMQA characterisation and local refine-ments via Lagrange-Padua interpolation are practicallyviable even for regular grids. Alternatively, building 2Dqubit arrays in Padua arrangements is equally attractive,as Padua-grid-configured hardware will not impede anyNMQA-based characterisation procedures. We look for-ward to greater exploration of how physical-layer quan-tum computer architectures may be impacted by the con-trol and characterization strategies in use.
ACKNOWLEDGMENTS
This work partially supported by the ARC Cen-tre of Excellence for Engineered Quantum SystemsCE170100009, the US Army Research Office underContract W911NF-12-R-0012, and a private grant fromH. & A. Harley.
CODE AND DATA AVAILABILITY
Access to the code-base and data required to repro-duce all figures is provided via http://github.com/qcl-sydney/nmqa.
AUTHOR CONTRIBUTIONS
R.S. Gupta performed all analysis, led technical ef-forts supporting the results presented, and co-wrote themanuscript. L.C.G. Govia proposed the Padua points forspatial mapping, provided technical guidance on interpo-lation and error analysis, and contributed to writing themanuscript. M. J. Biercuk set the original research di-rection, provided guidance on analysis, and co-wrote themanuscript.
Appendix A: Grid geometries for numericalsimulation
In Fig. 4 we illustrate the grid arrangements that havebeen used throughout the analysis of the main text.The top row of Fig. 4 shows four possible orienta-tions for the “Standard” local neighbourhood assignmentmethod that have been averaged over in all the data-figures presented. Local neighborhood assignments can-not generally be defined uniquely in 2D. For a regulardata-qubit grid defined on a unit square, we choose thefirst four nearest-neighbour data-qubits for every sensor-qubit, thereby maximizing the number of first nearestneighbours in regular lattice. This procedure can con-tinue at the interior points of any regular data-qubit gridwhen the number of rows and columns are even. For anodd row or column, the first nearest neighbours are nec-essarily equivalent to 1D case where two data-qubits areallocated to each sensor. The intersection of odd rowsand columns is a vertex about which no neighbours aredefined. We consider this asymmetric case for a 5 × κ = 1 , κ ≤ Appendix B: Ideal interpolation and approximationtheory
The details for the formalism of ideal interpolation andapproximation theory is the subject of this section. Werestate the theoretical formalism following [30] and es-tablish notation to present the relevant results of [13–16, 31, 32].
1. Least-squares criterion
Let C (Ω) be the space of continuous functions of twovariables defined on an open, connected and boundedsubset Ω ⊂ R [30]. The geometry of Ω affects theinterpolation problem and we focus on the case whereΩ := [ − , takes the form of a unit square. Our truefunction is denoted by f ∈ C (Ω) and f is the subject ofour approximation reconstruction over the region Ω.We now wish to define two members in the space ofbi-variate polynomials. Let P κ be the space of bi-variatepolynomials of degree at most κ with dimension N . Letthe points χ := { x i } Li =1 denote the position of L knownvalues of f , and define a set of associated positive weights as W := { w i } Li =1 ⊂ R + \ { } . We define the best polyno-mial approximation to f as p ∗ , and the minimum leastsquares polynomial as p L satisfying:min p L (cid:107) f − p L (cid:107) = min p L (cid:118)(cid:117)(cid:117)(cid:116) L (cid:88) i =1 w i | f ( x i ) − p L ( x i ) | (B1)This minimum least-squares operator is often re-statedin approximation problems as the interpolant L thatmaps f to p L , L : C (Ω) → P κ . The interpolant, L , isthus a linear operator that depends on the point-set χ ,the weight-set W and the approximation space P κ pa-rameterized by the value κ .The least squares approximating polynomial p L can befound by writing L in an appropriate polynomial basis,and solving for p L . We recap the specific inverse problemand the error properties of L in terms of χ and κ below.Let { p j } Nj =1 be a polynomial basis for P κ and V χ :=[ p j ( x i )] be the L × N Vandermonde matrix for this basisusing the point-set χ . The interpolant at some arbitrarylocation x ∈ Ω is written in terms of coefficients in thisbasis: L ( f )( x ) := ρ ( x ) T c = N (cid:88) j =1 p j ( x ) c j , x ∈ Ω (B2) c := (cid:2) c c . . . c N (cid:3) T (B3) ρ := (cid:2) p p . . . p N (cid:3) T (B4)Here, c represents N scalar coefficients; ρ represents an N -polynomial basis with each element corresponding toa basis-polynomial p j . The coefficients c are solutions toa linear system of equations W V χ c = W F χ , where: W := diag( w . . . w L ) (B5) F χ := (cid:2) f ( x ) f ( x ) . . . f ( x L ) (cid:3) , { x i } Li =1 ∈ χ (B6)The coefficients are thus determined via optimisation orusing the pseudo inverse c = ( W V χ ) † W F χ , yielding: L ( f )( x ) := ρ ( x ) T ( W V χ ) † W F χ , x ∈ Ω (B7)Here, the matrix W is invertible by construction. Thepseudo-inverse governs error properties of the interpolant L .Following [30], the statement that the solution c isunique is equivalent to the condition that V χ has fullrank. If additionally L = N , then the point set χ iscalled uni-solvent for which V χ is invertible and the in-terpolant simplifies to: L ( f )( x ) := ρ ( x ) T V − χ F χ , x ∈ Ω (B8)
2. Conditioning of the interpolant
Using notation in the previous section, we now exam-ine least squares solutions and add perturbations to the R egu l a r P adua S t anda r d Sensor qubitData qubit
FIG. 4. Placement of sensor and data-qubits on the unit square [ − , for Standard local neighborhood value-assignment(top), regular (middle) and Padua (bottom) configurations. A 5 × d × d sensors, with d = 2 , , ,
6. For d = 9, sensor-qubits that overlap exactly with data-qubits are removed, leading to a total of 54 sensor-qubits.Bottom: Padua locations for order κ = 1 , , , , ,
10, and associated generating curve (red dashed). linear system of the previous section. In practical appli-cation of interpolation theory, these perturbations arisefrom inaccuracies in our choice of basis V χ or the point-set χ that lead to inaccuracies in the solution c . Forthe purposes of this section, we present the concept ofa condition number η ( V ) for a choice of a basis for thematrix V . The condition number is unity for an optimalbasis, which means that any errors in the matrix V arenot magnified as errors in the solution c . These errors areto be contrasted with the effect of measurement noise orerrors in the data-set, which will be presented in futuresubsections, and thus, F χ is noiseless in the discussionbelow.Let the following equations represent a true linear sys-tem of equations and a perturbed linear system, denotedvia ˜ as: W V χ c = W F χ (B9) W ˜ V χ ˜ c = W F χ , (B10)Next, we impose the ∞ or the uniform norm which sat-isfies consistency when applied to matrices and vectorssuch that (cid:107) Ax (cid:107) ≤ (cid:107) A (cid:107) (cid:107) x (cid:107) , where A and x represent ar-bitrary matrices and vectors. Let E be an error matrix such that ˜ V χ = V χ + EW F χ = W ( V χ + E )˜ c = W V χ c (B11)= ⇒ ( W V χ + W E )˜ c = W V χ c (B12)= ⇒ ( I + ( W V χ ) † W E )˜ c = c (B13)= ⇒ ˜ c − c = ( W V χ ) † W E )˜ c (B14)= ⇒ (cid:107) ˜ c − c (cid:107) ≤ (cid:13)(cid:13) ( W V χ ) † W E ) (cid:13)(cid:13) (cid:107) ˜ c (cid:107) (B15)In the above simplification, we assume ( W V χ ) † is a leftpseudoinverse; and we apply consistency of the ∞ -norm.If the quantity (cid:13)(cid:13) ( W V χ ) † W E ) (cid:13)(cid:13) <
1, then the error in theestimated interpolation coefficients ˜ c relative to the true(unknown) coefficients c can be written as the followinginequality [33]: (cid:107) ˜ c − c (cid:107)(cid:107) c (cid:107) ≤ (cid:13)(cid:13) ( W V χ ) † W E ) (cid:13)(cid:13) (1 − (cid:107) ( W V χ ) † W E ) (cid:107) ) (B16)To interpret this inequality, we explicitly introduce thecondition number by expanding and re-writing the term0 (cid:13)(cid:13) ( W V χ ) † W E ) (cid:13)(cid:13) in terms of the condition number η ( V ): (cid:13)(cid:13) ( W V χ ) † W E ) (cid:13)(cid:13) ≤ (cid:13)(cid:13) ( W V χ ) † (cid:13)(cid:13) (cid:107) W E (cid:107) (B17)= (cid:13)(cid:13) ( W V χ ) † (cid:13)(cid:13) (cid:13)(cid:13)(cid:13) W ˜ V χ − W V χ (cid:13)(cid:13)(cid:13) (B18)= η ( V ) (cid:13)(cid:13)(cid:13) W ˜ V χ − W V χ (cid:13)(cid:13)(cid:13) (cid:107) W V χ (cid:107) (B19) η ( V ) := (cid:13)(cid:13) ( W V χ ) † (cid:13)(cid:13) (cid:107) W V χ (cid:107) (B20)In this derivation, we apply the consistency of the uni-form norm and we substitute the expression for E . Thenorm of the matrix W V χ is added as a dummy variable,yielding the final expression for the condition number.When (cid:13)(cid:13) ( W V χ ) † W E ) (cid:13)(cid:13) <
1, we can write the followingapproximate inequality that links the condition numberof V χ as the magnification of true errors when any errorsin the matrix V χ are passed on to the solution for c : (cid:107) ˜ c − c (cid:107)(cid:107) c (cid:107) ≤ η ( V ) (cid:13)(cid:13)(cid:13) W ˜ V χ − W V χ (cid:13)(cid:13)(cid:13) (cid:107) W V χ (cid:107) (B21)(B22)It is well known that the optimal condition number isone, and in general, η ( V ) ≥ ∞ -norm ofthe interpolant can be seen by re-arranging terms to ob-tain (cid:107)L ( f ) (cid:107) ≤ η ( V ) (cid:13)(cid:13) ρ ( x ) T (cid:13)(cid:13) (cid:107) c (cid:107) . Note that the conditionnumber does not depend on f but only the structure ofthe linear system which we are trying to solve. Hence,for a general point set χ and choice of polynomial basisfor P κ , η ( V ) can be large.Under an optimal choice of basis, however, the condi-tion number satisfies η ( V ) = 1 i.e. errors are not mag-nified in the inversion process [30]. This optimal basisfor the least squares problem is in fact an ortho-normalpolynomial basis with respect to a discrete inner productover χ [30]: (cid:104) p i , p j (cid:105) := L (cid:88) k =1 w k p i ( x k ) p j ( x k ) (B23)= δ i,j , { x k } Lk =1 ∈ χ, { p i } Ni =1 (B24)The geometry of Ω influences whether an optimal basiscan be found for a given application.
3. Lebesgue constant
The ∞ -norm of the interpolant (cid:107)L ( f ) (cid:107) is typically usedin approximation error analysis - in particular, to an-swer the question “how good is L ( f ) as an approxima-tion to f ?” To address this question, one seeks the so-called Lebesgue constant, Λ L , to provide the lowest upperbound on the norm of the interpolant:Λ L := min { λ ≥ (cid:107)L ( f ) (cid:107) ≤ λ (cid:107) f (cid:107) , ∀ f ∈ C (Ω) } (B25) Here, the constant Λ L is independent of the form of thefunction f or the construction of L ( f ). Hence, the be-haviour of Λ L enables a definition of optimal across dif-ferent interpolation strategies or different geometries.In addition to providing a bound on the ∞ -norm ofthe interpolant, we can also use Λ L to provide an upperbound on approximation error. To do this, we considerthree points in the space of continuous polynomials andapply the triangle inequality [30]. Since P κ ⊂ C (Ω), wemay consider f, L ( f ) and the best approximating poly-nomial p ∗ as points in the space of continuous functionson Ω. One applies the triangle inequality to obtain: (cid:107) f − L ( f ) (cid:107)≤(cid:107) f − p ∗ (cid:107) + (cid:107) p ∗ − L ( f ) (cid:107) , ∀ f ∈ C (Ω)(B26)Here, the best approximating polynomial is optimal withrespect to the ∞ -norm (cid:107) f − p ∗ (cid:107) , and thus, it should bedistingushed from the least-squares criterion used to con-struct the interpolant L ( f ).Using L ( p ∗ ) := p ∗ and the linearity of the interpolant,one substitutes: (cid:107) p ∗ − L ( f ) (cid:107) = (cid:107) L ( p ∗ − f ) (cid:107) (B27) ≤ Λ L (cid:107) p ∗ − f (cid:107) (B28)= ⇒ (cid:107) f − L ( f ) (cid:107) ≤ (1 + Λ L ) (cid:107) f − p ∗ (cid:107) (B29)The last step combines the two inequality relations andenables the interpretation of the Lebesgue constant as ameasure of how much worse an interpolant performs withrespect to some best approximating polynomial p ∗ . Appendix C: Sources of Error1. Measurement errors in true f In this section, we analyse the effect of errors in thefunctional values at the point-set χ . Let ˜ f denote theperturbed function for the true f , with the error in func-tional values expressed as: (cid:15) ( x ) := f ( x ) − ˜ f ( x ) , f, ˜ f ∈ C (Ω); x ∈ χ (C1)The error in the interpolant is thus linear by the linearityof the operator L : L ( f ) − L ( ˜ f ) := L ( f − ˜ f ) (C2)= L ( (cid:15) ) (C3)Letting F χ and ˜ F χ denote the corresponding vectorswhere each element is computed using f and ˜ f respec-tively, we obtain an expression for the error in the inter-polant: L ( (cid:15) )( x ) := ρ ( x ) T ( W V χ ) † W ( F χ − ˜ F χ ) , x ∈ Ω (C4)If we assume further that (cid:15) ∈ C (Ω), we can use the defi-nition of the Lebesgue constant: (cid:107)L ( (cid:15) ) (cid:107) ≤ Λ L (cid:13)(cid:13)(cid:13) f − ˜ f (cid:13)(cid:13)(cid:13) (C5)1. Under the specific condition that (cid:107)L ( f ) (cid:107) ≥ (cid:107) f (cid:107) , theinequality above can be cast to establish (cid:107)L ( f ) (cid:107) as thecondition number for passing on errors in functional val-ues to the interpolant: (cid:107)L ( (cid:15) ) (cid:107)(cid:107)L ( f ) (cid:107) ≤ Λ L (cid:13)(cid:13)(cid:13) f − ˜ f (cid:13)(cid:13)(cid:13) (cid:107) f (cid:107) (C6)In 1D, this condition is met for the Lagrange polynomials,leading to the straightforward interpretation that inter-polation strategies with minimal Λ L reduce sensitivityof the resulting interpolant to errors in functional val-ues (see [31] for 1D case). The condition (cid:107)L ( f ) (cid:107) ≥ (cid:107) f (cid:107) is generally not true for an arbitrary polynomial basisfor an interpolation strategy in 2D. For the specific casewhere χ represents the Padua point set on a unit square,one can derive the optimality relations discussed in sub-sequent sections.
2. Perturbed point-set χ In this section, we assess the impact of a perturbedpoint-set χ . For a chosen χ and polynomial basis for P κ ,finite size effects mean that accessing exact x ∈ χ duringpractical applications is impossible.There is very limited research on this subject (see [34]for the case of 1D continuous function on an equispacedgrid with a fixed perturbation); and as-yet no analysis forbi-variate interpolation problems on specific geometries(e.g. the unit square).In lieu of a formal derivation of error bounds, we con-sider the case that each point x in the selected point-set χ is weakly perturbed by a fixed displacement (cid:15) → (cid:15) := ˜ x − x, ∀ x ∈ χ (C7)The introduction of these errors now result in two non-linear perturbations: firstly, in the elements of V χ → V χ,(cid:15) and similarly, F χ → F χ,(cid:15) , where the subscript χ,(cid:15) denotesthat the matrix or vector elements are being computedusing the perturbed points ˜ x .For weak noise, the size of the perturbations dependon the derivatives of the function f and the chosen poly-nomial basis { p j } for the space P κ .A first order Taylor expansion recasts the effect of theseperturbations on the functional values F χ,(cid:15) as approxi-mately linear: f (˜ x ) := f ( x ) + (cid:15)f (cid:48) ( x ) + O ( (cid:15) ) (C8)= ⇒ F χ,(cid:15) := F χ + (cid:15)F (cid:48) + O ( (cid:15) ) (C9) F (cid:48) := (cid:2) f (cid:48) ( x ) . . . f (cid:48) ( x L ) (cid:3) (C10)Thus, to first order, F χ,(cid:15) will manifest as measurementerrors of the previous section. Higher orders may be con-sidered depending on the specific form of f . The strengthof the errors will depend on the derivatives of f . Similarly, first order expansions of the basis polynomi-als lead to an expression for the perturbed Vandermondematrix: p j (˜ x i ) := p j ( x i ) + (cid:15)p (cid:48) j ( x i ) + O ( (cid:15) ) (C11)= ⇒ V χ,(cid:15) := V χ + (cid:15)V (cid:48) + O ( (cid:15) ) (C12) V (cid:48) := [ p (cid:48) j ( x i )] (C13) j = 1 , . . . , N (C14) i = 1 , . . . , L (C15)We now observe that the quantity V (cid:48) will have a first col-umn consisting of all zeros; and the remaining elementsform the basis for the polynomial space P κ − . Hence V (cid:48) is singular, but a pseudo-inverse for the quantity W V (cid:48) exists. We substitute these new matrices into the inter-polant formula to get: L ( (cid:15) )( x ) := ρ ( x ) T ( W V χ + (cid:15)W V (cid:48) ) † W ( F χ + (cid:15)F (cid:48) ) . x ∈ Ω(C16)We use the properties of the Moore-Penrose pseudoin-verse for any pair of m × n matrices A and B [35],( A + B ) † = 12 (cid:2) I n I n (cid:3) (cid:20) A BB A (cid:21) † (cid:20) I m I m (cid:21) (C17)to re-write the pseudo inverse of the first and second order L × N Vandermonde matrices: V ( (cid:15) ) := ( W V χ + (cid:15)W V (cid:48) ) † (C18)= 12 (cid:2) I N I N (cid:3) (cid:20) W V χ (cid:15)W V (cid:48) (cid:15)W V (cid:48) W V χ (cid:21) † (cid:20) I L I L (cid:21) (C19)= 12 (cid:15) (cid:2) I N I N (cid:3) (cid:20) W V χ /(cid:15) W V (cid:48) W V (cid:48)
W V χ /(cid:15) (cid:21) † (cid:20) I L I L (cid:21) , (cid:15) (cid:54) = 0(C20)where the last line is obtained by using the basic propertythat the pseudoinverse of any non-zero scalar multiple ofmatrix A is the reciprocal multiple of the pseudoinverse A † satisfying ( (cid:15)A ) † = (cid:15) − A † , for some (cid:15) (cid:54) = 0.Under the weak perturbation approximation, and tofirst order in (cid:15) , we see that the combined effect of a per-turbed χ manifests as both errors in V ( (cid:15) ) and in func-tional evaluation: L ( (cid:15) )( x ) =12 ρ ( x ) T (cid:2) I N I N (cid:3) (cid:20) W V χ /(cid:15) W V (cid:48) W V (cid:48)
W V χ /(cid:15) (cid:21) † (cid:20) I L I L (cid:21) W ( F χ (cid:15) + F (cid:48) ) , (C21) (cid:15) (cid:54) = 0In this form, it appears that the total error will be me-diated by the condition number η ( V ( (cid:15) )) and additionalmagnification of order | (cid:15) | in functional evaluation in theweak error limit | (cid:15) | →
0. For f approximately slowly-varying or constant, we can ignore errors in functional2evaluation caused by the perturbed grid and set F (cid:48) ≈ (cid:15) inside the pseudo-inverseterm and the functional evaluation term do not easilycancel out in general and only cancel if the off-diagonalfirst order terms additionally satisfy V (cid:48) ≈
0. It remainsan open question to see if a perturbed grid can be recastas functional evaluation errors under some other consid-erations.
Appendix D: Bivariate polynomial interpolation atthe Padua points1. Optimality on the unit square
For the unit square, Ω := [ − , , it was discoveredthat the product Chebshev polynomials evaluated at theso-called Padua point-set χ p , and with the correct theo-retically derived weights W p form an optimal basis withrespect to the discrete inner product. Here, the sub-script p denotes Padua-based interpolation strategies.Further, the interpolation problem presented in earliersections turns out to be unisolvent, enabling a uniqueleast squares solution for the interpolants c p and the ma-trix W p V χ p turns out to have a unity condition number[30].For the set of 2D problems mappable to the unitsquare, bi-variate Lagrange interpolation at the Paduapoints also yields the slowest-growing error bound forapproximation error as Λ L p ∼ O (log ( κ )) where κ de-notes the order of the Padua points to interpolate f atmost degree κ [13–15]. In particular, for some constant a ( f, k ) that depends on f ∈ C k (Ω) and its k continuousderivatives [16]: (cid:107) f − L ( f ) (cid:107) ≤ (1 + Λ L ) (cid:107) f − p ∗ (cid:107) ≤ a ( f, k ) log ( κ ) κ k (D1)Here, κ represents both the space of polynomials of de-gree at most κ as well as the order of the Padua pointset; k denotes the number of continuous derivatives and a is a constant that depends on both the function f and k .
2. Lagrange-Padua interpolation
Padua points can be generated in 3 equivalent ways:(a) via the use of a generating curve approach [13]; theideal theory approach [14] and through the merger oftwo Chebychev-like grids [15]. There are four familiesof Padua points and we focus on the first family for theequations below.Let x ∈ χ be a point in the set of Padua points χ oforder κ over the unit square Ω := [ − , . The numberof Padua points depends on the order, κ , as: | χ κ | := ( κ + 2)( κ + 1)2 , κ > γ κ ( t ) on the unit square, Ω.The intersection of this generating curve with itself, theedges of the square or its vertices yields the Padua pointset. The generating curve can be defined as: γ κ ( t ) := ( − cos(( κ + 1) t ) , − cos( κt )) (D3)On this curve, the Padua points are equispaced along t and indexed by j, m as: t ( j,m ) := jκ + m ( κ + 1) κ ( κ + 1) π, j, m ≥ , j + m ≤ κ (D4)The set of Padua points are classified as interior, bound-ary or vertex points. Two vertex points occur at (1 , − κ , ( − κ +1 ) and edge points occur on theboundary of the square. The curve of Eq. (D3) is consis-tent with [15].The so called cubature weight w x of the point x de-pends on its classification: w x := 1 κ ( κ + 1) · / , x vertex1 , x boundary2 , x interior (D5)The cubature weights above agree with [15].The formulae above enable interpolation for any true f using Lagrange polynomial interpolation of degree atmost κ . The interpolation formula is written as [13, 14]: L ( f )( x (cid:48) ) := (cid:88) x ∈ P κ f ( x ) w x ( K κ ( x, x (cid:48) ) − T κ ( x [0]) T κ ( x (cid:48) [0]))(D6)= (cid:88) x ∈ P κ f ( x ) l ( x, x (cid:48) ) (D7) l ( x, x (cid:48) ) := K ∗ ( x, x (cid:48) ) K ∗ ( x, x ) (D8) K ∗ ( x, y ) := K κ ( x, y ) − T κ ( x [0]) T κ ( y [0]) (D9) K κ ( x, y ) := κ (cid:88) j j (cid:88) i T i ( x [0]) T j − i ( x [1]) T i ( y [0]) T j − i ( y [1])(D10) w x := 1 K ∗ ( x, x ) (D11)where ( x [0] , x [1]) are the coordinates of x in 2D, T κ ( · )are the Chebychev polynomials of order κ , K κ ( x, y ) is areproducing kernel for the space of bivariate polynomialson the unit square with degree at most κ [14, 15]. TheLagrange basis polynomials form an orthonormal basisover the Padua points, having the property l ( x , x ) =1 ⇐⇒ x = x and zero otherwise, for any two Paduapoints x , x ∈ χ κ .We supplement the geometric picture of Padua pointswith an alternative formulation that enables rapid calcu-lation. In this alternative picture, the Padua points are3a subset of a grid of Chebyshev-Gauss-Lobatto points, χ κ ⊂ C κ +1 × C κ +2 . The interpolant can be written interms of a matrix of coefficients , C ( · ), and a rectangularChebyshev matrix T ( · ): L ( f )( X ) := (cid:0) ( T ( X )) t C ( f ) T ( X ) (cid:1) t (D12) T ( S ) := ˆ T ( s ) . . . ˆ T ( s m )... . . . ...ˆ T κ ( s ) . . . ˆ T κ ( s m ) (D13) S =[ s , . . . , s m ] , ˆ T y ( · ) := √ T y ( · ) (D14)Here, X := X × X is a discrete Cartesian grid, X i isa vector of the i -th coordinate of all test points for theinterpolation of the function f . The notation t denotesa matrix transpose. The Chebyshev matrix T ( X i ) hasdimensions κ × dims( X i ) and has elements given by thescaled Chebyshev polynomials √ T y ( · ) of order y .The matrix of coefficients, C ( f ) is computed as essen-tially the left upper triangular component of a ( κ + 1) × ( κ + 1) square matrix C ( f ): C ( f ) := T ( C κ +1 ) G ( f )( T ( C κ +2 )) t (D15)There is one modification: the last element of the firstcolumn of C ( f ) is multiplied by a factor of one-half. Theconstruction of the remaining matrices will be discussedbelow.The rectangular Chebyshev matrices are now definedby a vectorised set of Chebyshev-Gauss-Lobatto pointsof order κ + 1: C κ +1 := { z κj = − cos(( j − π/κ ) , j = 1 , . . . , κ + 1 } (D16)In the above, the negative sign is required so that thePadua points, generating curve and the Chebyshev gridsyield the same point-set.Lastly, the ( κ + 1) × ( κ + 2) matrix G ( f ) incorporatesthe effect of the values of f evaluated on Padua pointset, as well as the Padua curbature weights. Entries ofthis matrix are non-zero only if the index of the matrixelement coincides with a Padua point: G ( f ) := ( g r,s ) (D17)= (cid:40) w x f ( x ) , if x = ( z κr , z κ +1 s ) ∈ χ κ , if ( z κr , z κ +1 s ) ∈ C κ +1 × C κ +2 \ χ κ (D18)The entries of G ( f ) which coincide with the Padua pointscan be quickly discovered by selecting only every otherpoint in flattened ‘meshgrid’ of ( κ + 1) × ( κ + 2) entries.The flattened mask is subsequently reshaped into a 2Dmatrix for odd values of κ as ( κ + 1) × ( κ + 2). For evenvalues of κ , the mask is reshaped to ( κ + 2) × ( κ + 1),followed by a matrix transpose. If this mask is appliedto the meshgrid of the Chebyshev-Gauss-Lobatto points C κ +1 × C κ +2 , the Padua points of the generating curveapproach in the previous section are recovered. Compu-tationally, it is easier to construct G ( f ) by generating andmerging the two Chebyshev grids, C κ +1 , C κ +2 to obtainthe Padua points.The connection with the geometric interpretation iseasier to see if we re-write the Padua points with slightlymodified index notation [14]: x [0] := cos kπκ , ≤ k ≤ κ (D19) x [1] := (cid:40) cos (2 j − π ( κ +1) , k evencos (2 j − π ( κ +1) , k odd j = 1 , . . . (cid:100) κ/ (cid:101) + 1(D20)The formulae above result in duplicate Padua points thatcan be removed by inspection; direct comparison with thedefinition of the Chebyshev-Gauss-Lobatto points con-firms that χ κ ⊂ C κ +1 × C κ +2 .Two additional modifications to our Python code-baseare required: firstly, the L ( f )( X ) matrix is flipped fromleft-to-right (corresponding to the use of the left uppertriangular matrix of coefficients). Secondly, L ( f )( X ) isglobally divided by a factor of four to compensate for thescaling factors in the rectangular Chebyshev matrix T ( · ).Both of these modifications yield the final result in thecorrect orientation consistent with [15].4 [1] John Preskill, “Quantum computing in the NISQ era andbeyond,” Quantum , 79 (2018).[2] Julian Kelly, Peter O’Malley, Matthew Neeley, Hart-mut Neven, and John M Martinis, “Physical qubit cal-ibration on a directed acyclic graph,” arXiv preprintarXiv:1803.03226 (2018).[3] Aaron D Tranter, Harry J Slatyer, Michael R Hush, An-thony C Leung, Jesse L Everett, Karun V Paul, PierreVernaz-Gris, Ping Koy Lam, Ben C Buchler, and Geoff TCampbell, “Multiparameter optimisation of a magneto-optical trap using deep learning,” Nature communica-tions , 1–8 (2018).[4] DT Lennon, H Moon, LC Camenzind, Liuqi Yu,DM Zumb¨uhl, GAD Briggs, MA Osborne, EA Laird,and N Ares, “Efficiently measuring a quantum device us-ing machine learning,” npj Quantum Information , 1–8(2019).[5] Paul B Wigley, Patrick J Everitt, Anton van den Hen-gel, John W Bastian, Mahasen A Sooriyabandara, Gor-don D McDonald, Kyle S Hardman, Ciaron D Quinlivan,P Manju, Carlos CN Kuhn, et al. , “Fast machine-learningonline optimization of ultra-cold-atom experiments,” Sci-entific reports , 1–6 (2016).[6] Sergey K. Tolpygo, Vladimir Bolkhovsky, Terence J.Weir, Leonard M. Johnson, Mark A. Gouker, andWilliam D. Oliver, “Fabrication process and proper-ties of fully-planarized deep-submicron Nb/Al AlO x /NbJosephson junctions for vlsi circuits,” IEEE Transactionson Applied Superconductivity , 112 (2015).[7] Daniel Koch, Brett Martin, Saahil Patel, Laura Wessing,and Paul M. Alsing, “Benchmarking qubit quality andcritical subroutines on IBM’s 20 qubit device,” (2020),arXiv:2003.01009 [quant-ph].[8] Frank Arute, Kunal Arya, Ryan Babbush, Dave Bacon,Joseph C. Bardin, Rami Barends, Rupak Biswas, Ser-gio Boixo, Fernando G. S. L. Brandao, David A. Buell,and et al., “Quantum supremacy using a programmablesuperconducting processor,” Nature , 505510 (2019).[9] Swarnadeep Majumder, Leonardo Andreta de Castro,and Kenneth R Brown, “Real-time calibration with spec-tator qubits,” npj Quantum Information , 1–9 (2020).[10] Shyam Shankar, Michael Hatridge, Zaki Leghtas,KM Sliwa, Aniruth Narla, Uri Vool, Steven M Girvin,Luigi Frunzio, Mazyar Mirrahimi, and Michel H De-voret, “Autonomously stabilized entanglement betweentwo superconducting quantum bits,” Nature , 419–422 (2013).[11] Olivier Landon-Cardinal, Luke C. G. Govia, andAashish A. Clerk, “Quantitative tomography for contin-uous variable quantum systems,” Phys. Rev. Lett. ,090501 (2018).[12] Marco Caliari, Stefano De Marchi, and Marco Vianello,“Bivariate polynomial interpolation on the square at newnodal sets,” Applied Mathematics and Computation ,261–274 (2005).[13] Len Bos, Marco Caliari, Stefano De Marchi, MarcoVianello, and Yuan Xu, “Bivariate Lagrange interpo-lation at the Padua points: the generating curve ap-proach,” Journal of Approximation Theory , 15–25(2006).[14] Len Bos, Stefano De Marchi, Marco Vianello, and Yuan Xu, “Bivariate Lagrange interpolation at the Paduapoints: the ideal theory approach,” Numerische Mathe-matik , 43–57 (2007).[15] Marco Caliari, Stefano De Marchi, Alvise Sommariva,and Marco Vianello, “Padua2DM: fast interpolation andcubature at the Padua points in Matlab/Octave,” Nu-merical Algorithms , 45–60 (2011).[16] Marco Caliari, Stefano De Marchi, and Marco Vianello,“Bivariate Lagrange interpolation at the Padua points:Computational aspects,” Journal of Computational andApplied Mathematics , 284–292 (2008).[17] Riddhi S. Gupta and Michael J. Biercuk, “Conver-gence analysis for autonomous adaptive learning appliedto quantum architectures,” (2019), arXiv:1911.05752[quant-ph].[18] Riddhi Swaroop Gupta, Alistair R. Milne, Claire L. Ed-munds, Cornelius Hempel, and Michael J. Biercuk, “Au-tonomous adaptive noise characterization in quantumcomputers,” (2019), arXiv:1904.07225 [quant-ph].[19] Zuzana Majdisova and Vaclav Skala, “Radial basis func-tion approximations: comparison and applications,” Ap-plied Mathematical Modelling , 728–743 (2017).[20] Michael L Stein, Interpolation of spatial data: sometheory for kriging (Springer Science & Business Media,2012).[21] Richard Franke,
A critical comparison of some meth-ods for interpolation of scattered data , Tech. Rep. (NavalPostgraduate School Monterey CA, 1979).[22] Ioan Dor´e Landau, Rogelio Lozano, Mohammed M’Saad,and Alireza Karimi,
Adaptive control: algorithms, analy-sis and applications (Springer Science & Business Media,2011).[23] Sebastian Thrun, Wolfram Burgard, and Dieter Fox,
Probabilistic robotics , Vol. 1 (MIT press Cambridge,2000).[24] B. Stanley, personal communication, in-person meetingwith WaveSense CTO at Greentown Labs, Somerville,MA (2019).[25] Steffen Michalek, Mirko Wagner, and Jens Timmer, “Anew approximate likelihood estimator for ARMA-filteredhidden Markov models,” IEEE Transactions on SignalProcessing , 1537–1547 (2000).[26] Carl Edward Rasmussen and Hannes Nickisch, “Gaussianprocesses for machine learning (GPML) toolbox,” Jour-nal of machine learning research , 3011–3015 (2010).[27] Jaakko Luttinen and Alexander Ilin, “Efficient gaussianprocess inference for short-scale spatio-temporal model-ing,” in Artificial Intelligence and Statistics (2012) pp.741–750.[28] Jouni Hartikainen, Jaakko Riihim¨aki, and Simo S¨arkk¨a,“Sparse spatio-temporal gaussian processes with generallikelihoods,” in
International Conference on ArtificialNeural Networks (Springer, 2011) pp. 193–200.[29] Kang Li and Yun Fu, “ARMA-HMM: A new approachfor early recognition of human activity,” in
Proceedings ofthe 21st International Conference on Pattern Recognition(ICPR2012) (IEEE, 2012) pp. 1779–1782.[30] Marc Van Barel and Matthias Humet, “Good point setsand corresponding weights for bivariate discrete leastsquares approximation,” Dolomites Research Notes onApproximation , 37–50 (2015). [31] Bayram Ali Ibrahimoglu, “Lebesgue functions andLebesgue constants in polynomial interpolation,” Jour-nal of Inequalities and Applications , 93 (2016).[32] Yurii Kolomoitsev, Tetiana Lomako, and J¨urgen Prestin,“On Lp-error of bivariate polynomial interpolation on thesquare,” Journal of Approximation Theory , 13–35(2018).[33] Gilbert W Stewart, Afternotes on numerical analysis , Vol. 49 (Siam, 1996).[34] Anthony P Austin and Lloyd N Trefethen, “Trigonomet-ric interpolation and quadrature in perturbed points,”SIAM Journal on Numerical Analysis , 2113–2122(2017).[35] Yongge Tian, “The Moore-Penrose inverse for sums ofmatrices under rank additivity conditions,” Linear andMultilinear Algebra53